EPA/600/R-09/121
September 2009
Americamysis bahia Stochastic
Matrix Population Model for
Laboratory Populations:
Technical Documentation
By
Glen B. Thursby
Atlantic Ecology Division
National Health and Environmental Effects Research Laboratory
Office of Research and Development
US Environmental Protection Agency
Narragansett, Rhode Island 02882
-------
Notice
The research described in this report has been funded wholly by the U.S. Environmental
Protection Agency. This report is contribution number AED-09-048 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 or commercial products does not
constitute endorsement or a recommendation for use.
Abstract
The population model described here is a stochastic, density-independent matrix model for
integrating the effects of toxicants on survival and reproduction of the marine invertebrate,
Americamysis bahia. The model was constructed using Microsoft® Excel 2003. The focus of the
model is on laboratory populations because neither biological variability of field populations nor
variability associated with time-varying toxicant concentrations, as might be expected in field
situations, are considered. The model employs several aspects of traditional population viability
analysis (PVA), establishing a dose-response relationship between exposure concentrations and
estimates of expected minimum population size. This documentation describes the model, and
also provides the justification for various default parameters for use when the ideal toxicity data
set may not be available.
Keywords: Americamysis bahia, PVA, matrix model, risk curves
Ackn owledgments
The work described builds upon earlier independent efforts by Drs. Anne Kuhn and Jason Grear
for the US EPA's Atlantic Ecology Division. The verification Visual Basic for Applications code
in the various macros and the verification that the Excel spreadsheets do what they were intended
to do was provided by the Software Engineering and Specialized Scientific Support (SES3)
contract under RSGIS Task Order 1502 to Raytheon and SAVCS3 Task Order 1503 to Computer
Sciences Corporation. However, any errors that remain are the sole responsibility of the author.
Denise Champlin and Dr. Anne Kuhn helped by providing access to mysid data. These data, with
the exception of the pesticide chlorpyrifos, were from AED's historical test data. Chlorpyrifos
data were provided by OPP. Some of the mysid reproduction data in section three were provided
by Dr. Sandy Raimondo from EPA's Gulf Ecology Division. Finally, Drs. Diane Nacci, James
Heltshe, Giancarlo Cicchetti and Wayne Munns reviewed the final manuscript and model.
11
-------
Preface
EPA's Office of Pollution Prevention and Toxic Substances (OPPTS) needs population modeling
methods to support assessments of the risks from toxic chemicals to non-target populations of
aquatic and terrestrial species ("wildlife"). ORD's Safe Pesticides/Safe Products (SP2) Multi-
Year Plan (EPA 2006) and NHEERL's Implementation Plan for SP2 Research (EPA 2005)
recommend population modeling approaches to address this need. The modeling effort described
herein is specifically targeted towards OPPTS' Office of Pesticide Programs (OPP) in support of
their effort to integrate population models into risk assessments associated with the pesticide
registration process. The work described in this plan contributes to the accomplishment of SP2
Long-term Goal 2's annual performance goal: provide methods for characterizing population-
level risks of toxic chemicals to aquatic life and wildlife.
The model is for the marine invertebrate, Americamysis bahia, a common toxicity test organism,
and one for which OPPTS has specific test guidelines for acute and chronic tests (OPPTS
850.1035 and OPPTS 850.1050). Currently, acute and chronic toxicity data are treated separately
within the risk assessment process for pesticide registration. The model allows the integration of
these kinds of data into a single risk assessment endpoint. Because the model uses laboratory
derived data for both its demographic parameters and the toxicity dose-response, the target
population is essentially a laboratory population.
There are several components to this modeling work. One such component is the mathematical
construction of the model which includes the selection of data used to derive the demographic
parameters for the base model. These data were taken from the control treatments of historical
tests conducted at EPA's Atlantic Ecology Division. An additional component of the modeling
effort provides default input values for needed toxicity information that may be missing for a
given pesticide. OPPTS has recently updated its toxicity data requirements for the registration of
pesticides. However, data that meet these new requirements may not be sufficient for all of the
model's needs. Therefore, if population modeling is to be incorporated reasonably soon into the
risk assessment process for pesticides, it needs to be done with existing data requirements.
The model is constructed using Microsoft® Excel 2003—filename MYSID 13x13 MATRIX—
LAB EXPOSURE 2009 Version 1.XLS. A copy of this file is embedded below. Excel's Solver
add-in may need to be activated for new users of the model.
Note
This PDF file has an attached Excel file at the end of this paragraph. You can either double left
click or single right click the icon. When you double left click the attached Excel file should
open. When you single right click on the icon a pop up menu should appear. From this pop up
menu choose either Open File or Save Embedded File to Disk,
in
-------
Table of Contents
Page
Abstract /'/'
Acknowledgements /'/'
Preface /'/'/'
Section 1: Introduction 1
Section 2: Model Structure—Summary of Demographic Parameters for
Americamysis bahia for the Matrix Population Model 5
Section 3: Minimum Toxicity Data Requirements for Americamysis bahia
Matrix Population Model—Justification for Default Parameters 11
Survival Dose-Response Curve 11
LC50 Versus Time 13
Reproduction Dose-Response Curve 18
Section 4: Sources of Stochasticity 21
Variability in Survival 21
Variability in Reproduction 26
Variability in Sensitivity to Toxicants 28
Environmental Stochasticity 30
Section 5: Model Validation of Default Values 31
Section 6: Description of Excel Version of Mysid Matrix Population Model 34
Sheet: Input Parameters 34
Sheet: Control Survival Curve 38
Sheet: Parameter Calculations 38
Sheet: Mysid Matrix 13x13 41
Sheet: Lambda Summary 41
Sheet: Quasi-Extinction 43
Sheet: Simpsons Rule 43
Sheet: PVA Solver 45
Sheet: Data Entry Error Check 45
Sheets Summary 45
Appendix A. Equations used to model dose-response data 46
Appendix B. Calculation of kinetic coefficient for time-to-effect curves 48
Appendix C. List of compounds used for Figure 8 in the main text 51
Appendix D. Graphs from chronic tests 52
Appendix E. Supplemental information for Weibull survivorship curve 56
Appendix F. Biological Significance 58
Appendix G. List of VBA macro code 60
References 69
IV
-------
Section 1: Introduction
The U.S. Environmental Protection Agency registers pesticides for intended use in the United
States under the authority of the Federal Insecticide, Fungicide, and Rodenticide Act (FIFRA). A
pesticide can be licensed for use if it will not cause "unreasonable adverse effects on the
environment," Section 2, Paragraph bb of FIFRA defines this as "any unreasonable risk to man
or the environment" or "a human dietary risk" from food residues. The focus of this document is
on the environmental portion of the above. Risk assessments for the registration of pesticides
have been based historically on the toxicological effects on organism-level vital rates such as
survival and reproduction. The intention, however, is to evaluate risks to populations, not to
individual organisms. In this context, a population-level endpoint may be more desirable. This
documentation shows a different approach for using data from traditional organism-level toxicity
tests to estimate population-level risks. The specific application being evaluated is risk to
populations of the marine invertebrate Americamysis bahia (Molenock, 1969). This is
accomplished through a population model which is a mathematical representation using vital
rates (e.g., survival and reproduction) to describe population dynamics, such as population
growth rate (A). Population risks from pesticides (and other chemicals and stressors) are
projected using a population model integrating both acute and chronic toxicity data from
standard organism-level toxicity tests into predictions of long term viability of a population. This
model demonstrates the potential for routine use of population level endpoints in pesticide risk
assessments, and paves the way for more detail modeling efforts that will be aimed at more field-
like situations. The document not only describes the inner workings of the model, but also
provides the necessary documentation of default model input parameters when such parameters
are not readily available from the toxicity data.
Population models are a convenient mathematical method to integrate multiple sources of effects
(e.g., effects on survival and reproduction) into a single endpoint. There are, however, additional
values to using population models instead of, or in addition to, standard toxicity test dose-
response data. When examining only organism-level data "it cannot be assumed that large effects
on the vital rates translate into large contributions to the effects on A" (Caswell 1989). In other
words, a compound cannot be assumed to cause the same degree of effect on population growth
rate as it does on survival or reproduction independently. Caswell (1996) further states:
The take-home message of every example that has been published so far is that it
is not safe to assume that the most obvious effect of a toxicant on the vital rates is
the source of that toxicant's effect on A. In each case there are large vital rate
effects that make only trivial contributions to effects on A and other, much smaller
vital rate effects whose contributions are much larger.
There are many examples of this in the literature. In a study on the effects of cadmium on
nematodes, reproduction was more sensitive to cadmium than was juvenile period (duration as
juvenile); however, change in juvenile period had a greater effect on population fitness than did
effects on reproduction (Kammenga et al. 1996). In addition, Gleason and Nacci (2001)
demonstrated that "populations of species that have different life history strategies can respond
differently to a stressor producing responses of similar type and magnitude at the individual
level." In summary, not only will population modeling endpoints integrate stress to multiple vital
-------
rates, but they may actually result in a different conclusion relative to which vital rate has the
most significant effect on the vitality of the population.
The particular mathematical model chosen for this product is a stochastic, density independent
matrix population model created with Microsoft® Excel 2003. The model evaluates exposure to a
constant concentration of a toxicant on laboratory populations. The normal operation of the
model includes a series of such concentrations, resulting in a dose-response for the probability of
decline in a population. The endpoint for the model derives from population viability analysis
(PVA) that is commonly used in conservation biology, and requires estimates of variability in
population growth rate. The PVA is based on the concept of quasi-extinction probabilities
introduced by Ginzburg et al. (1982). As explained by Morris and Doak (2002):
The basic goal of PVA is to predict the future with some reasonable degree of
assurance. We can do this much better if we don't try to predict when the very last
desert tortoise in Las Vegas County may die, but, rather, when the population will
reach a small enough number that many additional genetic and ecological
problems will further threaten it, making it perilously at risk.
The population threshold at which a group of organisms is at risk is referred to as the quasi-
extinction threshold. In this context, risk is described as the probability of a population declining
below this threshold.
PVA analysis requires incorporating estimates of vital rates and their variability into projections
of quasi-extinction. The source of the variability for the mysid model comes from the variance in
mysid vital rates under unstressed conditions (i.e., laboratory survival and reproduction) and the
uncertainty associated with responses to toxicants. This version of the model does not
incorporate biological variability within field populations or variability associated with time-
varying concentrations of a toxicant. It is expected that the output from this model will augment
the current use of mysid chronic test data in the FIFRA risk assessment process. Although the
model was developed primarily for use within EPA's Environmental Fate and Effects Division,
Office of Pesticides Programs' (OPP) risk assessments for registration of pesticides, its use is not
restricted to that program.
The model described herein provides a mechanism for expressing the toxicity of single
compounds as a series of risk curves. The risk curves are built using PVA. These also are
referred to as "risk of decline curves" (Ak9akaya 2000). The concept was originally described by
Ginzburg et al. (1982), but was not specifically referred to as risk curves. Simply put, a risk
curve is a plot of the probability of decline relative to the magnitude of the decline. The
probability is greatest for a small threshold of decline and least for a very large threshold. Risk
curves can be expressed either as a threshold population size—estimating the probability of
declining below that threshold (e.g., McCarthy 1996, McCarthy and Broome 2000, and Taylor et
al. 2003)—or as percentage decline (e.g., Inchausti and Weimerskirch 2004, and Schtickzelle et
al. 2005). The latter is the way the risk curves are expressed in this mysid model. With either
method, however, a judgment has to be made to determine the threshold or percentage decline of
concern. This is not always an easy judgment to make. As Morris and Doak (2002) point out,
"the size of the quasi-extinction threshold is often the subject of argument." This issue is
eliminated if risk is evaluated by comparing the area between different risk curves, each
representing a different environmental scenario (Burgman et al. 1993, see Figure 3.17). In the
case of the mysid model presented here these different scenarios represent different toxicant
-------
concentrations. McCarthy and Thompson (2001) took this one step further showing that the area
to the left1 of a single quasi-extinction curve was equivalent to the expected minimum population
size. This expected minimum is expressed as either an absolute number of organisms or as a
percentage of the current population size—depending on how the threshold is expressed. Thus,
using a series of toxicant concentrations in a model will allow the construction of a dose-
response curve with expected minimum population size as the response.
Americamysis bahia has been a standard toxicity test species for several decades, and the utility
of population level toxicity endpoints for this species was demonstrated almost from the
beginning (Gentile et al. 1982, 1983). In spite of this, the application of population models in the
interpretation of the effects of toxicants on this species has not been routinely used. Recently
there have been some attempts at developing matrix population models for A. bahia. Kuhn et al.
(2000, 2001) developed a deterministic age class matrix model with a daily time step using
demographic data from standard 28-day chronic early life-cycle toxicity tests. The matrix model
covered approximately one-third of the mysid's expected 90 d life history—using 28 separate
daily age classes. The authors demonstrated that the traditional chronic test results were likely to
be protective of mysid population growth. A different daily time step matrix model for A. bahia
was created by Raimondo et al. (2005a,b). These authors also used demographic data from
traditional chronic toxicity tests. However, they developed a stage class model with seven
different stages based on age. The oldest stage in the model was represented by late breeders (25
to 29 d). The authors of this stage class model also differed from Kuhn et al. (2000, 2001) by
incorporating stochastic simulations so that quasi-extinction estimations could be made, but
complete risk curves were not created. Both of the above matrix models were successfully used
to demonstrate the utility of population models for mysids, but neither fit the exact needs of the
current model.
The desired matrix model covers the entire life-cycle of A. bahia—approximately 90 days. The
models by Kuhn et al. (2000, 2001) and Raimondo et al. (2005a,b) could have been modified to
cover this time frame. The former was not used for this modification because daily age classes
would have made the model size too large for ease of use. The latter could have been modified
by adding one or a few additional stages to cover the remainder of the life cycle. This was not
done because it would still result in a model with a daily time step. In the next phase in the
development, the mysid model presented here will become part of a periodic series covering an
entire year. This means that a separate matrix model will be needed for each time step within the
year. Covering an entire year with a series of matrices will allow OPP to evaluate annual
pesticide concentration time series for their risks to populations. Currently OPP only uses a few
single concentration values (e.g., peak values and 21 d averages) from their annual estimated
exposure scenarios in the derivation of their risk quotients. Using a series of matrices will allow
the entire concentration time series to be used, and will allow the incorporation of time-varying
pesticide concentrations into the modeled effects. This in turn will allow exposure scenarios that
are more similar to what is expected in nature. To avoid using 365 separate matrices associated
with a daily time step, a model using a longer time step will be needed. A monthly time step is
convenient, because a year is easily divided into thirteen 4-week "months". This time step,
1 To the right of the risk curve applies in situations where the degree of population change is expressed as
percentage decline—as with this mysid model. To the left applies for those curves expressed as population
thresholds.
-------
however, will be too long because the month would represent an exceedingly large portion of the
mysid life cycle. A weekly time step lets the life history be divided into 13 separate weekly age
classes. Although this will require 52 separate matrices when the annual version of the mysid
model is developed, it is a reasonable compromise, and certainly more manageable than 365
daily matrices.
The rest of this document is divided into five sections. The next section describes the model
structure, and includes a summary of the baseline data used to develop the biological parameters
for the model. The third section details the minimum toxicity data requirements needed for the
model. Data submitted to OPP must meet certain minimum acceptability requirements. These
data, however, may not provide all of the information needed for a successful use of the
population model. Therefore, this third section also provides the technical documentation for the
selection of default values when desired data are absent. A summary of the data used to
determine the variability in model parameters (these are the source of stochasticity needed for
PVA analysis) is provided in the fourth section. The fifth section is an example showing a
comparison among model runs using data from toxicity tests conducted with the insecticide
endosulfan. This data set is sufficient to provide all of the necessary parameters for running the
model. The model was also run assuming that the data were not sufficient and that default values
for survival kinetics and reproduction dose-response curves were needed. This fifth section
provides an example of the difference in the endpoint relative to the concentration-response
curve when using a complete data set versus using default values for survival kinetics or
reproduction dose-responses. The sixth, and final, section contains a detailed description of each
Excel spreadsheet used in this version of the model.
-------
Section 2: Model Structure—Summary of Demographic Parameters for
Americamysis bahia for the Matrix Population Model
The model is an age class, density independent, stochastic matrix population model. The life
history of Americamysis bahia is divided into 13 one-week age classes, covering the expected
duration of the life cycle. The model time step (projection interval) is also one week. This means
that at each step of a model simulation an individual either dies or survives and grows into the
next age class. Density independent means that the survival and reproduction rates are not
influenced by the size of the population. Although density dependence has been incorporated
into mysid models for various reasons (e.g., Raimondo et al. 2005a), the relationship between
density and vital rates has not been quantified. The reference given for the existence of density
dependence is Lussier et al. (1988). That reference, however, only has two minor statements
about density dependence. The first:
Population density within a culture affects reproduction. A culture that is too
densely populated will cease reproduction.
Followed latter by:
A steady decline of a culture is more difficult to ascertain and may result from a
combination of factors such as inadequate nutrition, over-harvesting, chronic
disease, contaminant organisms, high mysid density, high male-to-female sex
ratio, or sublethal toxicity from algae, Artemia cysts or dilution water.
There is insufficient information to accurately portray density dependence; therefore, it was not
incorporated into this model. Finally, the need for stochasticity is because the model output relies
on an application population viability analysis (PVA) within risk curves. This means that the
output is based on the probability or the likelihood that a population will decline. The
underpinning of this type of analysis requires an estimate of variability in population growth rate.
This variability is described in Section 4: Sources of Stochasticity.
The matrix population model developed for A. bahia is shown in Figure 1. A Pt value is the
probability the ith age class survives during a time step (one week), therefore moving into age
class /'+!. The Ft values represent the fecundity for the ith age classes, which is the number of
young produced and surviving the time interval. The calculations for these parameters are shown
below in Equations 1 and 2. The grey area in the upper left of the matrix represents the portion of
the life cycle for which measured values are available. The matrix age classes, however, cover
the entire life history of A. bahia. This is needed because the sensitivity to a toxicant is in part a
function of the duration of exposure. Representing the entire life history in the matrix allows
tracking the potential effect of a toxicant on individuals of different ages no matter how old they
are.
-------
Age
1
2
3
4
5
6
7
8
9
10
11
12
13
10
11
12 13
Pi
0
0
0
0
0
0
0
0
0
0
0
2
0
p*
0
0
0
0
0
0
0
0
0
0
c
r3
0
0
P3
0
0
0
0
0
0
0
0
0
c
r4
0
0
0
P4
0
0
0
0
0
0
0
0
c
rs
0
0
0
0
P5
0
0
0
0
0
0
0
c
re
0
0
0
0
0
Pe
0
0
0
0
0
0
7
0
0
0
0
0
0
PT
0
0
0
0
0
8
0
0
0
0
0
0
0
PS
0
0
0
0
9
0
0
0
0
0
0
0
0
P9
0
0
0
10
0
0
0
0
0
0
0
0
0
Pio
0
0
0
0
0
0
0
0
0
0
0
0
Pn
0
0
0
0
0
0
0
0
0
0
0
0
Pl2
0
0
0
0
0
0
0
0
0
0
0
0
Figure 1. The 13x13 matrix population model used for /Amer/camys/s bahia. P, and F, are defined with
Equations 1 and 2. The grey area represents the extent of the standard 28d chronic test.
The survival and reproduction (maternity) data used to calculate baseline (or control) survival
probabilities and fecundities are shown in Table 1. These data were summarized from ten, 28-
day chronic tests conducted at the Atlantic Ecology Division (AED) between the years of 1987
and 1993. To be included in the dataset for estimating the demographic parameters, tests had to
meet reproductive and temperature requirements of the current OPPTS test guidelines2. This
meant that tests from the late 1970s and early 1980s were not included because they were
conducted at a lower temperature, which delays reproduction.
2 Ecological Effects Test Guidelines. OPPTS 850.1350. Mysid Chronic Toxicity Test. Office of Prevention,
Pesticides and Toxic Substances.
-------
Table 1. List of chemicals used in the evaluation of control demographic parameters for the mysid
Americamysis bahia. The six digit number following the chemical name refers to the AED experiment
number which takes the form of YYMM##. Note, there is no experiment number for triethylene
glycol/acetone.
Compound
Acenaphthene 890703
Acenaphthene 890805
Aniline 871 001
Carbofuran 921001
Dichlorvos 900101
Flouranthene 910702
Propoxur 900601
Pyrene 920701
Triethylene glycol/acetone
Thallium 880210
Mean
SD
7(l)a
0.967
0.933
0.950
1.000
0.917
0.883
0.867
1.000
1.000
0.950
0.947
0.0023
1(2)
0.933
0.833
0.933
1.000
0.900
0.867
0.800
1.000
0.933
0.883
0.908
0.0043
7(3)
0.700
0.633
0.783
0.900
0.867
0.833
0.783
0.850
0.833
0.717
0.790
0.0070
7(4)
0.567
0.450
0.683
0.767
0.750
0.700
0.683
0.767
0.733
0.617
0.672
0.0103
b
7M3
6.086
4.127
3.195
3.100
3.875
0.625
2.938
1.769
1.533
0.600
2.785
2.913
ff?4
5.417
7.825
6.500
4.091
6.479
2.545
5.813
3.423
2.114
3.643
4.785
3.589
a l(x) = the probability of survival from birth to age x.
b mx = maternity rate. Average total number of young per female at age x.
Americamysis bahia is a birth-flow (Caswell 2001) population, meaning births occur
continuously over the projection interval. Because of this, the age of individuals within a given
age class is not precisely known. For example, those individuals within age class 1 are between 0
and 1 week old; those within age class 2 are between 1 and 2 weeks old; and so on. The
probability of an individual in a given age class surviving to the next age class (P,) depends on
the age of that individual; therefore, there needs to be some sort of average for each age class3.
This probability is estimated using Equation 2.24 from Caswell (2001):
M±ffi±2>
Likewise the number of births in a projection interval has to be averaged in some way because
individuals born near the beginning of the interval have to survive the entire interval before
being counted and those born near the end of the interval have a higher probability of surviving
until the end of that interval. Caswell' s equations 2.34 and 2.35 are used to estimate the average
fecundity for a given age class
Fi = 7(0. 5) - '>>^ Equation 2
3 See Section 2.4.1 of Caswell (2001) for a more complete explanation.
4 Note: Equation 2.34 from Caswell has 2 in the denominator instead of 4. This is because Caswell's maternity rate
represents the number of females per female. Since here the maternity rate is the total number of young, this change
converts the results to number of females per female.
-------
Where 7(0.5) is approximately equal to (7(0)+/(l))/2. Note: when expressed in terms of
probability 7(0) is 1.0.
The survival rates and fecundity values for the first 4 age classes can be calculated directly from
the available data. These parameters for the remainder of the age classes, however, have to be
estimated. One option for these latter survivals is to assume that the survival rate for weeks 4
through 12 are the same and use a survival probability of 0.0 for week 13 so that life ends.
Survival, however, is a function of time and most, if not all species fit the model of either a Type
I, Type II or Type III survivorship curve. Knowing this a survivorship curve can be fit to the
averages from the data in Table 1 and estimates made of the survival probabilities for weeks 4
through 13. Survivorship curves can typically be fit with a Weibull function. There are several
versions of the Weibull equation; the one used herein is:
5(0 = e- Equation 3
This is Equation 6.2.3 from Lee and Wang (2003), where t is time and S(t) is the probability of
survival at time t. The constants k\ and &2 are shape parameters. If k\ is greater than 1.0, then the
survivorship curve is Type I, if equal to 1.0, the curve is a Type II, and if less than 1.0, a Type
III. The longer the life span of the organism, the smaller the value of k^. Figure 2 shows a
Weibull fit to the average survival values from all 10 tests (Table 1). The shape parameters are k\
= 2.045 and ki = 0.162. Appendix E gives a more complete explanation of how these parameters
were derived, as well as the general relationship between the two shape parameters.
We have no maternity rate data beyond that for week 4. Maternity rates for weeks 5 though 13,
therefore, are assumed to be the same as week 4. Note, however, that because the equation for
fecundity includes a survival rate the fecundity rates decline as the mysids age. Table 2 shows
the numbers used to set the baseline parameters for the 13x13 matrix. These parameters are those
used to calculate the stable age distribution5 that is in turn used as the default t-zero population in
the model simulations. The end user, however, can substitute any values in the t-zero vector.
5 Calculated using PopTools version 2.7.5 available on the internet. URL http://www.cse.cisiro.au/poptools.
-------
Estimated Survivorship Curve Using Mysid Control Averages
1.00
0.90
0.10
0.00
-O—Type I
• Mysid control data
0 1 2 3 4 5 6 7 8 9 10 11 12 13
Time (weeks)
Figure 2. Weibull function fit to the first 4 weeks of survival data for Americamysis bahia.
-------
Table 2. Summary of life table parameters for the matrix population model used for Americamysis bahia. The
probability of surviving from birth to a given age is l(i). The probability of surviving from a given age to the next
age is P/. The maternity rate for a given age is m/. All ages from week 4 on have the same maternity rate.
Finally, F, is the fertility of age class ;. Note that a week 14 is added in order to correctly calculate P, for week
13 (we assumed at the maximum age of 13 weeks that 1% of the initial population would still be alive—see
Appendix E for a more complete explanation).
Age / (weeks)
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
l(i)
1.000
0.976
0.905
0.795
0.661
0.521
0.388
0.273
0.182
0.114
0.068
0.038
0.020
0.010
0.005
Pt
-
0.952
0.904
0.856
0.812
0.769
0.727
0.688
0.651
0.615
0.582
0.547
0.517
0.500
m,
-
0
0
2.785
4.785
4.785
4.785
4.785
4.785
4.785
4.785
4.785
4.785
4.785
Ft
-
0
0.622
1.700
2.141
2.091
2.041
1.995
1.951
1.909
1.870
1.829
1.793
1.182
10
-------
Section 3: Minimum Toxicity Data Requirements for Americamysis bahia
Matrix Population Model—Justification for Default Parameters
The use of the Americamysis bahia population model for assessing population level effects from
toxicants assumes that the minimum data available from the acute test will be a 96 h acute LC50
(along with its 95% confidence interval) and an LOEC (Lowest Observable Effect
Concentration) for reproduction from a 28 d chronic6. If these data are not available, then the
model cannot be used. No other data besides those provided in the OPPTS test guidelines are
required. Besides these data, the following additional information is desirable, but defaults are
available and explained below: survival dose-response curve; LC50 vs time data; and
reproduction dose-response curve.
Survival Dose-Response Curve
Equation 4 is used to model survival dose-response curves. Information about this equation and
other forms of it are shown in Appendix A. Ideally, there would be enough data to fit this curve
and get an estimate of the slope. If enough data are not available, then the exposure response
uses a default slope of 3.3. This corresponds to the default used by the Office of Pesticide
Programs (OPP) for the probit, 4.5 (Urban and Cook 1986). Both Equation 4 and the probit
method (Equation 5) yield essentially the same dose-response curve. Figure 3 shows the
relationship between the probit slope and the slope of the logistic equation for a range of dose-
response curves. To achieve the data set eighteen probit curves were created using slopes ranging
from 1.5 to 10 with an LC50 of 2.5 for each data set. A best fit for these data to the logistic
equation using least squares regression (adjusting only the slope) was then made.
Survival Probability = ;— Equation 4
^ / \ slope l
I cone
LC50
Z = (X — jit) • probit _ slope Equation 5
where Z is the standard score (or the z-score); Xis the logio of a concentration; and n = logic of
the LC50. The probability of survival is the cumulative probability of the standard normal
distribution up to the value Z . This is easily calculated with the Excel function NORMSDIST(Z).
Based on Office of Prevention, Pesticides and Toxic Substances (OPPTS) test guidance numbers 850.1035 and
850.1350 for acute and chronic tests with mysids, respectively, these data should be readily available if acceptable
tests were performed.
7 Probits are actually Z+5 since the early user of this technique did not like to deal with negative numbers (see Bliss
1934, 1935, Fmney 1971, Hubert 1980).
11
-------
Probit vs Logistic
3.0 -i
7.0 -
6.0 -
2. s.o H
_o
U)
g 4.0 ^
|S
5 3.0 H
2.0 -
1.0 -
0.0
y = 0.737x
R2 = 1.000
Probit default: 4.5
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
Probit Slope
Figure 3. Relationship between the probit slope and the slope (shape factor) for the logistic equation.
Figure 4 shows a sample set of dose-response curves using an LC50 of 2.5 and OPP's default
probit slope of -4.58 and the corresponding logistic slope of -3.3 (calculated from the equation in
Figure 3). The dose-response curves are essentially identical
' Note: the negative of the default slope was used so that survival (rather than mortality) was calculated.
12
-------
LC50 = 2.5, probit slope = -4.5
1.00 -
0.80 -
> 0.60 H
3
w
0.40 -
0.20 -
0.00
o o
Concentration
Figure 4. Default survival curves using probit default OPP slope (-4.5) and the logistic equivalent (-3.3).
LC50 Versus Time
The mysid model describe herein has a time-step of 7 days with 13 different weekly cohorts.
Therefore, an estimate is needed of how the LC50 changes with increasing exposure time. Figure
5 shows an idealized relationship between the LC50 for a given toxicant and the duration of
exposure. Equation 6 is a mathematical expression often used to model this relationship. An
LC50 at any t can be calculated, if the LC50 with infinite exposure9 is known (approached as an
asymptote) and the kinetic parameter, k.
LC50t =
Equation 6
Both the LC50 at infinity and the kinetic parameter can be estimated if there are enough data
from the acute test (e.g., 24, 48, 72 and 96 h LC50s), the chronic test (7, 14, 21 and 28 d), or a
combination of both. However, the likely scenario is that there will not be enough survival data
to directly use Equation 6. But a default value can be derived and used to estimate the LC50s for
various exposures using historical data comparing LC50s for 28 d with those for 4 d. This is
explained in detail below and in Appendix B.
' Also called the asymptotic LC50 or incipient LC50. The latter appears to have been introduced in Sprague (1969).
13
-------
LC50 Kinetics
8
o
14 21 28 35 42 49 56
Exposure Duration (da)
63 70 77 84 91
Figure 5. Idealized relationship of lethal concentration to exposure duration. Curve was created using
Equation 6 with the asymptotic LC50 (LC50«,) set at 1.0 and the kinetic parameter, k set at 0.3045
Figure 6 shows three different hypothetical curves relating LC50 to duration of exposure using
Equation 6. For convenience the LC50 is set to 1.0 on day 1 for all three curves. The curves
represent asymptotic LC50s reached in 5d, lOd and 50d. Two time durations have been marked,
4d and 28d. It is clear from the figure that as the time to reach an asymptotic LC50 increases, the
difference between the LC50 after 4 days of exposure and that after 28 days get larger. Figure 7
shows this relationship in another manner. Here the x-axis is the time to reach the asymptotic
LC5010. The y-axis is the ratio of theLC50s for days 4 and 28.
Figure 8 shows a cumulative distribution for LC50 ratios from toxicity tests conducted with 24
compounds using Americamysis bahia. The actual data used are listed in Appendix C. The 95th
percentile for the ratio of LC50s on days 4 and 28 is 5. This is used to determine the default
number of days to reach an asymptotic LC50, as shown in Figure 9. The horizontal line labeled
"A" is drawn at the above ratio of 5. A vertical line labeled "B" is drawn where "A" crosses the
curve for the ratio of days 4 and 28. This second line crosses the x-axis at 100 days. Thus 100
days becomes the default time to reach an asymptotic LC50. Using Equation B3 (Appendix B),
this yields a kinetic parameter (K) for Equation 6 of 0.03045—which becomes the default kinetic
value.
10 This is actually set as the time to reach within 5% of the true asymptotic LC50. Because the true value is never
reached an approximate value was selected that was close enough to this true value (see Appendix B for further
explanation).
14
-------
ro 0.40
8 0.20 -
0.10 -
0.00
100
Figure 6. Three different curves representing the idealized relationship between LC50 and duration of
exposure. In the upper curve the asymptotic LC50 is reached in 5d. In the middle curve, 10d. In the lower
curve, 50d. The vertical lines showtime at days 4 and 28.
7.0 n
6.0 -
5.0 -
O
O
4.0 -
10 100
Time to get to 5% of LC50 infinity (da)
1000
Figure 7. Relationship between the ratio of LC50s on days 4 and 28 to time to reach the asymptotic LC50
(represented by a value within 5% of the true asymptote). Data from the 3 curves in Figure 4 are indicated
with solid circles.
15
-------
Ratio of 96h & 28d LCSOs
ro
>
'>
D
CO
o
1O
o
0
0?
o
'c
o
£
o
3
15 -,
14 -
13 -
12 -
11 -
10 -
9 -
8 -
7 -
6 -
4 -
3 -
2 -
1 -
n
u
C
, .
I csdmium f\
• Non-Metals u
O Metals
ammonia
<5.14
•_ .
.o-°
0
-••••O«OOO
• ••00*
• • •
) 10 20 30 40 50 60 70 80 90 100
% Rank
Figure 8. Cumulative distribution of survival acute/chronic ratios for mysids. Data are from tests listed in
various EPA water quality criteria, and are shown in Appendix C. The horizontal dashed line is the 95th
percentile.
o
^
o
O
10 100
Time to get to 5% of LC50 infinity (da)
1000
Figure 9. Relationships between LC50 ratios and time to reach asymptotic LC50. The horizontal line
labeled "A" is drawn at the ratio of 5. A vertical line labeled "B" is drawn where "A" crosses the curve.
16
-------
There is an alternate method to estimate the kinetic parameter for Equation 6. If LC50 values are
available for two different periods of time (in days), then Equation B4 (Appendix B) can be used.
You solve iteratively for the value of k that yields the correct ratio of the two LC50 values.
Knowing the kinetics of how the LC50 changes with time and assuming the "slope" in Equation
4 remains constant, a series of dose-response curves can be created for each age class (e.g.
Figure 10). Therefore, the cumulative survival for any single exposure concentration can be
calculated (e.g. Figure 11). This in turn allows the calculation of the probability of survival due
to toxicity for each age class to the next by a simple ratio of two consecutive cumulative survival
rates. The use of these values is explained in the Section 6 describing the Excel spreadsheets,
Sheet: Parameter Calculations.
Sample Dose-Response Curves
2
Q.
"55
CO
Increasing Exposure
,' Duration from 1 to 5 wks
6789
Concentration
10 11 12 13 14 15
Figure 10. Sample set of dose-response curves for the first five weeks of exposure.
17
-------
Exposure Concentration = 2.5
.a
TO
.Q
2
Q.
CO
0 1 2 3 4 5 6 7 8 9 10 11 12 13
0.
Figure 11. Survival probability based on a single exposure concentration, but allowing for the shift in
dose-response curves as the duration increases.
Reproduction Dose-Response Curve
Table 3 shows the summary toxicity data from dose response curves for nine chemicals for
which data were available for acute and chronic tests (survival and reproduction—number of
young per available female reproductive day). The full dose-response curves for each chemical
are shown in Appendix D. There is a relationship between the EC50/LC50 ratio and the slope of
the dose-response curve for reproduction (Figure 12). There is also a relationship between the
EC75/LC50 ratio and the reproduction effect slope (Figure 12—the EC75 was calculated from
the equation for the dose-response for reproduction; Equation Al in Appendix A). Figure 13
shows the relationship between the slope for the dose-response for reproduction and a range of
ECx/LC50 ratios (with 'x' ranging from 10% to 75% in 5% effect increments. Each of the curves
in Figure 13 can be described with the same general equation:
repro _ slope = Y0+ aRb
Equation 7
where R is the ECxILCSO ratio, YO is the y-intercept, and a and b are shape factors for a particular
curve. These latter three values are calculated using Equations 8-10, where 'x' is the % of the
control response for a given concentration. Plots showing the relationship between each of these
three factors and % of the control are shown in Appendix D.
18
-------
70 = 2.20 - 0.09 • ln(x - 8.62) Equation 8
a = 16 + 0.6 • JK° 3 Equation 9
6 = 6.25-1.10-ln(jc-9.00) Equation 10
Why is the above important? In many data sets a dose-response curve for reproduction may not
be available. However, the data requirements for a mysid chronic test will at least yield an LOEC
value. The associated percent difference from the control for this value becomes the 'x' in
Equations 8-10. These in turn give the values to use in Equation 7 for the slope of the
reproduction dose-response curve. The EC50 for reproduction is then calculated using Equation
4 (for % reproduction instead of survival). (Note: the ECx used to calculate R in Equation 7 is the
LOEC associated with the % difference used in Equation 8-10.)
Table 3. Summary data from dose response curves in Appendix D.
Compound
Acenaphthene
Aniline
Carbofuran
Chlorpyrifos
Dichlorvos
Fluoranthene
Propoxur
Pyrene
Thallium
Acute LC50
460
2015
5.29
42
13.9
29.3
36
31.7
4830
Repro EC50
35
791
2.25
3.5
9.7
2.4
15
4.84
48.6
EC50/LC50
0.076
0.393
0.425
0.083
0.698
0.082
0.417
0.153
0.010
Repro Slope
2.75
5.2
4
2
10
1.5
4
2
1.5
19
-------
18 -,
16 -
o 14
-------
Section 4: Sources of Stochasticity
The mysid model is a stochastic matrix model. During a model run there are several sources of
Stochasticity that are used to adjust the survival and reproduction parameters of the projection
matrix. This includes variability in survival probability, maternity rate and sensitivity to a
toxicant. Variability also can be categorized by demographic and environmental Stochasticity, as
well as what can be called toxicological Stochasticity11. Demographic Stochasticity is variation in
vital rates (survival or reproduction) due to independent "chance variation in the actual fates of
different individuals" within a time step of interest (Morris and Doak 2002, page 22; see also
Lande et al. 2003, page 6). Environmental Stochasticity, on the other hand, "describes temporal
variation in vital rates driven by changes in the environment that are inherently erratic or
unpredictable (unlike, for example, predictable seasonal changes)" (Morris and Doak 2002, page
18; and Lande et al. 2003, page 6). The term toxicological Stochasticity describes the uncertainty
associated with the true value of a toxicity endpoint—LC50 for example. This latter uncertainty
could be further influenced by demographic or environmental factors, or both.
Variability in Survival
Table 4 shows the probability of survival for each of the first three age classes of mysids from
the control treatments for each of the ten chronic tests used. These data yield three pairs of mean
and variance values which were used to estimate the variance associated with any survival value.
This variability comes from two main sources. The first is traditional demographic variance
which can be based on the variance of the binomial distribution (Ak9akaya 2002), which is
var = />•(!- p)ln Equation 11
where p is the probability of survival and n is the number of individuals (which in the case of the
mysid chronic test is 60). Demographic variance is usually an issue only for "fairly small
populations" (Morris and Doak 2002, pages 23-25)12. This is demonstrated in Figure 14 for
population sizes ranging from 10 to 100 individuals. The equation for any of these curves is
y = ^max •4' (* ~ *2) Equation 12
Where y is the binomial variance, ymax is the maximum variance (defined as 0.25/w), and x is the
mean survival rate. This equation simplifies to
y
= (x - x2)/ n Equation 13
11 Because the model parameters are based on laboratory data the occurrence of what has been called observational
or measurement error is assumed to be zero.
12 Morris and Doak state: "As a rough guideline, we suggest that if there are at least 100 individuals in the
population for a count-based PVA, or at least 20 individuals in the most important life stages in a demographic PVA
(often the reproductive adults), it is usually safe to ignore demographic stochasticitiy."—page 25.
21
-------
The other source of variability is associated with the fact that the tests were conducted over a
period of several years. Thus even though the same test procedures were used and the tests were
conducted in the same laboratory, there are some expected differences in the water source, mysid
population, and personnel conducting the tests. This might be thought of as the "environmental
stochasticity" associated with laboratory tests. Table 4 lists the survival probabilities for the first
three age classes from the chronic tests listed in Table 1 (the fourth age class is not included
because, based on Equation 2, measured data is needed from the fifth age class, which is not
available).
Table 4. Individual P, values for each compound, calculated from the individual l(i) values presented in
Table 1. Note: the mean P, values in Table 4 do not exactly match the corresponding values from Table 2
because of Jensen's inequality (see Ruel and Ayres 1999).
Compound
Acenaphthene 890703
Acenaphthene 890805
Aniline 871 001
Carbofuran 921001
Dichlorvos 900101
Flouranthene 910702
Propoxur 900601
Pyrene 920701
Triethylene glycol/acetone
Thallium 880210
Mean
Var
Pi
0.966
0.914
0.966
1.000
0.948
0.929
0.893
1.000
0.967
0.940
0.952
0.0012
P2
0.860
0.830
0.912
0.950
0.972
0.971
0.950
0.925
0.914
0.873
0.916
0.0023
Pa
0.776
0.739
0.854
0.877
0.915
0.902
0.926
0.874
0.887
0.833
0.858
0.0037
The mean and variance pairs for the Pt values in Table 4 are plotted in Figure 15 along with the
binomial variance curves for n = 60 and n = 35. The variance in these three pairs includes both
demographic and "environmental" stochasticity. The portion of the variance associated with
demographic stochasticity could be subtracted—as has been recommended by some (Kendall
1998, Ak9akaya 2002)—because its inclusion can result in an overestimate of risks of population
decline. However, this source of variance has not been eliminated for two reasons. First, the
laboratory population size is large enough that this "extra" variance may not be significant.
Second, the environmental variance that is currently represented by the data in Table 4 is likely
to be small relative to the true variation expected in the field13. Leaving the demographic
stochasticity in the model is not likely to be a significant overestimate of the true risk.
For model simulations, the relationship between the total laboratory variance and mean survival
is described by Equation 11 with n = 35. Thirty-five was selected because it results in a curve
that reasonably includes the three measured points. Although how accurately this curve
represents the true relationship cannot be known, the variance has to be 0 when the survival rate
13 This variance also is likely to be smaller than that associated with conducting the tests in different laboratories—
recall that all of the data herein are from the same laboratory.
22
-------
is 0 (all dead) and when survival rate is 1 (all alive). This relationship is assumed to hold no
matter what the reason for a change in survival rate (e.g., age or toxicity).
0.025 -,
00
0.000
0.2
0.3
0.4 0.5 0.6
Survival Probability
0.7
0.8
0.9
Figure 14. Demonstration of the binomial variance for different population sizes. These are, from top to
bottom,10, 20, 30, 40, 50, 60, 100 and 500 individuals. The darkest line is the line for the typical mysid
chronic treatment (60).
23
-------
0.008 n
0.007 -
0.006 -
0.005 -
in 0.004 -
0.003 -
0.002 -
0.001 -
0.000
0.2
0.3 0.4 0.5 0.6 0.7
Mean Probability of Survival
0.8
0.9
Figure 15. Plot of binomial variance for n = 35 and n = 60 along with the three pairs of mean and
variance data points from Table 4.
The parameter for the binomial distribution is p (the probability of a successful event) which is
usually fixed for a given situation—for example probability of a tossed coin being heads.
However, in the modeling effort with the mysids p is assumed not to be fixed. At each time step
of the model simulation survival values for each of the age classes is randomly chosen. A
common continuous distribution for p values of a binomial distribution is the beta distribution.
And the distribution recommended by Morris and Doak (2002, page 275) for survival rates is the
beta distribution. All that is needed to establish the appropriate parameters for the beta
distribution is the mean and variance of the survival rate for a given age class. There are two
shape factors for the beta distribution called alpha and beta in order to be consistent with
®
Microsoft Excel. They are related to the mean and variance (Morris and Doak 2002, page 265;
Evans et al. 2000) by the following equations
alpha = x
x- (i-x,
-i
beta
-(,-,-(£±^-1]
L s
Equation 14
Equation 15
where x = mean survival rate ands2 = variance of the survival rate. Beta distributions for three
different mean survivals are presented in Figure 16 (variance was calculated using Equation 13).
24
-------
Beta Distribution Examples
', mean survival =0.10
mean survival = 0.90
mean survival = 0.50
-------
Variability in Reproduction
Table 5 shows additional data for the maternity rate during chronic tests using the mysid
Americamysis bahia. These data are from ORD's Gulf Ecology Division (GED—provided by
Sandy Raimondo) and are summarized by different age classes than data from AED. These four
means and variances along with the two means and variances from Table 2 are plotted in Figure
17. There are two significant conclusions from Figure 17. First, the variance increases as the
mean increases. And second, the mean and variance are approximately equal. One appropriate
distribution to use for maternity rate is the stretched beta distribution (Morris and Doak 2002,
page 281). A stretched beta distribution is similar to a "regular" beta distribution except that it
allows upper and lower limits other than 0 and 1. Equations 14 and 15 are modified slightly to
accommodate these limits by substituting the following for the mean and variance (Morris and
Doak 2002, page 281)
x- lower
M = Equation 16
upper - lower
2
g
V = Equation 17
(upper - lower}
where M and V are the "new" mean and variance used in Equations 13 and 14 to calculate the
shape factors for the beta distribution.14 Figure 18 shows several examples of stretched beta
distributions based on the mean maternity rate for age class 4 from Table 1. Because of the
relationship shown in Figure 17 the mean and variance are assumed equal. In the model
simulations this will be true not only for the controls, but also for the means adjusted due to
toxicity.
Table 5. Summary of mysid reproduction data in number of young per female of the given age class
(original data provided by Sandy Raimondo from EPA's Gulf Ecology Division). Data column headings are
the age classes that GED uses in their mysid chronic tests.
Chemical
DEF
Endosulfan
Fenthion
Methoprene
Silver Nitrate
Thiobencarb
mean
variance
1 6 days
1.44
0.50
0
0
0
0
0.32
0.34
20 days
1.36
3.67
2.375
1.122
1.13
0.67
1.72
1.23
24 days
1.71
3.81
6.82
4.222
2.90
2.07
3.59
3.44
28 days
0.00
5.64
4.22
2.1
3.75
3.67
3.23
3.79
14 Note, minimum value of 0.3 is used for both shape factors.
26
-------
Mysid Reproduction
6.0 -,
5.0 -
4.0 -
HI
u
| 3.0 -
>
2.0 -
1.0 -
0.0
y = 0.9136x
R2 = 0.843
0.0
1.0
2.0 3.0 4.0
Mean number young per female
5.0
6.0
Figure 17. Relationship between mean and variance for among test variability for mysid reproduction.
Open circles are data from AED (Table 1) and closed circles are data from GED (Table 5).
Stretched Beta Distributions
control
si
s
o
a.
34567
Maternity Rate (# young per female)
10
Figure 18. Sample of stretched beta distributions used for maternity rates for mysids. Control is for mean
of 4.785. For the other curves the mean was reduced to 0.75, 0.50 and 0.25 of the control. Variance was
assumed to be equal to the mean.
27
-------
Variability in Sensitivity to Toxicants
When assessing the risk of toxicants to population endpoints, the uncertainty associated with
toxicity data can be as great as or greater than any other source of variability (Barnthouse et al.
1990). Minimum toxicity data requirements for the mysid acute test include the calculation of
both the 96h LC50 and its 95% confidence interval. This is sufficient information to create a
distribution from which to incorporate the variability of toxicity into the stochastic mysid model.
Median lethal values (LC50s) are often assumed to have a lognormal distribution (Stephan et al.
1985). All that is needed is an estimate of the mean and standard deviation of the log (LC50s).
As a start, an estimate of the arithmetic mean and its standard deviation are already available —
reported 96h LC50 and its 95% CI, from which the standard deviation can be estimated. Because
the LC50 has a lognormal distribution the 95% CI is often asymmetrical. A typical symmetrical
95% CI is calculated by subtracting and adding the following from the mean
1.96-cr/V« Equation 18
where o is the standard deviation and n is the sample size15. Even though the 95% CI is
asymmetrical the confidence interval width is used to estimate the standard deviation by setting
two times Equation 17 equal to this confidence interval width (CW) and solving for sigma
CW = 2(1.96- a 1 4n) Equation 19
a = Cwjn 1 4 Equation 20
With an estimate of the arithmetic mean and its standard deviation, they are converted into mean
and standard deviations of the associated lognormal distribution. This is done using the following
relationships (Evans et al. 2000)
E(X}2
= \n(E(X)) -- ^ - \-LJ- Equation 2 1
O = Li 1 + Var(X) I Equation 22
11 I E(Xf)
where (i and o are the mean and standard deviation of the lognormal distribution and E(X) is the
expected value of the arithmetic mean and var(X) is its variance. Figure 19 shows three examples
of a lognormal distribution. Each example has the same arithmetic mean (10.0—units are not
relevant), only the arithmetic standard deviation is different. The sample size for Equation 20 is
12 since this is the minimum from the OPP guidance (two replicates, five treatments and a
control).
15 Sample size for the acute test is the number of replicates per treatment times the number of treatments (including
controls).
28
-------
246 8101214161820222426283032343638404244464850
LC50
246 8101214161820222426283032343638404244464850
LC50
5- 0.15 -
246 8101214161820222426283032343638404244464850
LC50
Figure 19. Sample lognormal distributions for
LC50 values. Arithmetic mean was 10.0 for
each with standard deviations of 3.0 (top), 6.0
(middle) and 12.0 (bottom). Confidence interval
widths were 3.47, 6.93 and 13.86, respectively.
29
-------
For the sake of simplicity the model couples the variability in reproduction (EC50s) directly to
that for survival (LC50s). This is done by simply maintaining a constant LC50IEC50 ratio. This
is analogous to using a fixed acute-chronic ratio (ACR) in the derivation of water quality criteria
(Stephan et al. 1985). Little, if any, data are available to describe the effect of exposure duration
on the EC50 for reproduction. Therefore, the EC50 is assumed to remain constant for each age
class.
Environmental Stochasticity
The target population for this model is a typical laboratory population of Americamysis bahia.
Therefore, this version of the model does not incorporate what might be considered traditional
environmental Stochasticity—variability associated with field populations having "good" or
"bad" weeks for survival or reproduction. This source of Stochasticity could be incorporated
easily into the model. However, as of now, there is insufficient information on the magnitude of
environmental Stochasticity for mysids.
30
-------
Section 5: Model Validation of Default Values
As stated in Section 3, the use of the Americamysis bahia population model assumes that the
minimum data available for a pesticide of interest will be the 96h acute LC50 (along with its
95% confidence interval) and an LOEC (Lowest Observable Effect Concentration) for
reproduction from a 28 d chronic—along with its response expressed a percentage of control.
Besides these data, the following additional information is desirable, but defaults are available:
survival dose-response curve (probit slope needed); LC50 vs time data; and reproduction dose-
response curve (an EC50 plus associated slope—based on Equation Al).
Data are available for the toxicity of the insecticide endosulfan to A. bahia that provide all of the
above information, including the additional desired data. Thus the model can be run for
endosulfan without needing to rely on the use of default parameters. This section compares the
model output using this complete data set with two other runs of the model where some of the
parameters were replaced with default values. Acute survival data were available from historical
records at EPA's Atlantic Ecology Division and were summarized as part of an interlaboratory
comparison (Schimmel 1981). Chronic survival and reproduction data were available from a
similar interlaboratory comparison of chronic tests (McKenney 1982). The data from "laboratory
4" were used from this later study. The acute and chronic survival data are presented in Figures
20 and 21. The reproduction dose-response data are presented in Figure 22.
100
0.0
24 h
1.0
2.0 3.0
Endosulfan (ug/L)
4.0
5.0
Figure 20. Dose-response curves from acute and chronic standard tests using Americamysis bahia.
Curves were fit to the original data using Equation 4.
31
-------
4.00 -i
3.50 -
3.00 -
_ 2.50 -
~ 2.00 -{
o
>o
O
~" 1 .50 -
1.00 -
0.50 -
0.00
10 15 20
Time (da)
25
30
Figure 21. Plot of the LCSOs from the data in Figure 20. Curve was fit using Equation 6.
100
0.0
2.0 3.0
Endosulfan (ug/L)
4.0
5.0
Figure 22. Dose response curve for reproduction (young per available female reproductive day)
expressed as percentage of the control response. Curve was fit using Equation A1.
32
-------
The 96 h LC50 for endosulfan is 1.29 (ig/L and the lower and upper 95% confidence limits are
1.00 and 1.75, respectively. The probit slope of the 96 h acute dose-response curve is 7.56.
Equation 6 was fit to the LC50 data in Figure 21 and the kinetic coefficient (k) is 0.27. Fitting
Equation Al from Appendix A to Figure 22 yields and EC50 of 0.89 (ig/L and a slope of 5.47.
TheLOEC was 1.26 (ig/L, representing a response that is 13% of the control.
The above endosulfan data were used for a model run using 1000 time steps and 10
concentrations plus a control. Two other model simulations were run using default values for
some of the data needs. In the first of these (default 1) OPP's default survival probit slope was
used (4.5); the reproduction dose-response curve was calculated using the LOEC and % of
control (EC50 =1.13 (ig/L and slope = 17.5); and a survival kinetic value calculated using a
second LC50 (2.43 (ig/L) at 2 days. The second default run (default 2) used the same information
as the first default run, except the survival kinetics were dictated by assuming the 100 day to
incipient LC50 default. The dose-response curves for all three model runs are shown in Figure
23. Note that, at least in this example, the use of default values results in more conservative
output. The concentrations causing a 30% decline in populations are 0.47, 0.32 and 0.27 (ig/L for
"all data", "default 1" and "default 2", respectively.
0.00 0.20 0.40 0.60 0.80
Endosulfan (ug/L)
1.00
1.20
1.40
Figure 23. Model output from three separate runs. The line labeled "all data" was the run within which all
of the input parameters were derived from the original toxicity tests data. The other two lines are from
runs within which some data were ignored and default values used instead. See text for explanation.
33
-------
Section 6: Description of Excel Version of Mysid Matrix Population Model
The platform for this version of the mysid model is a series of seven Microsoft Excel
spreadsheets. A limited amount of Visual Basic for Applications (VBA) code is used in various
macros (Appendix G) and serves to automate the processing of the various spreadsheets.
Technically one could run the model without this code; it just would be more tedious and time
consuming. The spreadsheets are explained in detail in the sections below.
Sheet: Input Parameters
The Input Parameter worksheet consists of three parts. The first is the area for inputs needed to
run the model (Figure 24). This is the place for user supplied toxicity data and information on
model runs to calculate the dose response curve.
% Completion of Current Run[
Run Time 2.03 Minutes
Figure 24. Area of the Input Parameter worksheet requesting model inputs.
34
-------
The cells in Figure 24 that are colored pale yellow are fields for data input. The green "RUN"
button begins the simulation run. Click on the "View Input Summary Graphs" box and the view
moves to the area that displays the survival and reproduction dose-response curves. The time and
% completion locations at the bottom are just for information purposes so that the user can
monitor the operation of a given simulation. Note that only the cells where data can be entered
can be accessed.
There are four main sections of the input screen. The first, Toxicant Information, is used to
identify the compound of interest. This information also is used to update the x-axis label in the
output section of the Input Parameter worksheet (see below).
The second, Survival, is for the data associated with survival from the standard 96h toxicity tests.
These include the LC50, the lower and upper 95% confidence limits, the probit slope16 of the
dose-response curve, and the sample size (n). The sample size, along with the confidence limits,
is used in Equation 20 to estimate the standard deviation of the LC50. The right side of the
Survival area has four "buttons". One of these needs to be selected to indicate which method is
being used to derive the kinetic parameter for the change in LC50 with time (the "&" in Equation
6). The first button provides the data needed to use Equation B4 from Appendix B to calculate
the kinetics factor based on known LC50s at two different durations of exposure. The next two
buttons provide the information needed to use Equation B3. The final button provides k
estimated directly from LC50s from a series of durations. The Survival Curve Factor is
explained below in the section entitled, Sheet: Survival Curve Graph.
The third area on the input sheet is the area for data on Reproduction effects. The LOEC and %
of the control response are needed if Equations 4 and 7 are going to be used to estimate the EC50
and slope. This information is only needed if the first button on the right is selected. If there are
enough data from the chronic test to directly calculate an EC50 and the slope of the dose-
response curve, then the second button should be selected and these data provided.
The last section on the input sheet is called Model Parameters for Multiple Exposure
Concentrations. The first cell indicates the number of weekly time steps that each exposure
concentration will be run. One hundred weeks is recommended for a "range finding" run to make
sure that the dose-response curve covers a sufficient range of concentrations. Once the range of
concentrations has been established, then the final run should be for at least 1000 weeks. At each
weekly time step the growth rate for that week (A) is calculated and stored (both as the original A,
and ln(A)). At the end of each concentration the mean and variance of ln(A) are calculated and
used in the PVA analysis (explained later). The next three cells in this area indicate how many
toxicant concentrations are used and what they are. The final number of data points (currently
fixed at 11) for the dose response will be the "# of Tox Cone" plus one for the control. The "Tox
Cone Increment" is the amount by which the exposure concentration increases (linearly) for each
set of exposures.
The second part of the Input Parameters worksheet displays the results from the most recent
model run (Figure 21). This area has one additional place for user input, "30% Decline at".
Thirty is the default, but any value can be inserted (only the number is needed, the "% Decline
16 The probit slope is entered, but the logistic slope (Equation 4) is actually used—calculated from the probit slope.
This was done because probit slope is what is traditionally supplied with test data received by OPP.
35
-------
at" is part of the cell's formatting). When a model run has been completed the concentration
calculated for the given effect is displayed. However, the percentage can be changed after the
model run and the calculation will be updated without rerunning the model. A "GREATER
THAN" or "LESS THAN" appears if the calculated value is outside the range of exposure
concentrations.
% Decline in Expected
Minimum Population Size
100
0
Population Growth Rate
(per week)
1.80
0.2
1
1.2
0.4 0.6 0.8
Endosulfan--ug/L
% Decline trendline —O— Growth Rate + minus 1 SD
0.00
1.4
30 % Decline at
0.47
ug/L
Figure 25. Example of output from a model run. The "% Decline at" box at the bottom is for user input of a
specific effect of interest. This is automatically calculated at the end of a run. However, it can be changed
without rerunning the model.
36
-------
Figure 25 displays four sets of data. The solid blue circles are the calculated "% Decline in
Expected Minimum Population" based on PVA analysis using changes in area under risk curves
(McCarthy and Thompson 2001). The solid blue trend line is the best fit of Equation A5 (a more
generic version of Equation 4). The solid red circles are the mean weekly population growth
rates. And the red plus markers are these growth rates minus one standard deviation. A 30%
decline in the expected minimum population size has been selected as a default cut off for
calculating an exposure concentration that results in a biologically significant difference. Thirty-
percent is the value used by the IUCN Red List (IUCN 2001) to categorize a species as
vulnerable to extinction, and this can be applied to local extirpation. Red List categories are
defined by consensus from the World Conservation Union. See Appendix F for more discussion
of biological significance.
There is a small area just below the "Print" and "Save" boxes where the default pathway for the
saved files to be placed. A new path may be entered into cell Q43 and this will update the
information in the macro that runs the save feature. If the pathway does not exist, then you will
be asked if you want it created.
The final area on the Input Parameter page (Figure 26) is an area to the right of the other two
areas that is used to display a summary of the dose response curves for the data that were entered
(top). This plot also shows a vertical line for the population "effect" so that it can be compared to
the original dose-response curves. The bottom plot displays the calculated LC50s for each of the
weeks associated with the 13 age classes of the model.
Dose Response
1.5
Endosulfan-ug/L
Kinetics of Toxicity
1.4-p
1.2-- -
1.0-
0.8-
10 11 12 13
Time (weeks)
Figure 26. Summary of the input dose-response curves (top) and the kinetics of the change in LC50 over
time for the current run of the model. These are based on the data entered into the Survival and
Reproduction areas of the Input Parameter section.
37
-------
Sheet: Control Survival Curve
At the request of OPP, a method to modify the baseline survival curve for mysids was provided.
This will allow a crude method to extrapolate the model from laboratory populations, with
relatively high survival during the early weeks (i.e., a Type 1 survivorship curve—see Figure
13), to a population that may more accurately represent the field situation with presumably
higher mortality during those early weeks (e.g., from predation). A survival curve factor of 1
gives the original data. A survival curve factor of 0.5 approximates a "Type 2" survivorship
curve. Factors greater than 0.5 result in "Type 1" curves, and factors less than 0.5 "Type 3"
curves. This worksheet is just for display purposes. The value for the factor is entered on the
Input Parameter worksheet.
Sheet: Parameter Calculations
The worksheet (Figure 27) has three primary areas used to calculate the Pt and Ft values for the
mysid matrix model THERE ARE NO DIRECT USER INPUT CELLS ON THIS
WORKSHEET. All of the needed inputs for this page are linked to the Input Parameter
worksheet; these are indicated by pale yellow file patterns in the appropriate cells. The three
main areas are Matrix: Survival, Matrix: Reproduction, and Toxicology. The first calculates
Pi values, the second calculates Ft values, and the third holds information related to dose-response
curves.
Matrix: Survival—There is a cell for the computer to update the current exposure concentration.
This value, along with the "slope", is used in Equation 4 to calculate the values in the Toxicity
Cumulative Survival column. Below this value are cells that contain the current 96h LC50 which
is sampled from a lognormal distribution using Excel's LOGINV function. The "&" is the kinetic
parameter for Equation 2 which calculates how the LC50 changes with time. The "LC50 inf is
the calculated incipient LC50 value—the LC50 at infinite exposure duration. This is the other
variable needed for Equation 6. It is calculated from a rearranged Equation 2 using the current 96
h LC50, "k," and t = 4 days. The time scale for the kinetic factor is a day, therefore when the
equation is used to fill in the values for the LC50 column the age in weeks is multiplied by 7.
The Weibull factors k\ and &2 are used in Equation 3 to fill out the column entitled, Control
Cumulative Survival—l(x). This column is used to calculate the Control Pi column using
Equation 1.
The column with the heading Toxicity Pi is calculated as the ratio of two consecutive values in
the Toxicity Cumulative Survival column. Product Pi values are the product of the pairs of
Control Pi and Toxicity Pi values for each age class. These values are used as the mean values
for the beta distribution. The paired variances also needed for the beta shape factors are
calculated using Equation 12. The shape factors, alpha and beta, are calculated with Equations
14 and 15 (the lower limit is set at 0.3 for both shape factors). These shape factors are in turn
used to select Random Pi values using Excels BETAINV function.
38
-------
Matrix: Survival
6,jMsu» 1.2S uij/L p(X\ m
weeks
0
1
2
3
LC5096H 2314 4
\t 0 27000 5
LC60 irf 1 ,5232 6
7
slope 5.57 8
9
Weibui factors 10
kl 2045 11
k2 0 1 623 12
13
Surmal Curve Factor 1 14
Matrix: Reproduction
Il0.fl 0504
Age 00
2
3
4
EC50 5
Based on LOEC 2 028 6
Based on Actual 1 597 7
8
fractional effect 0 79 9
10
11
12
13
Control
Cumulative
Survival -IM C
1 0000
09760
09047
0.7949
0.6614
05208
03878
0.2730
0.1816
0.1141
0.0677
00379
0 0201
0.0100
0.0047
total
maternity-
mix)
0
2 206
3791
3.791
3.791
3791
3791
3791
3791
3.791
3.791
3791
ontrol Pi
1 000
0.952
0.904
0.357
0812
0769
0727
0688
0.650
0615
0581
0.549
0518
0.490
var
0
2 2063
3.7908
3.7908
37908
3.7908
3.7903
37908
3.7908
37908
37908
3.7908
LC50
1 800
1.564
1.533
1.529
1.528
1 528
1528
1.528
1.528
1.528
I 528
1 528
1.528
1.528
lower limit
0
0.0300
0.0000
00000
00000
00000
0.0000
0.0000
0 0000
00000
0.0000
0.0000
Tci'ddty
Cumulative
Survival
I 0000
0884 1
0.7770
0.7575
0.7545
0.7540
07539
0 7539
0.7539
0.7539
0.7539
07539
07539
0.7539
0 7539
upper limit
0
6.0000
8.0000
80000
80000
80000
8.0000
8.0000
8.0000
8 0000
8.0000
8.0000
rt'jucity Pi
1 .0000
0.8841
08788
09749
09960
09994
0.9999
1.0000
1 0000
1 0000
1 .0000
1 .0000
1 .0000
I 0000
Adjusted
Mean
0.368
0.474
0474
0474
0474
0.474
0.474
0.474
0474
0.474
0.474
Product P>
1000
0841
0.794
0.835
0.809
0.768
0727
0688
0.650
0.615
0.581
0549
0518
0.490
Adjusted
Var
0.06129
0.05923
0.05923
0 05923
0.05923
0.05923
0.05923
0-05923
0 05921
0 05923
0.05923
variance
0.0000
0.0038
00047
00039
0.0044
00051
00057
0.0061
00065
00068
00070
0.0071
00071
0007I
alpha
03
52
52
.52
.52
52
52
52
52
,52
52
alpha
28,610"
27002
28400
27490
26 115
24,725
23 389
22116
20905
19755
18663
17627
16645
feta
1 77
1.69
1.69
1 69
1.69
1.69
1 69
1.69
1 69
1 69
1.69
beta
5389
6997
5.599
6.510
7685
9274
10.610
1 1 .383
13094
14.244
15336
16372
17.354
beta
upper
bound
1 77
1 69
1.69
1.69
1 69
1 69
1 69
1 69
1.69
1.69
1 69
PI
0.808
0851
0851
0705
0612
0.736
0.725
0.59?
0.644
059J
0.517
0.445
0346
random P;i*
rn(x) m(s-H)
0 2.481
2916 3.134
3 684 0.885
1 .255 3 520
5.751 5153
7003 2139
2950 2.123
3 556 0.038
0 059 4.356
7.3S2 3749
7.253 1 620
3 839 0.000
LC50
kinetic
parameters
006145 FALSE "k" based on 2 ICSQs
#Diy/OI FALSE "k" based on days entered
0.03045 FALSE "k" based on defalut 100 days
0 27000 TRUE Use "k" that is entered
Random
Fi
0.561
1.367
1.032
1-079
2464
2.066
1.146
0.812
0-998
2510
2005
0.822
Toxicology
086 0071 0367 0.606
Actual 5470 0890 0690
2"d time 9611 LC50 2"* LC50 Ratio
129 243 1.88372
FALSE Use ECSO based en LOEC
TRUE Use EC50 based on actual data
Figure 27. Screen "shot" of Calculate Parameters worksheet.
OJ
-------
The last information shown in the Matrix: Survival area is a boxed area to the right. The
TRUE/FALSE values in this box correspond to the kinetic buttons on the Input Parameter
worksheet. Whichever button is selected the appropriate cell in this boxed area is labeled
TRUE17. The corresponding "K" value is selected for the LC50 kinetics parameter mentioned
above.
Matrix: Reproduction—The 7(0.5) value is calculated from the survival area and is used in
Equation 8 for the Ft values. There are two EC50 values listed. The EC50s are based on the
randomly selected LC50 and the appropriate EC50ILC50 ratio from the Toxicology area. The one
used in the fractional effect calculation depends on which button is selected on the Input
Parameters worksheet.
The total maternity—m(x) column values within the reproduction area are multiplied by the
fractional effect calculated from the dose-response curve. The maternity rate is the total number
of young produced per female. As discussed in Section 4 under Variability in Reproduction, the
variance is assumed to be equal to the mean. The lower limit and upper limit columns are the
ranges used for the stretched beta distribution. The Adjusted Mean and Adjusted Var columns are
calculated using Equations 16 and 17. The alpha and beta shape factors are from Equations 14
and 15—as with survival these have a lower bound of 0.3. The column entitled, beta upper
bound, is included because on some occasions very large beta values are calculated causing an
error in the BETAINV function. The upper bound is the same as the value in the beta column
unless it exceeds 1000. In this case the upper bound is set to 1000. The random m(x) column are
maternity rates (total young) selected using Excel's BETAINV function. The Px*m(x+l) column
is a portion of Equation 2 (for convenience in calculation validation). Finally, the Random Ft
column is based on Equation 2 (the units are now the number of females per female).
The third section on this worksheet is used to calculate parameters associated with the toxicology
of the compound of interest. This area is subdivided into survival, reproduction and LC50 versus
time. All of the pale yellow cells are data linked to the Input Parameter worksheet, and should be
self explanatory. The last four cells in the survival row are calculated using Equations 20, 21 and
22, respectively. Within the reproduction row the LOEC/LC50 value is the ratio (R) used in
Equation 7. The values YO, a and b are from Equations 8, 9 and 10, respectively. The slope is
calculated using Equation 7. TheEC50 is calculated by rearranging Equation Al, solving for "c"
(EC50). The last cell in the row is the ratio of this calculated EC50 and the current LC50. Just
below the end of the reproduction row are cells containing data for an EC50 and its slope based
on direct calculation from the actual test data. Finally, the TRUE/FALSE to the right of the
reproduction data indicate which button was selected on the Input Parameter worksheet. The
"TRUE" indicates which of the two EC50s and slopes will be used to calculate the fractional
effect in the Matrix: Reproduction area.
The LC50 versus time area of the Toxicology section is the location where the LC50 kinetics
parameter is calculated only when LC50s from two different exposure durations are used. The
value for 2" time can be less than or greater than 4 days (the duration for the 96h LC50—which
17 Note: the error "#DIV/0!" appears when any needed data are blank on the Input Parameter page.
40
-------
will always be available). The values for 2nd time, 96h LC50 and 2nd LC50 are linked to the Input
Parameter worksheet. The Ratio value is the largest LC50 (shortest duration) divided by the
smaller LC50 (longer duration). The Calculated Ratio uses Equation B4 and iteratively adjusts
the "k" using Excel's Goal Seek command until this ratio is close to the actual ratio of the
LC50s. If the ratio of theLC50 values is 1.05 or less (i.e., they are within 5% of each other), then
the "k" defaults to 3.045 (see text for Equation B3) because the incipient LC50 is reached
quickly. If the ratio of the LC50s is greater than the ratio of t2/tl, then the default "k" is 0.001
(see text for Equation B4). Note that tl is always the earliest time. Whether the calculated "k" or
one of the defaults is used is determined in the appropriate cell within the "boxed" area to the
right of the Matrix: Survival area described above (using a formula in cell Q4).
Sheet: Mysid Matrix 13x13
This worksheet contains the actual matrix model and an area for the summary of the current
concentration's simulation data. Figure 28 shows a screen shot of the model area. The cells of the
matrix area are linked to the Pi and Fi columns just below the matrix. These two columns are
linked to the Parameter Calculations worksheet and are reproduced on this sheet just for
convenience. The area just to the right of the matrix has three columns of numbers, t, t+1 and
normalized (t+1). The "t" column is the column vector of the current relative size of each of the
age classes. The "t+1" column is the product of the matrix and the t vector—which results in the
vector of the relative size of each age class at the next time step (one week). At the bottom of
each of these two columns are their respective sums. The ratio of these sums is the growth rate
(lambda) for the current weekly time step. This growth rate (along with In (lambda)) is recorded
for each time step. The "normalized (t+1)" vector is each value from the "t+1" column divided
by the "t+1" total. At each new time step this normalized vector is copied and pasted to become
the new "f vector. Since this is a density independent model, this normalization step prevents
the population size from exceeding the capacity for Excel to display numbers, allowing
essentially an unlimited number of time steps.
The area entitled, Summary information for current concentration, is a location where the
information at the end of each exposure is displayed. This information is copied and pasted into
another area to the right of the matrix in Figure 28 (not shown). This latter area contains the
summary information for the entire run. Some of this is used to update the endpoint graph on the
Input Parameter worksheet. In addition to this summary information, there is also a location (not
shown) on this worksheet where the growth rate for each time step for the current exposure
concentration is recorded, and the mean and variance of In (lambda) are calculated.
Sheet: Lambda Summary
At the end of each exposure concentration for a given run of the model the growth rates
mentioned above, along with their associated mean and variance, are copied and pasted into the
Lambda Summary worksheet. This is done because these data are cleared from the Mysid Matrix
13x13 worksheet at the start of each new concentration. The Lambda Summary worksheet
provides a summary of the "raw" data for any give run of the model. Note: these data are
cleared at the beginning of each model run. It must be saved or printed if a record is
desired.
41
-------
Age Class
1
2
3
4
5
6
7
8
9
10
11
12
13
1 2 3 4 5 6 7 8 9 10 11 12 13
0
0.808
0
0
0
0
0
0
0
0
0
0
0
0.561
0
0.851
0
0
0
0
0
0
0
0
0
0
1.367
0
0
0.851
0
0
0
0
0
0
0
0
0
1.032
0
0
0
0.705
0
0
0
0
0
0
0
0
1.079
0
0
0
0
0.612
0
0
0
0
0
0
0
2.464
0
0
0
0
0
0.736
0
0
0
0
0
0
2.066
0
0
0
0
0
0
0.725
0
0
0
0
0
1.146
0
0
0
0
0
0
0
0.597
0
0
0
0
0.812
0
0
0
0
0
0
0
0
0.644
0
0
0
0.998
0
0
0
0
0
0
0
0
0
0.592
0
0
2.510
0
0
0
0
0
0
0
0
0
0
0.517
0
2.005
0
0
0
0
0
0
0
0
0
0
0
0.445
0.822
0
0
0
0
0
0
0
0
0
0
0
0
Age
Class
1
2
3
4
5
6
7
8
9
10
11
12
13
Survival Fecundity
Random Random
Pi Fi
0.808
0.851
0.851
0.705
0.612
0.736
0.725
0.597
0.644
0.592
0.517
0.445
0.346
0.56053
1 .36704
1 .03234
1 .07901
2.46381
2.06571
1.14632
0.81200
0.99794
2.51043
2.00490
0.82225
t-zero
0.44490
0.26139
0.14583
0.07704
0.03860
0.01832
0.00822
0.00349
0.00140
0.00053
0.00019
0.00006
0.00002
t
0.01948
0.00130
0.16606
0.00000
0.00663
0.00024
0.00031
0.00578
0.25499
0.03858
0.45169
0.02199
0.03296
1.000
lambda
In lambda
t+1 normalized (t+1)
1.69
0.02
0.00
0.14
0.00
0.00
0.00
0.00
0.00
0.16
0.02
0.23
0.01
2.290
2.290
0.828
0.7396
0.0069
0.0005
0.0617
0.0000
0.0018
0.0001
0.0001
0.0015
0.0717
0.0100
0.1020
0.0043
1.000
Summary information for current concentration
mean In var In Mean
cone lambda lambda lambda
1.25 -0.3659 0.3654 0.694
Figure 28. Screen "shot" of Mysid Matrix 13x13 worksheet. Note, the "summary" area displays the values from the end of the previous run.
-------
Sheet: Quasi-Extinction
This is a hidden worksheet. It can be viewed by selecting Sheet from Excel's Format menu. This
sheet is called Quasi-Extinction because this is the phrase often used by conservation biologists
in population viability analysis (PVA) to describe the risk of decline to a small population size,
first used by Ginzberg et al. (1982). This sheet calculates the probability density function (pdf)
and the cumulative density function (cdf) for the probability of a population falling below a
given threshold during a given length of time. The threshold population size is expressed as a
proportion of the initial population size (i.e., it can range between 0 and 1). These are based on
Equations 3.4 and 3.5 from Morris and Doak (2002). The model currently only records the cdf
•I Q
value at 30 weeks . However, both the pdf and cdf are calculated for time steps 1 through 250.
This was done in part so that mean and variance values from Morris and Doak (2002) could be
inserted in order to recreate their Figures 3.4 and 3.5 as a quality assurance step for the checking
the accuracy of the Excel version of their equations19.
Sheet: Simpsons Rule
This worksheet is named Simpsons Rule because this is the mathematical rule used for
approximating the areas under the risk curves. It also is a hidden worksheet. The equation can be
found in most any introductory calculus textbook. Even though the risk curves can account for
threshold proportions from 0 to 1 (which has a maximum area of 1.0), the model only uses 0.05
to 0.95 in 0.05 increments. Therefore, the calculated area is normalized as a percent of the
maximum area of 0.90.
In traditional PVA, a population threshold is selected for which the risk (or probability) of the
current population declining to that threshold size is calculated. This risk is a function of this
threshold (the smaller the change from the current population size, the greater the probability of
observing that change), as well as the mean growth rate and the variability in that growth rate20
(the greater the variability the greater the probability of decline). One problem with this approach
has been the selection of a population threshold. This is overcome by the use of risk curves in
which a range of population thresholds is used and the area between such curves calculated
(Burgman et al. 1993).
For quasi-extinction risk curves, this area represents the change in the expected
minimum population size. This is because quasi-extinction risk curves are
cumulative probability distributions of the minimum population size, and the area
to the left21 of the curve is equal to the mean of the distribution (i.e. the expected
18 The probability of falling to or below a particular threshold has to be within a given time duration. Ten
generations was chosen (30 weeks in the case of the mysidAmericamysis bahid).
The legends for Figure 3.4 and 3.5 in Morris and Doak (2002) both contain typos. The Equation numbers referred
to should be Equation 3.4 and Equation 3.5, respectively.
20 Actually what is used is the mean and variance of the In of the growth rates.
21 Note McCarthy's x-axis is the percentage of the initial population size. My calculations use an x-axis of % decline
from initial population size (as does Inchausti and Weimerskirch 2004, and Akfakaya 1992), therefore my areas are
those to the "right" of the curve rather than to the "left" of the curve—i.e., the x-axis is reversed. And the value is
the expected minimum as a percentage of the current value.
43
-------
minimum population size). Therefore, the area between two risk curves is
equivalent to the change in the mean (McCarthy 1996, page 48, bottom of second
column).
The area to the left (in the case of McCarthy 1996), or to the right in the mysid, case represents
the expected minimum population size. If you take the area under the curve you have the
proportional decline in this expected minimum—this is how the final results of the mysid model
are expressed (it is analogous to the difference between expressing a survival curve as % survival
or % mortality). Figure 29 shows some examples of risk curves for a control (curve on the
bottom left) and 6 different hypothetical exposure concentrations. The greater the concentration
the more the curves move towards the right, and the greater the percentage decline in the
expected minimum population size. "Expected minimum" is referred to because by definition a
risk curve is the cumulative distribution function for the minimum population size observed
within a given period of time (McCarthy and Thompson 2001). The change in the expected
minimum population size is equal to the area between two risk curves (McCarty and Thompson
2001). Therefore, the total area under a given curve represents the expected percentage decline in
the expected minimum populations. This is graphed in the summary plot on the Input Parameter
worksheet of the model.
Sample Risk Curves
10 20 30 40 50 60 70 80
% Decline from "Current" Population Size
90 100
Figure 29. Sample risk curves. Bottom curve is the lowest exposure concentration and the top curve is
the highest exposure concentration. As the magnitude of the decline increases, the probability for that
amount of decline decreases.
44
-------
Sheet: PVA Solver
The hidden worksheet named PVA Solver uses Excel's Solver function to make the best fit curve
to the dose-response curve for the percentage decline in the expected minimum population size.
This curve is plotted on the summary graph on the Input Parameter worksheet. The Solver
function selects the 50% effect concentration and slope for Equation A5 giving the best fit to the
above dose-response curve (i.e., a non-linear least squares regression). Equation A5 is
rearranged, using the above 50% effect concentration and slope, and solved for the concentration
that yields the "% Decline at" from the Input Parameter worksheet. If the calculated value is
greater than the highest concentration "tested," then a "GREATER THAN" appears. Likewise, a
"LESS THAN" appears if the value is less than the lowest concentration.
Sheet: Data Entry Error Check
This is another hidden worksheet. It is linked to the Input Parameters worksheet, and checks for
data entry errors. These included checking for missing data (blanks), and checking for non-
numbers or negative numbers. It verifies that the confidence limits and the 96 h LC50 are in the
correct numerical order. This sheet also verifies that when the "two LC50 method" is used to
calculate the kinetics of the LC50 that the LC50 for the shorter duration is greater than the LC50
for the longer duration. If errors in data entry are found, then an appropriate message is displayed
and the model run terminated.
Sheets Summary
The information in Section 6 summarizes the contents of all nine worksheets that are part of the
Excel model's file—four hidden sheets and five visible sheets. Only the Input Parameter sheet
remains open to the user. The other four visible sheets have to be selected to see their contents,
but this is not necessary to either run the model or save the information at the end of a simulation
run. Likewise the hidden sheets do not need to be seen for the model to be used. If you want to
view the contents of any of the hidden sheets, then just select the Format menu, followed by
Sheet, and then Unhide. You will be presented with a list of the hidden sheets; select the one
you wish to view. Any unhidden sheets will be hidden again during the next model run—so it is
not necessary to use the Format menu to hide the sheets.
45
-------
Appendix A: Equations used to model dose-response data.
There are a variety of sigmoid equations that have been used to model dose-response data
(including both survival and reproduction endpoints). The one that is chosen here is a fairly
simple one:
a . .
y = 7- Equation Al
Where:
y = proportional response of interest — for example survival or reproduction.
x = stressor concentration.
a = upper limit — for example 100 for percentage survival or percentage response relative
to the control; or it can be the actual control response for reproduction or growth.
b = slope (or shape factor) — describes the steepness of the dose-response curve. If b is
positive, then y decreases as x increases. If b is negative, then you get the opposite; y
increases as x increases.
c = is the concentration that results in a 50% response — the inflection point of the
sigmoid curve.
99
Because of the "a, b, and c", this is often called a three-parameter logistic equation . It
frequently is reported in other forms:
a
Equation A2
- - r / N / VK
1 + exp{6[ln(jt)- InfcjJ}
a- Xb
J\ ~T~ C-
~b Equation A3
If you apply mathematical rules for exponents and logarithms, Equation A2 simplifies to
Equation Al because:
Equation A3 is the "familiar" Hill equation used in biochemistry23. If you divide the numerator
and denominator of Equation A3 by xb (same as multiplying by 1), the equation becomes:
There are other equations that are also called logistic equations and they are probably just as useful for fitting
dose-response data as Equation Al above. However, I prefer Equation Al because it incorporates EC 5 Os (or
LCSOs)—for which data are often available.
46
-------
a
Equation A4
Equation A4 is slightly different from Equation Al because the Hill equation is most often used
in situation where y increases as x increases. Mathematically Equations Al and A4 are
equivalent; however, in Equation A4 unlike Equation Al if b is positive, then y increases as x
increases. If b is negative, then you get the opposite; y decreases as x increases. This is true
because:
The range for y in the above equations is from 0 to a. If there is a different lower limit, then
Equation Al is simply modified into a four-parameter logistic equation:
, a
y~d+ -T-^T Equation A5
Where d is the lower limit for y.
23 Mouton, JW, MN Dudley, O Cars, H Derendorf and GL Drusano. 2005. Standardization of
pharmacokinetic/pharmacodynamic (PK/PD) terminology for anti-infective drugs: an update. J. Antimicrobial
Chemotherapy 55:601-607.
47
-------
Appendix B. Calculation of kinetic coefficient for time-to-effect curves.
Traditionally changes in the median lethal effect as exposure duration increases have been the
primary concern. However, the mathematics involved can be applied to any effect (ECx). The
following is based on the assumption that the kinetics of time-to-effect can be described by:
yt = ^ Equation Bl
_L t-
Which is the same as Equation 4 in the main text, but where:
yt = a concentration that causes a specific percent effect (lethal or otherwise) at a given time t in
days,
y^ = a concentration for that same specific effect, but for a sufficiently long time such that no
additional adverse effects occur (the threshold or asymptotic concentration at "infinity"),
k = the kinetic coefficient that determines how quickly the effect concentration changes with
time.
We can rearrange Equation Bl using simple algebraic steps and solve for k.
yt yt yt yt
Ji-^l
-> -k = —± ^- Equation B2
Note: as t approaches °°, yt approaches yx , and y^/yt approaches 1.
Let's assume that when t is "close enough" to °° thatj^ = 1.05- yx . In other words, when yt is
within 5% of yx it is close enough to the asymptote for toxicological purposes24. Substituting
this value for yt in the equation above results in:
ml-1
1.05 j ln(0.0476)
k = — Equation B3
24 The equation for k can be solved using any percentage. Because the asymptotic LC50 is just that, an asymptote,
5% gets close enough.
48
-------
This means that if the approximate time in days at which the effect concentration nears the
threshold is known the kinetic coefficient for the time-to-effect curve can be calculated.
Sufficient data to estimate the number of days to the threshold effects concentration will not
likely be available. It can be estimated, however, from a limited amount of data that are likely to
be available. In fact the "&" used in Equation B3 can be estimated given a few more
mathematical relationships.
For any given kinetic coefficient yt can be calculated for pairs of time values (e.g., days 4 and
28) and the ratio of those y values will be a constant no matter what the value of ym. This is
shown below. For y at tj and ^ the following two equations result from Equation Bl:
and v, = ^i—
\-e*
The ratio of yt lyt is:
The yx 's cancel out:
The fractions rearrange, and the final ratio becomes:
Equation B4
- e~kt2
The value of this ratio is independent of y^ . If a standard set of time interval pairs for the ratio is
chosen, then the value of the ratio only depends on k. The relationship between this ratio and k
(or between the ratio and time-to-threshold since k = 3.045/0 can be examined. Figure Bl shows
the complete relationship between the ratio of effects on day 4 and day 28 and the exposure
duration needed to reach within 5% of the asymptotic LC50. Similar curves could be easily
created for any pair of durations. The main text sets the default ratio at 5, this corresponds to an
exposure duration of 100 days. The upper limit for the ratio in Equation B4 is the ratio t2/t}.
49
-------
o
u>
O
o
u>
O
7.0 -i
6.0 -
5.0 -
4.0 -
10 100
Time to get to 5% of LC50 infinity (da)
1000
Figure B1. Relationship between the time to the asymptotic LC50 and the ratio of the LC50s for 4 and 28
day exposure durations. This is the same figure as Figure 7 in the main text.
50
-------
Appendix C: List of compounds used for Figure 8 in the main text.
Chemical
Acenaphthene
Ammonia
Aniline
Arsenic III
Atrazine
Cadmium
Carbaryl
Carbofuran
Chlorpyrifos
Chromium VI
Copper
Dichlorvos
Fluoranthene
Mercury
Nickel
Phenanthrene
Phenol
Propoxur
Pyrene
Selenium IV
Silver
Thallium
1 ,2,4-Trichlroobenzene
Zinc
Chronic
(28 da)
LC50
211
>331
1547
818
232
8.04
9.7
2.0
0.020
>909
126.00
8.3
27
>2.51
117
7.2
3924
14.2
23
376
64
2050
235
229
Acute
(96 hr)
LC50
460
1700
1930
1740
1000
110
9.0
5.3
0.042
2030
181
13.9
31
3.50
508
27.1
8500
36.0
32
1500
141
4830
490
499
Ratio Comment
2.18
<5.14 EC40at260
1.25
2.13
4.32
13.68
0.92 Essentially a ratio of one
2.65
2.10 Chronic estimated: test lasted 35 da
< 2.23 NOEC at 909, EC10 at 68, no EC15 or above
1.44
1.67
1.15
< 1.39 EC40at2.11
4.36
3.76
2.17
2.54
1.39
3.99
2.22
2.36
2.08
2.18
References
Thursbyetal. 1989
Cardinl986
Thursby and Berry 1987b
Lussier, et al. 1985
Ward and Ballantine 1985
Lussier, et al. 1985; Gentile, et al. 1982
Thursby and Champlin 1991
Champlin and Poucher 1992b' acute
day 4 of chronic
OPP 1993, 1994
Lussier, et al. 1985
Lussier, et al. 1985
Thursby, et al. 1990
Champlin and Poucher 1991a,b
Lussier, et al. 1985
Lussier, et al. 1985
Kuhn and Lussier 1 987
Kuhn et al. 1987; Thursby and Berry
Thursby, et al. 1990
Champlin and Poucher 1992a; acute
day 4 of chronic
Ward, etal. 1981
McKenney 1982
Thursby and Berry 1988
Springborn 1988a,b
Lussier etal. 1985
from
1987a
from
-------
Appendix D. Graphs from chronic tests used to establish the relationship between
acute toxicity and the slope of the reproduction dose response curve. Equation 4
was fit to each set of data.
400 600 800 1000
Acenaphthene (ug/L)
1000 1500 2000
Aniline (ug/L)
Carbofuran (ug/L)
20 40 60 80
Chlorpyrifos (ng A.I.IL)
20 30
Dichlorvos (ug/L)
40 60
Fluoranthene (ug/L)
Figure D1. Dose response data from chronic tests using the mysid Americamysis bahia. Open circles are
survival data from 96h acute tests, closed circles are 28d chronic survival, and open squares are
reproduction data from 28d tests. All of the data, with the exception of chlorpyrifos, are from AED's
historical test data. Chlorpyrifos data were provided by OPP.
52
-------
60 80 100
Propoxur (ug/L)
120 140
15 20 25 30
Pyrene (ug/L)
4000 6000 8000
Thallium (ug/L)
10000 12000
Figure D1 continued.
53
-------
2.30 -i
2.20 -
b2.10
C
o
S 2.00 -
57
UJ
f 1.90
u>
£
a. 1.80
HI
u
a
•S 1.70
1.60 -
1.50
Y0 = 2.20 - 0.09*ln(x-8.62)
0 10 20 30 40 50 60 70 80 90 100
% of Control
Figure D2. Relationship between the y-intercept for Equation 7 and the % of the control response for a
given concentration.
18.6 -i
18.4 -
1 18.2 H
-------
7 n
6 -
.2 5
•s
£
_o
W
&
+J
0
2 -
1 -
b = 6.25-1.10*ln(x-9.00)
0 10 20 30 40 50 60 70 80 90
% of Control
Figure D4. Relationship between the shape factor 'b' for Equation 7 and the % of the control response for
a given concentration.
55
-------
Appendix E. Supplemental information for Weibull survivorship curve for
Americamysis bahia.
Data in Figure El represent the general relationship between the two shape parameters of the
Weibull equation (Equation 3). If the lifespan is fixed at a given length, then there is a
mathematical relationship between the two coefficients (k\ and &2). These values are easily
derived by selecting a series of k\ values25 and then, through an iterative process, selecting a k2
value for each k\ so that the survivorship reaches the desired minimum at the end of the lifespan.
This can be "zero", or any other value. The equation in Figure El represents a value of 0.01 (1%)
at the end of a lifespan of one (this can represent 1 day, 1 week, 1 month, etc.).
Survivorship Weibull Coefficients for Lifespan of "One"
1000
100
10
y = e
1.5270
1
0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00
Figure El
The equation for Figure El can be modified to allow for any lifespan duration. The most general
form of the &2 equation is:
25 kj greater than 1 is a Type I survivorship, equal to 1 is Type II, and less than 1 is Type III.
56
-------
,1.5270/fc!
Where omega is the lifespan duration.
Equation El
Equation El assumes that 1% of the population will still be alive at the end of the normal
lifespan. The equation can be easily adjusted to accommodate other probabilities of survival by
just changing the value of the constant. Table El lists several other constants.
Table E1. Constants for Equation E1 using different end of lifespan (w = 13) probabilities of survival.
Probability of Survival Constant
0.00001 2.438
0.00005 2.292
0.0001 2.220
0.0005 2.028
0.001 1.932
0.005 1.668
0.01 1.527
57
-------
Appendix F. Biological Significance
We have selected a 30% decline in the expected minimum population size as our default cut off
for calculating an exposure concentration that results in a biologically significant difference.
Thirty percent is the value used by the IUCN Red List (IUCN 2001) to categorize a species as
vulnerable to extinction, and this value can be applied to regional extirpation. Red List categories
are defined by consensus from the World Conservation Union.
Most scientists recognize that the goal in any experimental design is not to find a "statistically
significant difference". If this were the case, then with a sufficiently large number of replicates
any difference can be statistically significant. Kirk (1996) summarizes:
As John Tukey (1991) wrote, "the effects of A and B are always different—in
some decimal place—for any A and B. Thus asking 'Are the effects different?' is
foolish." Because the null hypothesis is always false, a decision to reject it simply
indicates that the research design had adequate power to detect a true state of
affairs, which may or may not be a large effect of even a useful effect.
What we really want to know is: "are the differences we observe among treatments, stations, etc.
different enough to be a cause of concern?" In other words, are they biologically significant (or
in the case of education or medical research—practically significant or clinically significant).
Biological significance is not easy to define, but we usually know it when we see it.26 The
decision of what constitutes a biologically significant difference requires judgment. And as Kirk
(1996) further wrote relative to practical significance:
...judgment inevitably involves a variety of considerations, including the
researcher's value system, societal concerns, costs and benefits, and so on.
The lure of "statistical significance" is that it gives the appearance of being totally objective,
turning a "continuum of uncertainty in to a dichotomous reject-do-not-reject decision" (Kirk
1996). I say "appearance of objectivity" because the common reliance on an alpha = 0.05 is
entirely arbitrary, used routinely only because of historical convention. Cohen (1990) wrote:
[It] is not a cliff but a convenient reference point along the possibility-probability
continuum. There is no ontological basis for dichotomous decision making in
psychological inquiry.
There have been numerous reviews of the fallacy of continued sole reliance on statistical
significance for evaluating results (e.g., Cohen 1990, Kirk 1996, Campbell 2005, Di Stefano
2005). Although there is some logic for abandoning—or at least supplementing—statistical
significance for a different evaluation of what magnitude of difference among treatments is of
value (e.g., practical significance or biological significance), there appears to be some reluctance
on the part of researchers to embrace any other evaluation other than statistical significance—
perhaps because the former requires the insertion of judgment. Kirk (1996) continues:
...focusing onp values and rejecting null hypotheses actually distracts us from our
real goals: deciding whether data support our scientific hypothesis and are
26 With apologies to Justice Potter Stewart
58
-------
practically significant or useful....It is true that an element of subjectivity is
introduced into the decision process when researchers make this kind of
judgment....It is a curious anomaly that researchers are trusted to make a variety
of complex decisions in the design and execution of an experiment, but in the
name of objectivity, they are not expected or even encouraged to decide whether
data are practically significant.
In fact an essential ingredient in the "research process is the judgment of the scientist" (Cohen
1990).
The issue of "what is 'biologically significant' is a major problem in conservation biology and
does not seem to have easy solutions" (Reed and Blaustein 1997). It is a "policy decision" and
should be selected by consensus (Reed and Blaustein 1997). We are recommending the Red List
definition of "vulnerable" as the default value to be used to define a biologically significant
difference. A 30% decline over some time interval of interest is the value used by the IUCN
Red List (IUCN 2001) to categorize a species as vulnerable to extinction. Red List categories are
defined by consensus from the World Conservation Union.
27 The Red List uses 10 years or 3 generations, which ever is longest. However, the model applies the 30% decline
over a duration of 10 generations—30 weeks in the case of mysids.
59
-------
Appendix G. List of VBA macro code
THIS FIRST SET ARE THE MACROS CONTAINED WITHIN THE MODULE
NAMED: Mysid_Model
Private Function StaticRandQ
'this is a function that allows the creation of a random number that does
'not change until the model calls for a new number—saw this in several places on web, but also in
'Walkenbach, John. 1999.Microsoft Excel 2000 Power Programming with VBA
'Published by Hungry Minds, Inc, New York, NY. 869 pp.
StaticRand = Rnd
End Function
Sub RUN_MYSID()
'this is the master macro for running the simulation
Application. ScreenUpdating = False
'This section checks for data entry errors
Sheets("Data Entry Error Check"). Visible = True
Sheets("Data Entry Error Check"). Select
Range("Il"). Select
If ActiveCell. Value > 0 Then
Date_Entry_Error_Check 'data entry error check macro
Sheets("Data Entry Error Check").Visible = False
Sheets("Input Parameters"). Select
Exit Sub
End If
Sheets("Data Entry Error Check").Visible = False
Application. ScreenUpdating = True
Sheets("Input Parameters"). Select
ActiveSheet.Unprotect
Range("I23").Select 'hides "selected" cell behind object—cosmetic only
'the "N" cells below are used to calculate time elapsed—font is white to hide these cells
Range("N3").FormulaRlCl = "=now()"
Range("N4").FormulaRlCl = "=now()"
Application. ScreenUpdating = False
Range("N3").Copy
Range("N3").PasteSpecialPaste:=xlValues '"freeze" start time
Dim NumSteps, NumConc As Integer
Dim InitConc, IncrConc
NumSteps = Range("D22") 'number of time steps in weeks
NumConc = Range("D23") 'number of treatments in addition to initial concentration
InitConc = Range("D24") 'usually zero
60
-------
IncrConc = Range("D25") 'the amount of increase for the concentration for each new exposure
ClearPrevious 'calls macro that clears the data from the previous run
LC50_kinetics_goalseek 'calls macro that finds survival kinetic parameter when 2 LC50s are
used
'LOOP FOR EXPOSURE CONCENTRATIONS
For k = 0 To NumConc
Sheets("mysid matrix 13xl3").Select
Range("Z2") = (IncrConc * k) + InitConc 'sets current exposure concentration
'set initial conditions
Range("F21 :F33").Copy 'copies the t-zero population numbers
Range("P3").PasteSpecialPaste:=xlValues
Application. CalculateFull 'new set of random numbers
Range("Q19").Copy 'starting In lambda
Range("X2").PasteSpecialPaste:=xlValues
'LOOP FOR NUMBERING TIMES TEPS
Range("W2").FormulaRlCl = "0"
For i = 1 To NumSteps
Range("W2").Offset(i, 0).Select
ActiveCell.Value = i
Next i
'LOOP FOR LN LAMBDAS
For j = 1 To NumSteps
Range("R3:R15").Copy 'copies the normalized t+1 population
Range("P3").PasteSpecial Paste: =xlValues 'paste the normalized t+1 into the "t" location
Application. CalculateFull 'new set of random numbers
Range("Q19").Copy 'copies the new lambda value
Cells(j + 2, 24).PasteSpecial Paste:=xlValues 'paste the new lambda into the data list
Nextj
'at end of cone loop, copy raw In lambda values
Columns("W:AA").Copy
Sheets("Lambda Summary").Select
Range("Al").Offset(0, (k * 5) + k).PasteSpecial Paste: =xlValues
Sheets("mysid matrix 13x13").Select
'copy summary statistics for cone
Range("B39:E39").Copy
Cells(k + 2, 30).PasteSpecial Paste: =xlValues
Range("U5") = (k + 1) / (NumConc + 1)
Sheets("Simpsons Rule"). Visible = True 'unhides the worksheet while it is needed
Sheets("Simpsons Rule").Select' to calculate the area below the risk curves
For p = 1 To 19 'the number of theta values
Range("A27"). Select
61
-------
ActiveCell. Value = p * 0.05 'updates current threshold value ("theta")
Range("A29").Copy 'copies the cdf values associated with the current theta
Cells(p + 1, 2).PasteSpecial Paste:=xlPasteValuesAndNumberFormats 'paste cdf into column B
Next p
Range("C25:D25").Copy 'copies current concentration and normalized area value
Cells(k + 2, 7).PasteSpecial Paste: =xlPasteValuesAndNumberFormats
Range("C27").Copy 'copies cone header
Cells(l, 13 + k).PasteSpecial Paste:=xlPasteValuesAndNumberFormats
Range("B2:B20").Copy 'copies threshold curve data
Cells(2, 13 + k). PasteS pecial Paste:=xlPasteValuesAndNumberFormats
Sheets("Simpsons Rule"). Visible = False 'hides the worksheet
Sheets("Input Parameters"). Select
Range("N4").FormulaRlCl = "=now()" 'updates the elapsed time at the end of each cone.
exposure
'Sheets("Input Parameters"). Select
Application. ScreenUpdating = True 'allows the screens to update between each exposure cone
Application. ScreenUpdating = False
Nextk
PVA_SOLVER 'calls macro that fits curve to dose-response data
Range("I23").Select 'hides "selected" cell, just cosmetic
ActiveSheet.Protect
End Sub
Sub ClearPreviousQ
'this macro clears the data from the previous run
Sheets("mysid matrix 13x13").Select
Dim Count
Count = Range("U2") 'the count of the number of lambda values from the previous run
Range("W3").Resize(Count, 2).ClearContents
Range("AD2:AH51").ClearContents 'clears previous lambda values
Range("U5").ClearContents 'zeros out the % complete datum
'these next lines enable the "risk curve area" data to be cleared from the previous run
Sheets("Simpsons Rule").Visible = True
Sheets("Simpsons Rule").Select
Range("B2:B20").ClearContents 'clears last curve's information
Dim Col_Count
Col_Count = Range("J3") 'the number of risk curves from previous run
Range("G2").Resize(Col_Count, 2).ClearContents 'clears concentration vs area information
Range("Ml").Resize(20, Col_Count).ClearContents 'clears all of the curves' data
Sheets(" Simpsons Rule").Visible = False
62
-------
Sheets("Lambda Summary"). Select
Cells. Clear 'clears all of the lambda summary data from the previous run
Sheets("Input Parameters"). Select
'Change line color of PVA trendline to match background
ActiveSheet.ChartObjects("PVAchart").Activate
ActiveChart. SeriesCollection(3). Select
With Selection.Border
.LineStyle = xlNone
End With
'ActiveSheet.ChartObjects("DR_Curve").Activate
'ActiveChart. SeriesCollection(3). Select
'With Selection.Border
' .LineStyle = xlNone
'End With
Range("P25").Font.ColorIndex = 35 'change PC20 color to match background
Range("P26").Font.ColorIndex = 2 'color to match background
Range("I23").Select 'hides selected cell, just cosmetic
Application. ScreenUpdating = True
Application. ScreenUpdating = False
End Sub
Sub PVA_SOLVER()
'this macro runs the solver function to fit the dose-response curve
Sheets("PVA-Solver").Visible = True 'unhides worksheet
Sheets("PVA-Solver"). Select
Range("Kl").FormulaRlCl = "1.00" 'makes sure each run starts from same value
Range("K2").FormulaRlCl = "-1.00" 'makes sure each run starts from same value
SolverOk SetCell:="$H$5", MaxMmVal:=l, ValueOf:="0", ByChange:="$K$l:$K$2"
SolverSolve (True)
'this portion of the macro reset colors of endpoints to make them visable
Sheets("Input Parameters"). Select
'Change line color of PVA trendline to match markers
ActiveSheet.ChartObjects("PVAchart").Activate
ActiveChart. SeriesCollection(3). Select
With Selection.Border
.Colorlndex = 5
.Weight = xlMedium
.LineStyle = xlContinuous
End With
Range("N4").FormulaRlCl = "=now()" 'sets the ending time so that duration of run
'can be calculated
Range("N4").Copy '"freeze" end time
Range("N4"). Select
63
-------
With Selection
.PasteSpecial Paste: =xl Values
End With
Range("P25:P26").Font.ColorIndex = 1 'change PCx font color to black
Sheets("PVA-Solver").Visible = False 'hides worksheeet
Sheets("Input Parameters").Select 'just to make sure that model ends on correct sheet
End Sub
Sub LC50_kinetics_goalseek()
'this macro sets the kinetic factor when two different LC50s are known
Sheets("parameter calculations").Select
Range("Q4").Select
If ActiveCell. Value = True Then
Range("B51"). Select
Range("B51").FormulaRlCl = "1.0"
Range("G53").GoalSeekGoal:=Range("G51"), ChangmgCell:=Range("B51")
End If
End Sub
THIS MACRO IS CONTAINED WITHIN THE MODULE: ERROR_Check
Sub Date_Entry_Error_Check()
Range("I9"). Select
If ActiveCell. Value = "ERROR" Then
MsgBox "There is at least one data entry error in the main set of five" _
& vbCr & "survival parameters. Make sure that all are positive numbers," _
& vbCr & "and that the lower CL and upper CL are less than and greater" _
& vbCr & "than the 96h LC50, respectively."
Exit Sub
End If
Range("Ill"). Select
If ActiveCell. Value > 0 Then
MsgBox "There is at least one survival kinetics data entry error." _
& vbCr & "Make sure there are no blanks and that all are positive numbers." _
& vbCr & "If you selected the 'two' LC50 method, then make sure that the" _
& vbCr & "earlier LC50 is greater than the later one."
Exit Sub
End If
64
-------
Range("I22"). Select
If ActiveCell. Value > 0 Then
MsgBox "There is at least one reproduction data entry error." _
& vbCr & "Make sure there are no blanks and that all are positive numbers."
Exit Sub
End If
Range("I3 5"). Select
If ActiveCell. Value = "ERROR" Then
MsgBox "There is at least one exposure parameter data entry error." _
& vbCr & "Make sure there are no blanks and that all are positive numbers."
Exit Sub
End If
End Sub
THESE MACROS ARE CONTAINED WITHIN THE MODULE: Move_Macros
Sub gotoDose_Response_Input_Graphs()
Application. ScreenUpdating = False
LC50_kinetics_goalseek 'resets kinetic parameter—only needed with "2 LC50" method
Sheets("Input Parameters"). Select
ZoomToRange ZoomThisRange:=Range("AA1: AN1"), PreserveRows:=False
Range("AAl:AN32").Select
ActiveWindow.Zoom = True
Range("Ac3"). Select
End Sub
Sub gotoInput_Parameter_Area()
Range("Al:S32"). Select
ActiveWindow.Zoom = True
Range("I23"). Select
End Sub
Sub ZoomToRange(ByVal ZoomThisRange As Range, _
ByVal PreserveRows As Boolean)
This procedure was adapted from "Zooming On A Range" —www/cpearson. com/excel/topic.aspx
The author of this code has declared it "Public Domain" (see
www/cpearson.com/Excel/LegaleseAndDisclaimers.aspx)
ActiveSheet.Unprotect
Dim Wind As Window
Set Wind = Active Window
Application. ScreenUpdating = False
65
-------
' Put the upper left cell of the range in the top-left of the screen.
i
Application. Go To ZoomThisRange(l, 1), True
With ZoomThisRange
If PreserveRows = True Then
.Resize(.Rows. Count, 1).Select
Else
.Resize(l, .Columns.Count).Select
End If
End With
With Wind
.Zoom = True
.VisibleRange(l, 1). Select
End With
ActiveSheet.Protect
End Sub
THESE MACROS ARE CONTAINED WITHIN THE MODULE: Print_Save_Macros
Function PathExists(MyPath) As Boolean
'returns TRUE if the path exists
'--saw this in several places on web, but also in
'Walkenbach, John.l999.Microsoft Excel 2000 Power Programming with VBA
'Published by Hungry Minds, Inc, New York, NY. 869 pp.
Dim x As String
On Error Resume Next
x = GetAttr(MyPath) And 0
If Err = 0 Then PathExists = True _
Else: PathExists = False
End Function
Sub printSummaryQ
'prints Input Parameter and endpoint graph to the current default printer
Application. ScreenUpdating = False
ActiveSheet.PageSetup.PrintArea = "A1:T31"
ActiveWindow.SelectedSheets.PrintOut From:=l, To:=l, Copies:=l, Collate_
:=True
ActiveSheetPageSetup.PrintArea =""
With ActiveSheet.PageSetup
.FitToPagesWide = 1
.FitToPagesTall = 1
End With
End Sub
Sub CopyPasteValuesQ
ActiveSheet.Unprotect
66
-------
' this macro clears out the extraneous materials from the parameter summary sheet
' which being saved as a summary of a model run
Cells. Select
Selection. Copy
Selection.PasteSpecial Paste:=xlPasteValuesAndNumberFormats
ActiveSheet.Shapes("ViewBox").Cut
ActiveSheet. Shapes("RunBox"). Cut
ActiveSheet.Shapes("ReturnBox").Cut
ActiveSheet.Shapes("PrintBox").Cut
ActiveSheet. Shapes(" SaveBox").Cut
ActiveSheet. ChartObjects("DoseResponse"). Activate
ActiveChart. ChartArea. Select
Active Window. Visible = False
Selection.Delete
ActiveSheet. ChartObjects("LC50kinetics").Activate
ActiveChart. ChartArea. Select
Active Window. Visible = False
Selection.Delete
Columns("Z:AK").Clear
Range("M3"). Select
End Sub
Sub CopyPasteChartQ
'the macro converts the endpoint chart to a picture
Dim TheChart As Chart
Application. ScreenUpdating = False
ActiveSheet. ChartObjects("PVAchart").Activate
Set TheChart = ActiveChart
TheChart. CopyPicture appearance: =xlScreen, Size:=xlScreen, Format: =xlPicture
Active Window. Visible = False
ActiveSheet. Paste
Selection.Left = TheChart. Parent. Left
Selection. Top = TheChart. Parent. Top
ActiveSheet. ChartObjects("PVAchart").Delete
Range("M3"). Select
End Sub
67
-------
Sub Save_Previous_Run()
Dim modelWB As Workbook
Set modelWB = Active Workbook This assigns the model filename to a variable—allows this
'macro to still opperate even if the filename is changed
Dim default_Path
default_Path = Range("Saving_Default")
This is the default path to which the results will be saved.
It's location is on the "Input Parameter" sheet off screen
'below the save box.
Dim MyFile, MyPath
MyPath = InputBox("Pathname for Summary Information", "Enter Pathname—don't forget the
backslash at the end", _
default_Path)
'You can accept the default path or enter a new one.
If PathExists(MyPath) = False Then
ans = MsgBox("Path does not exist. Create?", vbYesNo, "Path Selected")
If ans = vbNo Then Exit Sub
If ans = vbYes Then MkDir MyPath 'If the path does not exist, it will be created.
End If
MyFile = InputBox("Filename for Summary Information", "Enter UNIQUE filename",
"PVA_Summary.xls")
Dim MySummary As Workbook
Set MySummary = Workbooks. Add 'adds a workbook to serve as the summary file
MySummary. SaveAs (MyPath & MyFile) 'names the new workbook
Application. ScreenUpdating = False
modelWB. Activate
Sheets("Lambda Summary").Copy before:=Workbooks(MyFile).Sheets(l)
modelWB. Activate
Sheets("Input Parameters").Copy before: =Workbooks(MyFile).Sheets("Lambda Summary")
Sheets("Input Parameters").Name = "Summary of Model Parameters" 'renames sheet
These next two macros clean up the summary file
CopyPasteValues
CopyPasteChart
Sheets("Sheetl"). Visible = False
Sheets("Sheet2"). Visible = False
Sheets(" Sheets"). Visible = False
MySummary. Save
MySummary. Activate
Application. ScreenUpdating = True
End Sub
68
-------
References
Ak9akaya, HR. 1992. Population viability analysis and risk assessment, p. 148-17 (in) DR
McCullough and RH Barrett (eds). Wildlife 2001: Populations. Elsevier Applied Science.
1163pp.
Ak9akaya, HR. 2000. Population viability analyses with demographically and spatially
structured models. Ecological Bulletins 48:23-38.
Ak9akaya, HR. 2002. Estimating the variance of survival rates and fecundities. Animal
Conservation 5:333-336.
Barnthouse, LW, GW Suter II and AE Rosen. 1990. Risks of toxic contaminations to exploited
fish populations: influence of life history, data uncertainty and exploitation intensity.
Environmental Toxicology and Chemistry 9:297-311.
Bliss, CI. 1934. The method of probits. Science 79:38-39 and correction on pages 409-410.
Bliss, CI. 1935. The calculation of the dosage-mortality curve. Annals of Applied Biology
22:134-167.
Burgman, MA, S Person and HR Ak9akaya. 1993. Risk Assessment in Conservation Biology.
Chapman & Hall. 314 pp.
Campbell, TC. 2005. An introduction to clinical significance: An alternative index of
intervention effect for group experimental designs. Journal of Early Intervention 27:210-
227.
Cardin, J. 1986. Ammonia acute and chronic tests with Mysidopsis bahia. Memorandum to D
Hansen. Octobeer 15. 5 pp. Available from the US EPA, Atlantic Ecology Division,
Narragansett, RI.
Caswell, H. 1989. Analysis of life table response experiments. I. Decomposition of effect on
population growth rate. Ecological Modelling 46:221 -237'.
Caswell, H. 1996. Demography meets ecotoxicology: Untangling the population level effects of
toxic substances, (in) MS Newman and CH Jagoe (eds) Ecotoxicology: A Hierarchical
Treatment. CRC Press. 441 pp.
Caswell, H. 2001. Matrix Population Models: Construction, Analysis, and Interpretation. 2nd Ed.
Sinauer Associates. 722 pp.
Champlin, D and S Poucher. 199la. Flow-through chronic toxicity of fluoranthene to Mysidopsis
bahia. Memorandum to DJ Hansen. Environmental Research Laboratory-Narragansett.
October 1. 2 pp. Available from the US EPA, Atlantic Ecology Division, Narragansett,
RI.
Champlin, D and S Poucher. 1991b. Acute toxicity test conducted with fluoranthene (dark).
Memorandum to DJ Hansen. Environmental Research Laboratory-Narragansett.
December 9. 2 pp. Available from the US EPA, Atlantic Ecology Division, Narragansett,
RI.
69
-------
Champlin, D and S Poucher. 1992a. Chronic toxicity of pyrene to Mysidopsis bahia in a 28 day
flow-through test. Memorandum to S. Lussier. Environmental Research Laboratory-
Narragansett. September 15. 4 pp. Available from the US EPA, Atlantic Ecology
Division, Narragansett, RI.
Champlin, D and S Poucher. 1992b. Flow-through chronic toxicity of carbofuran to Mysidopsis
bahia. Memorandum to S. Lussier. Environmental Research Laboratory-Narragansett.
September 15. 2 pp. Available from the US EPA, Atlantic Ecology Division,
Narragansett, RI.
Cohen, J. 1990. Things I have learned (so far). American Psychologist 45(12)1304-1312.
Di Stefano, J, F. Fidler, G. Gumming. 2005. Effect size and confidence intervals: An alternative
focus for the presentation and interpretation of ecological data. pp. 71-102. (in) AR Burk
(ed). New Trends in Ecological Research. Nova Science, New York.
Evans, M, N Hasting and B Peacock. 2000. Statistical Distributions. 3r Edition. John Wiley &
Sons. 221 pp.
Finney, DJ. 1971. Probit Analysis, 3r Ed. Cambridge Univ. Press. 333 pp.
Gentile, JG, SM Gentile, NG Hairston and BK Sullivan. 1982a. The use of life-tables for
evaluating the chronic toxicity of pollutants to Mysidopsis bahia. Hydrobiologia 93:179-
187.
Gentile, JG, SM Gentile, G Hoffman, JF Heltshe and NG Hairston. 1983. The effects of a
chronic mercury exposure on survival, reproduction and population dynamics of
Mysidopsis bahia. Environmental Toxicology and Chemistry2:6l-68.
Gentile, SM, JH Gentile, J Walker and JF Heltshe. 1982b. Chronic effects of cadmium on two
species of mysid shrimp: Mysidopsis bahia and Mysidopsis bigelowi. Hydrobiologia
93:195-204.
Ginzburg, LR, LB Slobodkin, K Johnson and AGBindman. 1982. Quasiextinction probabilities
as a measure of impact on population growth. Risk Analysis 2:171-181.
Gleason, TR and DE Nacci. 2001. Risks of endocrine-disrupting compounds to wildlife:
Extrapolation from effects on individuals to population response. Human and Ecological
Risk Assessment. 7(5): 1027-1042.
Hubert, JJ. 1980. Bioassay, 2nd Ed. Kendall/Hunt Pub. 180 pp.
Inchausti, P and H Weimerskirch. 2004. Wandering albatross (Diomedea exulans chionoptera) in
the southern oceans, pp. 421-430. (in) HR Ak9akaya, MA Burgman, O. Kindvall, CC
Wood, P Sjogren-Gulve, JS Hatfield and MA McCarthy (eds). Species Conservation and
Management: Case Studies. Oxford Univ. Press. 533 pp.
IUCN (World Conservation Union). 2001. IUCN Red List categories and criteria: version 3.1.
IUCN Species Survival Commission. IUCN, Gland, Switzerland.
Kammenga, JE, M. Dusschers, NM Van Straalen, PC Jepson, and J Bakker. 1996. Stress induced
fitness reduction is not determined by the most sensitive life-cycle trait. Functional
Ecology 10:106-111.
70
-------
Kendall, BE. 1998. Estimating the magnitude of environmental stochasticity in survivorship
data. Ecological Applications 8:184-193.
Kirk, RE. 1996. Practical significance: A concept whose time has come. Educational and
Psychological Measurement 56:746-759.
Kuhn, A and S Lussier. 1987. Phenanthrene results of acute and chronic tests (flow-thru) with
Mysidopsis bahia. Memorandum to DJ Hansen. Environmental Research Laboratory-
Narragansett. August 3. 3 pp. Available from the US EPA, Atlantic Ecology Division,
Narragansett, RI.
Kuhn, A, S Lussier and GB Thursby. 1987. Memorandum to DJ Hansen. Results of chronic life
cycle, mysid test conducted with phenol at ERL, Narragansett. August 3. 3 pp. Available
from the US EPA, Atlantic Ecology Division, Narragansett, RI.
Kuhn, A, WR Munns, S Poucher, D Champlin and S Lusssier. 2000. Prediction of population-
level response from mysid toxicity test data using population modeling techniques.
Environmental Toxicology and Chemistry. 19(9):2364-2371.
Kuhn, A, WR Munns, D Champlin, R McKmney, M Taghabue, J Serbst and T Gleason. 2001.
Evaluation of the efficacy of extrapolation population modeling to predict the dynamics
of Americamysis bahia populations in the laboratory. Environmental Toxicology and
Chemistry. 20(1 ):213-221.
Lande, R, S. Engen and B-E Sasther. 2003. Stochastic Population Dynamics in Ecology and
Conservation. Oxford Univ. Press. 222 pp.
Lee, ET and JW Wang. 2003. Statistical Methods for Survival Data Analysis. 3rd Edition. Wiley
Series in Proability and Statistics. Wiley & Sons, Inc. 513 pp.
Lussier, SM, JH Gentile and J Walker. 1985. Acute and chronic effects of heavy metals and
cyanide onMysidopsis bahia (Crustacea: Mysidacea). Aquatic Toxicology!: 25-35.
Lussier, SM, A Kuhn, MJ Chammas and J Sewall. 1988. Techniques for the laboratory culture of
Mysidopsis species (Crustacea: Mysidacea). Environmental Toxicology and Chemistry
7:969-977.
McCarthy, MA. 1996. Red kangaroo (Macropus rufus] dynamics: effects of rainfall, density
dependence, harvesting and environmental stochasticity. Journal of Applied Ecology
33:45-53.
McCarthy, MA and LS Broome. 2000. A method for validating stochastic models of population
viability: a case study of the mountain pygmy-possum (Burramys parvus). Journal of
Animal Ecology 69:599-607.
McCarthy, MA and C Thompson. 2001. Expected minimum population size as a measure of
threat. Animal Conservation 4:351-355.
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.
71
-------
Morris, WF and DF Doak. 2002. Quantitative Conservation Biology: Theory and Practice of
Population Viability Analysis. Sinauer Associates. 480pp.
OPP. 1993. US Environmental Protection Agency, Office of Pesticide Programs. Data
Evaluation Record. MRID No. 421449-06. Mysid Chlorpyrifos Acute
OPP. 1994. US Environmental Protection Agency, Office of Pesticide Programs. Data
Evaluation Record. MRID No. 426649-01. Mysid Chlorpyrifos Chronic
Reed, JM and AR Blaustein. 1997. Biological significant population declines and statistical
power. Conservation Biology ll(l):281-282.
Raimondo, S and CL McKenney, Jr. 2005a. Projected population-level effects of thiobencarb
exposure on the mysid, Americamysis bahia, and extinction probability in a concentration-
decay exposure system. Environmental Toxicology and Chemistry. 24(3):564-572.
Raimondo, S and CL McKenney, Jr. 2005b. Projecting population-level responses of mysids
exposed to an endrocrine disrupting chemical. Integrative and Comparative Biology 45:151-
157.
Ruel, JJ and MP Ayers. 1999. Jensen's inequality predicts effects of environmental variation.
Trends in Ecology and Evolution. 14(9): 3 61 -3 66.
Schimmel, SC. 1981. Results: Interlaboratory Comparison—Acute Toxicity Tests Using
Estuarine Animals. U.S. EPA Report # EPA-600/4-81-003/
Schtickzelle, N, J Choutt, P Goffart, V Fichefet and M Baguette. 2005. Metapopulation
dynamics and conservation of the marsh fritillary butterfly: Population viability analysis and
management options for a critically endangered species in Western Europe. Biological
Conservation 126:569-581.
Sprague, JB. 1969. Measurement of pollutant toxicity to fish—I. Bioassay methods for acute
toxicity. Water Research 5:793-821.
Springborn. 1988a. Acute toxicity of 1,2,4-trichlorobenzene to mysid shrimp (Mysidopsis bahia)
under flow-through conditions. Report submitted to the Chemical Manufacturers
Association, Washington, DC. June 28. 36 pp.
Springborn. 1988b. Chronic toxicity of 1,2,4-trichlorobenzene to mysid shrimp (Mysidopsis
bahia}. Report submitted to the Chemical Manufacturers Association, Washington, DC.
July 26. 57 pp.
Stephan, CE, DI Mount, DJ Hansen, JH Gentile, GA Chapman and WA Brungs. 1985.
Guidelines for deriving numerical national water quality criteria for the protection of
aquatic organisms and their uses. EPA Report 822/R-85-100 (NTIS Report PB85-
227049). 98 pp.
Taylor, RJ, T Regan, M Burgman and K Bonham. 2003. Impacts of plantation development,
harvesting schedules and rotation lengths on the rate snail Tasmaphena lamproides in
northwest Tasmania: a population viability analysis. Forest Ecology and Management
175:455-466.
72
-------
Thursby, GB and WJ Berry. 1987a. Results of acute toxicity tests conducted with phenol at ERL,
Narragansett by the URI cooperative agreement. Memorandum to DJ Hansen. July 31.9
pp. Available from the US EPA, Atlantic Ecology Division, Narragansett, RI.
Thursby, GB and WJ Berry. 1987b. Acute and chronic toxicity of aniline to Mysidopsis bahia:
Flow-through. Memorandum to DJ Hansen. September 19. 6 pp. Available from the US
EPA, Atlantic Ecology Division, Narragansett, RI.
Thursby, GB and WJ Berry. 1988. Acute and chronic toxicity of thallium to Mysidopsis bahia:
Flow-through. Memorandum to J Scott (SAIC) and DJ Hansen (EPA). May 20. 7 pp.
Available from the US EPA, Atlantic Ecology Division, Narragansett, RI
Thursby, GB, WJ Berry and D Champlin. 1989. Flow-through acute and chronic tests with
acenaphthene using Mysidopsis bahia. Memorandum to DJ Hansen. September 19. 5 pp.
Available from the US EPA, Atlantic Ecology Division, Narragansett, RI.
Thursby, GB and D Champlin 1990. Acute and chronic toxicity of dichlorvos to Mysidopsis
bahia. Memorandum to DJ Hansen. May 9. 5 pp. Available from the US EPA, Atlantic
Ecology Division, Narragansett, RI.
Thursby, GB and D Champlin. 1991. Flow-through acute and chronic toxicity of carbaryl to
Mysidopsis bahia. Memorandum to DJ Hansen. Environmental Research Laboratory-
Narragansett. June 13. 3 pp. Available from the US EPA, Atlantic Ecology Division,
Narragansett, RI.
Thursby, GB, D Champlin and WJ Berry. 1990. Acute and chronic toxicity of propoxur to
Mysidopsis bahia. Memorandum to DJ Hansen. Environmental Research Laboratory-
Narragansett. September 16. 4 pp. Available from the US EPA, Atlantic Ecology
Division, Narragansett, RI.
Urban, DJ and NJ Cook. 1986. Hazard Evaluation Division Standard Evaluation Procedure:
Ecological Risk Assessment. EPA Report 540/9-85-001. 96pp.
US EPA. 2005. Draft NHEERL SP2 Implementation Plan. Appendix B: SP2—LTG2 Program
Project Descriptions. August 30.
US EPA. 2006. Multi-Year Plan for Safe Pesticides/Safe Products: 2007-2015. Office of
Research and Development. December.
Ward, GS and L Ballantine. 1985. Acute and chronic toxicity of atrazine to estuarine fauna.
Estuaries 8:22-27.
Ward, GS, TA Hollister, PT Heitmuller and PR Parrish. 1981. Acute and chronic toxicity of
selenium to estuarine organisms. Northeast Gulf Science 4:73-78.
73
------- |