FDA Logo U.S. Food and Drug AdministrationCenter for Food Safety and Applied Nutrition
U.S. Department of Health and Human Services
horizontal rule

July 19, 2005

horizontal rule

Quantitative Risk Assessment on the Public Health Impact of
Pathogenic Vibrio parahaemolyticus in Raw Oysters

Table of Contents

Appendix 5: Details of the Data Analysis for the Exposure Assessment Component of the Vibrio parahaemolyticus Risk Assessment Model

Relationship between Levels of Vibrio parahaemolyticus in Oysters At-Harvest and Environmental Conditions

There have been a number of extensive studies conducted over a wide range of geographic locations showing the relationship of environmental factors and total V. parahaemolyticus levels in water and oysters. These studies were reviewed and evaluated here with regard to their utility for developing or estimating an appropriate predictive relationship between total V. parahaemolyticus densities in oysters at the time of harvest and environmental conditions; specifically water temperature and salinity. Most of the older studies did not provide sufficient information with respect to a quantitative relationship, primarily because these studies were either limited to specific seasons with correspondingly little variation of environmental parameters, measured V. parahaemolyticus levels in water or sediment rather than oysters or reported little quantitative data on densities per se. The following sixteen studies that were evaluated are listed below:

Harvest Module

Water Temperatures
With the exception of the Pacific Northwest region, distributions of regional/seasonal water temperatures were developed based on accumulated records from selected coastal water buoys maintained by the National Buoy Data Center (NBDC). Hourly water temperature measurements were generally available from 1984 through 1998 from several buoys in each region. However, given intermittent records and lack of both water temperature and air measurements for some buoys, a single representative buoy was selected for each region that had both water and air temperature measurements. The data from these selected buoys were analyzed to determine an appropriate summary distribution for temperature in the Monte Carlo simulation model. After examination of these data, implementation of temperature distributions in the Monte Carlo simulation by resampling from empirical distributions was determined to be overly cumbersome. Although there is some error associated with simpler distribution summaries that have been used, and are discussed here, the differences appear to be minor in consideration of the natural variation of V. parahaemolyticus levels at any given temperature and the other factors/uncertainties identified in the risk assessment process.

Although there is a diurnal cycle in water (and air) temperature, the effects of hourly changes in water temperatures were not considered in predicting V. parahaemolyticus levels at the time of harvest. Examination of selected NBDC buoy datasets indicated that the hourly water temperature variations were minor in comparison to the variations across days or weeks. This is illustrated in Figure A5-1, which shows the mean, the 2.5% and 97.5%-percentiles of the hourly water temperature measurements recorded at the NBDC Dauphin Island, Alabama buoy during the summer of 1997. As is evident in the figure, the variation of mean water temperature by time of the day is only slightly greater than 1 °C; much less than the variation across different days or weeks as indicated by the percentiles.

The day-to-day variation in temperature is temporally correlated as weather patterns determining air and water temperatures persist over time spans varying from several days to several weeks. Figure A5-2 shows the temporal pattern of daily (midday) water temperatures recorded at Dauphin Island in the summer of 1997. Figure A5-3 shows a histogram plot of the same temperature data with an approximate normal distribution summarizing the variation about the mean. A normal distribution was fit to the data by the method of "moments". As evident in Figure A5-3, the actual temperature data are skewed and the normal distribution summary does not capture this facet of the data since only the 1st two moments of the fitted distribution match that of the empirical distribution.

Image of Figure A5-1. Mean and Percentiles (2.5% and 97.5%) of Hourly Water Temperature Profile for Dauphin Island, AL (NBDC, July - Sept 1997)

Figure A5-1. Mean and Percentiles (2.5% and 97.5%) of Hourly Water Temperature Profile for Dauphin Island, AL (NBDC, July - Sept 1997)


Image of Figure A5-2. Temporal Pattern of Day-to-Day Variation of MiddayWater Temperature Profile for Dauphin Island, AL (NBDC, July - Sept 1997)

Figure A5-2. Temporal Pattern of Day-to-Day Variation of MiddayWater Temperature Profile for Dauphin Island, AL (NBDC, July - Sept 1997)


Image of Figure A5-3. Histogram of Day-to-Day Variation of Midday Water Temperature Profile and Approximating Normal Distribution Summary for Dauphin Island, AL (NBDC, July - Sept 1997)

Figure A5-3. Histogram of Day-to-Day Variation of Midday Water Temperature Profile and Approximating Normal Distribution Summary for Dauphin Island, AL (NBDC, July - Sept 1997)

Overall, considering the patterns in the observed distributions of (midday) water temperatures in the other regions and in other seasons for all years of data available from NBDC (1984 - 1998), the normal distribution approximation (as shown in Figure A5-3) was judged to be a reasonable summary for within season water temperature variation. This summary distribution was chosen for simplicity and in consideration of larger determinant factors and uncertainties in the risk assessment model. Within season air temperature, and consequently water temperature of shallow water bodies, is known to be slightly skewed, primarily as a consequence of precipitation patterns. This is reflected in more complex temporal modeling of weather patterns and their effect on temperature (Richardson, 1981). The skewness in the NBDC temperature data for other regions and seasons, when present in a given year, is typical of that shown in Figure A5-3. For these data the approximating Gaussian distribution underestimates the median temperature slightly (approximately 1 °C). Within the context of the risk assessment model, the effect of this bias would be to underestimate levels of total V. parahaemolyticus at harvest.

In addition to the variation of daily water temperature in the NBDC data, the variation of water temperature distributions across multiple years was also evaluated for incorporation as a factor in the risk assessment. Figure A5-4 shows a plot of the means and standard deviations of within year and season daily (midday) water temperature distributions across multiple years of data for the Dauphin Island buoy. The data represented here are from 1989 through 1998 (with 1995 excluded due to instrument malfunction). There are four clusters of points based on the definitions of the four seasons used to categorize the data. As evident in Figure A5-4, temperature distributions were more variable in the spring and fall (middle two clusters of points in the plot) and the least variable in the summer.

Image of Figure A5-4. Interannual Variation of the Mean and Standard Deviation of Within Season Water Temperature Distributions for Dauphin Island, AL (NBDC, 1997).

Figure A5-4. Interannual Variation of the Mean and Standard Deviation of Within Season Water Temperature Distributions for Dauphin Island, AL (NBDC, 1997).

The data suggest that a relationship between the mean and variance of daily water temperatures applies not only to comparison of temperature variation in one season versus another within the same year but also to the variation of temperature in the same season across different years. That is, there appears to be a tendency for a warmer than average summer to be less variable day-to-day than a cooler than average summer. Rather than using an approximating relationship to model or summarize this characteristic of the data (e.g., such as by a quadratic), the year-to-year variation of the means and standard deviations of the within year seasonal distributions were summarized by means, variances and correlations. These summaries are shown in Table IV-2.

The number of years of water temperature data available from NDBC sites was limited to at most 15. The statistical precision associated with estimates of the variation and correlation between the mean and standard deviation across different years is low. Nevertheless, some degree of correlation seems reasonable a priori, and therefore the apparent correlations were used in Monte Carlo simulations as part of evaluating the effect of year-to-year variations in temperature on the total illness rate predictions. To obtain model predictions, the parameters for within year and season water temperature distributions were obtained by sampling from the bivariate normal distributions with means, standard deviations and correlation specified in Table IV-2. The parameters of the bivariate normal distributions were obtained by a method of moments fit to the relevant summary statistics.

NDBC water temperature data were not available for the Pacific Northwest region. However, the same analysis was conducted using a set of approximately 50,000 water temperatures measurements collected from 1988 through 1999 during the course of various harvest water monitoring programs (e.g., for fecal coliforms, vibrios, etc). These data were made available to the V. parahaemolyticus QRA Team by the Washington State Department of Health (1999). The temperature measurements were generally taken at the time of sample collection (e.g., for water or shellfish). The monitoring programs from which these data were abstracted were conducted in multiple oyster harvesting areas of Washington State. Frequently, the dataset had multiple measurements from different sites within the same estuary on the same day. The time of temperature measurement was not reported. Prior to developing the summary statistic described above, the data were averaged over measurements in the same estuary on the same day. It was assumed that variation of time of water temperature measurement across the data set was of minor concern (i.e., in consideration of the discussion of Figure A5-1 above).

Total V. parahaemolyticus per gram versus Water Temperature Relationship
Given water temperature distributions, the prediction of the distribution of the density of total V. parahaemolyticus in oysters at the time of harvest was obtained based on fitted regression relationships between V. parahaemolyticus density and water temperature. Data from three sources were evaluated to determine the relationship. These data were considered the most appropriate for determining the regression relationship because oyster sampling was conducted year round. Two of the studies (DePaola et al., 1990, FDA/ISSC, 2001) were nationwide studies with samples being obtained from geographically diverse harvest areas.

Additional data, specific to the Pacific Northwest, were made available to the V. parahaemolyticus QRA Team. These data were collected during Washington State monitoring programs at selected times from 1997 through 2001 ((Washington State Department of Health, 2000; 2001). Because the ecology of vibrios is notably different in the Pacific Northwest (e.g., possibly as a consequence of higher salinities) and this area is underrepresented in the two studies above, a separate analysis of the V. parahaemolyticus versus water temperature relationship was considered appropriate for this region. However, it was apparent that data collected by Washington State were influenced by the nature of commercial harvesting (e.g., intertidal collection versus dredged). Examination of the Washington State data shows high levels of V. parahaemolyticus/g for areas of Hood Canal and South Puget Sound during the summer. These areas are known for intertidal harvesting and the high observed densities are not appropriately predicted based on water temperature alone.

With respect to other regions, the risk assessment model is structured to define "at-harvest densities" as those occurring when oysters are submerged in water. This definition is less obvious for the Pacific Northwest but, to maintain consistency, "at-harvest densities" are defined to be the same as that of the other regions. By this definition, elevation in V. parahaemolyticus/g due to intertidal collection is a Post Harvest effect and is most appropriately modeled as separate from the effect of water temperature. Thus, an estimate of the relationship of V. parahaemolyticus/g versus water temperature in "at-harvest" (i.e., submerged) oysters was considered a necessary component of the risk assessment model construction.

In this regard, it was determined that subsetting of the data from Washington State was appropriate by excluding from "at-harvest" estimation any data that was likely to have been collected intertidally. Areas of Hood Canal and Southern Puget Sound in the Washington State dataset were considered predominantly intertidal harvest areas and were therefore excluded from "at-harvest" density estimation. The remaining data collected from harvest areas in Willapa Bay and Northern Puget Sound were considered to be predominantly dredged and thus the most appropriate with respect evaluating the relationship between water temperature and "at-harvest" densities (i.e., in order to predict V. parahaemolyticus/g in submerged oysters). Appropriate modeling of the effect of the intertidal collection on V. parahaemolyticus/g is discussed further under the sections of this appendix pertaining to modeling of the Post Harvest module.

All three of the data sets that were evaluated are subject to some limitations of measurement. When water temperatures are low, total V. parahaemolyticus levels are generally below the limit of detection of currently available methods (MPN procedure, direct plating). Thus, microbiological analysis of the oyster samples obtained during the winter frequently yields an outcome that is "nondetect" or nondetectable (i.e., the microorganism is not found in the samples). Statistically, the outcome of such measurements are said to be left censored at the limit of detection (LOD), since failure to isolate V. parahaemolyticus from a relatively small analytical portion of a sample places an upper limit on the density in the sample rather than an estimate of that density per se. Depending on the size of analytical portion relative to the sample, a conclusion that the microorganism is not present in the sample is generally not warranted. Appropriate statistical analysis of data sets with censored values depends on the overall extent or proportion of samples that are censored. When the extent of censoring is relatively small (<10% of samples), a mean imputation of half the LOD is a commonly used strategy and has been shown not to overly bias parameter estimates. However, when the extent of censoring is large (e.g., >40%), it is generally accepted that mean imputation is not appropriate (EPA, 2000). The degree of censoring is approximately 40% for all three of the data sets considered here. Thus, given that mean imputation is questionable, alternative methods of analysis were considered (see for example, Carabin et al., 2001).

The method of analysis used here to determine the regression relationship is the censored or Tobit regression method (Tobin, 1958). The method is appropriate when censoring occurs in the context of regression with normally distributed and homogeneous variance about the mean. Following the observation that the distribution of V. parahaemolyticus densities at a given temperature is positively skewed and asymmetric, the log10 transformation of densities is approximately normally distributed. This is not particular to V. parahaemolyticus but is common to exposure assessment of other microbial pathogens, possibly due to the exponential characteristics of birth and death of microbial populations. Thus, the Tobit regression method was found to apply to an analysis of log10 V. parahaemolyticus/g versus water temperature and was chosen because it is a relatively simple procedure.

Assuming a linear relationship between mean log10 V. parahaemolyticus/g and water temperature, both above and below the censoring point or limit of detection (LOD), the regression model assumed for parameter estimation is of the form

mean log10 V. parahaemolyticus/g = α + β *WTEMP + ε

where WTEMP is the water temperature, α and β are parameters for the linear relationship of log10 V. parahaemolyticus /g versus temperature, and ε is a normally distributed variate with mean zero and variance σ2. Estimation of parameters in the Tobit regression method is commonly obtained by the frequentist procedure of maximizing the likelihood function, which models the probability of obtaining nondetectable outcomes as well as quantified values. Specifically, when an observation is found to be quantifiable with value Xi (i.e., a value above the LOD) at a water temperature WTEMPi the contribution to the likelihood is

Image of the probability density function

and when an observation is nondetect (at a water temperature WTEMPi) the contribution to the likelihood is

Image of the cumulative distribution function

where φ is the probability density function and Φ is the cumulative distribution function of the standard normal. The total likelihood is the product of the likelihood components for each datum. Maximum likelihood estimates (MLEs) of the parameters are those values of the parameters that maximize this function (i.e., the values of the parameters for which the observed data are most probable in the context of the model).

The SAS (SAS Institute, 1999) procedure LIFEREG was used to obtain the estimates of the regression parameters based on these criteria for each of the three data sets being analyzed. The parameter estimates obtained are given in Table A5-1. The variance-covariance matrix of the MLEs obtained is given in Table A5-2. These variance-covariance estimates were used to construct an approximate uncertainty distribution for the parameter estimates. Based on asymptotic normal distribution theory, the parameter uncertainty was taken to be multivariate normal with mean equal to the MLEs obtained and variance-covariance matrix equal to that estimated based on the data.

Table A5-1. Maximum Likelihood Estimates of the Parameters of the Tobit Regression of log10 Vibrio parahaemolyticus per gram Versus Water Temperature (MLEs and 95% Confidence Intervals of Parameters of Temperature-only Regression)
Study α β σ
DePaola et al., 1990 -1.03 (-2.14,0.08) 0.12 (0.072,0.17) 1.07 (0.83,1.37)
FDA/ISSC, 2001 -0.63 (-0.87,-0.39) 0.10 (0.092,0.11) 0.76 (0.71,0.82)
Washington State DOH, 2000; 2001 -4.32 (-5.77,-2.88) 0.24 (0.16,0.32) 0.78 (0.61,0.99)

Table A5-2. Estimated Variance-Covariance Matrix for MLEs of log10 Vibrio parahaemolyticus per gram Versus Water Temperature Regression Parameters
Study   α β σ
DePaola et al., 1990 α 0.32 -0.014 -0.032
β -0.014 0.00062 0.0012
σ -0.032 0.0012 0.019
FDA/ISSC, 2001 α 0.015 -0.00063 0.0012
β -0.00063 0.000028 0.000044
σ 0.0012 0.000044 0.00081
Washington State DOH, 2000; 2001 α 0.543 -0.031 -0.029
β -0.031 0.0018 0.0015
σ -0.029 0.0015 0.0092

Alternative methods of parameter estimation that could have been applied to appropriately estimate parameter values in the presence of censoring include logistic regression of data classified as >LOD versus <LOD and multiple imputation of left censored data according to estimated distributions. Analyses of the three data sets by these methods do not give substantially different parameter estimates than that obtained by the Tobit regression method.

As assessed by likelihood ratio statistics and measures of goodness-of-fit appropriate to the Tobit regression model (Cameron and Windmeijer, 1996, 1997), the fits of the models to the data were good with temperature being a highly significant effect (p<0.001). Based on McFadden's R2 (a likelihood-based extension of the usual R2) which is appropriate to the Tobit model, the proportion of the variance in log10 V. parahaemolyticus densities which is explained by the effect of temperature is approximately 50%.

The effect of the uncertainties in the parameter estimates obtained by these regression analyses was incorporated into risk assessment evaluation. Although the uncertainty is ostensibly an uncertainty with respect to the relationship existing at the time samples were collected, this was assumed a reasonable surrogate for the potential variation of the relationship across different years due to the possibility of changes in other environmental conditions affecting V. parahaemolyticus densities (e.g. oyster physiology, oyster disease, nutrient levels). The effect of regression parameter uncertainty was implemented in the risk assessment by using a multivariate normal approximation for parameter uncertainty (i.e., asymptotic normality of MLEs with sufficiently large sample size). A multivariate approach was necessary due to the fact that the parameter estimates for the slope and intercept of the regressions were highly correlated. The effect of this uncertainty was implemented in Monte Carlo simulations by taking a sample of 1,000 sets of parameters from the uncertainty distributions (a multivariate normal with mean equal to the MLEs and variance-covariance equal to the estimated variance-covariance matrix).

Independent estimates of method error were then used to correct the estimated variance about the regression lines to predict the population variance of the density per gram. The FDA/ISSC (Cook et al., 2002b) study utilized a direct plating procedure with DNA probes. The method error variance associated with this method has been estimated to be 0.03 based on the difference between counts on replicate analyses of sample aliquots (Ellison et al., 2001). A method error variance of 0.03 for log10 V. parahaemolyticus (Vp)/g corresponds to a standard deviation of 0.17 log10 Vp/g between replicate analyses. The FDA-BAM method with 3 tubes per dilution was the standard method of analysis for the Washington State data (Garthwright, 1995). The FDA-BAM method has a method error variance of 0.35. In the DePaola et al. study (1990) the HGMF procedure as developed by Watkins et al. (1976) and later revised by Entis and Boleszczuk (1983) was used. When all suspect colonies are tested for confirmation, the precision of the hydrophobic grid membrane filtration (HGMF) procedure has been shown to be somewhat greater than the 3 tube MPN (most probable number) procedure (Entis and Boleszczuk, 1983; Watkins et al., 1976). In the DePaola et al. study (1990), enumeration of V. parahaemolyticus colonies was based on testing of five suspect colonies. Consequently, enumeration was not as precise as possible and overall method error associated with estimating V. parahaemolyticus densities may have been more comparable to that of a 3 tube MPN procedure. Therefore the method error variance of the FDA-BAM was considered a reasonable estimate of the method error for the HGMF method used in the DePaola et al. study (1990).

Analysis of Effects Other than Water Temperature on Mean log10 Vibrio parahaemolyticus per Gram
The effect of salinity was considered as a potential predictor of V. parahaemolyticus densities in addition to water temperature. For the three datasets considered here, a regression model that is linear in the effect of temperature and quadratic in the effect of salinity was fit to estimate the additional effect of salinity. This regression model is of the form:

log10(V.parahaemolyticus/g) = α + β *WTEMP + γ1*SAL+ γ2*SAL2 + ε

where WTEMP denotes water temperature in °C and SAL denotes salinity in parts per thousand (ppt). The parameters α and β are the regression parameters for the temperature effect, γ1 and γ2 are parameters for the salinity effect, and e is a random normal deviate with zero mean and variance σ2 corresponding to the combined effects of population and method error variation. The maximum likelihood estimates for the fit of this model to the data are shown in Table A5-3.

Table A5-3. Maximum Likelihood Estimates of Tobit Regression of log10 Vibrio parahaemolyticus per gram Versus Water Temperature and Salinity (MLEs and 95% Confidence Interval of Parameters of Temperature and Salinity Regression)
Study α β γ1 γ2 σ
DePaolaet al.,1990 -2.63(-2.14,0.08) 0.12(0.075,0.17) 0.18(0.016,0.34) -0.0042(-0.0084,0) 1.00(0.78,1.28)
FDA/ISSC, 2001 -2.05(-2.76,-1.34) 0.097(0.087,0.11) 0.20(0.13,0.27) -0.0055(-0.0073,-0.0038) 0.73(0.68,0.79)
Washington State DOH, 2000; 2001 -1.02(-34.3,35.0) 0.30(0.18,0.42) -0.39(-3.0,2.2) 0.0084(-0.04,0.06) 0.87(0.64,1.16)

For the DePaola et al. (1990) data set the effect of temperature is highly significant (p<0.0001) and the effect of salinity is marginally significant (p=0.03 for the linear term and p=0.05 for the quadratic term). The MLE of the optimal salinity level (-γ1/(2*γ2)) based on this model and data set is 21.4 ppt.

For the FDA/ISSC (2001) data set, the effects of both temperature and salinity were highly significant (p<0.0001). This is a consequence of the much larger sample size of this study compared to that of DePaola et al. (1990). The MLE of the optimal salinity level was 18.1 ppt. For the Washington State data, the effect of water temperature was also highly significant (p<0.0001) but salinity was not a significant effect. The range of salinities associated with samples in this data set was much narrower compared to the other two data sets and this would appear to be the most obvious reason for lack of significance.

The added value of prediction based on salinity as well as temperature is shown in Figure A5-5 as estimated by both the temperature-only and temperature/salinity regression fits to the FDA/ISSC (2001) data. The relative difference plotted in Figure A5-5 is the difference in the prediction based on temperature and salinity versus temperature-only divided by the prediction based on temperature alone. A relative difference greater than zero indicates that predictions based on water temperature and salinity are higher than based on water temperature alone. When salinity is in a nominal range of 10 to 25 ppt and water temperature is high (>25 °C), the relative difference in predicted values of mean log10 V. parahaemolyticus density is relatively small (i.e., an absolute difference of <10%). However, when water temperature is low the difference in predictions is more substantial (up to 40% at 15 °C). Water salinities in harvest areas are more variable than water temperature and no sufficiently comprehensive data sources were identified with respect to including this as a predictive factor in the assessment. Figure A5-5 indicates that the effect on model predictions of neglecting salinity effects is likely to be minor when water temperature is high (e.g., Gulf Coast summer). Furthermore, salinity may not be a strong effect in estuaries of the Pacific Northwest due to the fact that the range of salinities is narrower there than in other harvest areas. Thus, although salinity was identified as a significant effect in the regression analysis, its impact on predicted risk was not judged to be substantial and as such was not included as a parameter/component of the risk assessment model based on these considerations.

For the Washington State data, the possibility of additional effects such as year-to-year differences or differences between sampling areas were also considered. There was an apparent difference in the estimated regression relationship in 1998 when water temperatures were warmer than average but the difference was not large and regression parameter estimates obtained by fit to all available data were used in the assessment.

Image of Figure A5-5. Effect of Salinity on Predicted Mean log10 Vibrio parahaemolyticus Density in Oysters Relative to Predicted Density at Optimal Salinity (22 ppt).

Figure A5-5. Effect of Salinity on Predicted Mean log10 Vibrio parahaemolyticus Density in Oysters Relative to Predicted Density at Optimal Salinity (22 ppt).

Percentage of Total Vibrio parahaemolyticus that are Pathogenic
Studies of the distribution of total and pathogenic V. parahaemolyticus at harvest (Kaufman et al., 2003; DePaola et al., 2002, DePaola et al., 2003a) provide the best information available on both the mean and the variation of the relative abundance of pathogenic (tdh+) strains across samples. The information available from older studies is less detailed because the proportion of pathogenic strains in individual samples was generally not reported. Only the average percentage pathogenic, aggregated across multiple samples, can be inferred from the data reported in the older studies that were identified.

The study by DePaola et al. (2003a) was a component of the collaborative ISSC/FDA V. parahaemolyticus harvest study (FDA/ISSC, 2001). This study collected a total of 156 oyster samples from two harvest areas in Alabama over a period of 14 months in 1999 and 2000. The study by Kaufman et al. (2003) was conducted in the summer of 2001 with samples taken from a single Gulf Coast harvest area. A total of 60 individual oysters were sampled with half of these being analyzed immediately after harvest and the other half analyzed after 24 hours of storage at 26 °C. The DePaola et al., (2002) study analyzed samples from selected areas of Hood Canal collected in August 2001. Approximately 60 samples were analyzed for pathogenic and total V. parahaemolyticus in this study. All three studies utilized a direct plating procedure and gene probes to obtain paired counts of both pathogenic (tdh+) and total V. parahaemolyticus in analytical portions of each sample. The paired analytical portions assayed for pathogenic versus total V. parahaemolyticus were not necessarily the same volume (or weight).

The Beta-Binomial model was assumed as a model for the distribution of observed counts of pathogenic V. parahaemolyticus in analytical sample portions. This model implies an overall average percentage pathogenic but the percentage pathogenic in individual samples is also assumed to vary about the average according to a Beta distribution. The extent of the variation about the average that would occur in samples of oysters is likely dependent upon the size of the sample. In the context of the risk assessment, it is the variation of percentage pathogenic between servings that is of particular interest. The serving size varies with 6 and 12 oysters being typical. Environmental studies typically composite oysters, with 6 or 12 oysters per composite for microbiological analysis. Consequently, there is reasonable agreement as to definition of sample versus serving and it was judged that there was little need to correct the distribution of percentage pathogenic on the basis of a substantial difference between the number of oysters per "sample" versus per "serving". Implicitly, it is assumed that the oysters in a typical serving are harvested from the same location (i.e., come from the same harvest collection).

The variation of percentage pathogenic across samples was not estimated as the distribution of the ratio x/n, where x is the observed count of pathogenic V. parahaemolyticus (i.e., in an analytical sample portion) and n is a corresponding observed or estimated count of total V. parahaemolyticus (in a comparable volume) from the same oyster sample because, for most samples, the pathogenic count was zero. Estimation based on summarizing the data in this manner would potentially bias the estimate of percentage pathogenic distribution towards zero. Instead the Beta-Binomial distribution was fit directly to model the counts of pathogenic in the sample data and obtain estimates of the parameters of the Beta distribution assumed for variation of percentage pathogenic across samples.

Estimation of the distribution of percentage pathogenic was obtained conditional on the observed number of total V. parahaemolyticus in each sample and the volumes of sample (i.e., size of analytical portions) examined for both pathogenic and total, respectively. Specifically, given an observation of "n" total V. parahaemolyticus in an analytical portion of a sample of volume "Vt" and a count of "x" pathogenic V. parahaemolyticus in a replicate analytical portion of volume "Vp", it was assumed that the pathogenic count corresponded to the outcome of a (random) Binomial trial of size "integer(n*(Vp / Vt))" with probability of success equal to the percentage pathogenic in the sample. That is, the density of total V. parahaemolyticus in the sample portion examined for pathogenic was assumed known and equal to the estimate obtained from the sample portion assayed for total V. parahaemolyticus. The variability of total V. parahaemolyticus across analytical portions taken from the same sample was not considered for the purpose of estimation of percentage pathogenic.

The SAS procedure NLP (NonLinear Programming) was used to obtain the MLEs of the Beta-Binomial model for each of the data sets considered. Under the Beta-Binomial model the likelihood of the data is:

Image of the Beta-Binomial model equation

where Fp is the assumed Beta distribution (with parameters α and β ) for variation of the relative abundance (p) of pathogenic strains in different samples.

Reparameterization facilitated a more numerically stable estimation of the parameters of the Beta-Binomial model by the procedure NLP (e.g., to mitigate effects of excessively large values in the Gamma functions). The reparameterized likelihood used to obtain parameter estimates is of the form:

Image of the Gamma function equation

where P is the average percentage pathogenic and θ is a transformation of the overdispersion parameter (φ) of the Beta-Binomial:

P = α
φ =
1
θ =
φ
α + β α + β + 1 1 + φ

The MLEs for P and θ were obtained by NLP subject to the constraints that 0<P<1 and 0<θ<1. The corresponding MLEs for the original parameters α and β were then obtained by the inverse of the defining transformations of the reparameterization. The estimates obtained are shown in Table A5-4.

Table A5-4. Estimates of the Distribution of Percentage Pathogenic Vibrio parahaemolyticus
Study Data Pa φ α β 95% Confidence Interval for P
DePaola et al., 2002 (Hood Canal)b 2.33% 0.076 0.283 11.86 (1.05%, 5.47%)
Kaufman et al., 2003 (0 hr) 0.51% 0.0114 0.442 86.2 (0.21%, 1.47%)
Kaufman et al., 2003 (24 hr) 0.08% 0.0011 0.681 907 (0.03%, 0.16%)
Kaufman et al., 2003 (all data)c 0.18% 0.0045 0.394 221 (0.09%, 0.44%)
DePaola et al., 2003a 0.44% 0.0146 0.297 67.2 (0.24%, 0.82%)

a α and β denote the parameters, φ denotes the overdispersion and P denotes the average of the assumed Beta distribution
bestimate used in the risk assessment model for the Pacific Northwest region
cestimate used in the risk assessment model for regions outside of the Pacific Northwest

Given the discrete nature of the observed count data, and small samples sizes relative to the mean percentage pathogenic being estimated, a parametric bootstrap procedure (Garren et al., 2001) was used to estimate an uncertainty distribution for the parameters α and β. The parametric model assumed was the Beta-Binomial with the parameter values equal to the MLEs obtained based on the observed data. With respect to each study, replicate bootstrap samples of the count of tdh+ colonies in analytical sample portions were generated by random sampling of percentage pathogenic from the fitted Beta distributions followed by random sampling from Binomial distributions with size parameter equal to the number of total V. parahaemolyticus observed (or estimated) in a volume of sample comparable to that assayed for tdh+ colonies. A total of 1,000 bootstrapped data sets were generated for each study and MLEs for the parameters α and β were obtained for each bootstrap by the NLP procedure as described above. On a few rare occasions the NLP procedure did not converge for a bootstrapped outcome, suggesting the lack of existence of an MLE. When this occurred the bootstrapped outcome was dropped from the set defining uncertainty in the estimates of α and β.

All model fits of the Beta-Poisson to the observed data were found to be adequate. Goodness-of-fit was assessed by the method of Brooks et al. (1997). Briefly, the maximum likelihood of the fit of the Beta-Binomial model to the observed data (of each study separately) was compared to a null distribution of maximum likelihood values generated by parametric bootstrapping of hypothetical outcomes and obtaining maximum likelihood values for each bootstrap by refitting the Beta-Binomial model. If the maximum likelihood value of the fit of the Beta-Binomial model to the observed data lies at either extreme of the null distribution obtained by bootstrapping then the fit is questionable. As shown in Table A5-5, all of the maximum likelihood values obtained with respect to the observed data are not extreme (e.g., not < 2.5%-tile or > 97.5%-tile). The nature of the variation of relative abundance of pathogenic strains could be other than Beta-Binomial, but the fit of this model to the datasets was not rejected or marginal in any way. The hypothesis that there is no extra-binomial variation in the data (i.e., fit of a binomial model to the data would be adequate) was rejected at the 95% confidence level for all three data sets based on the deviance statistic of the fits.

Table A5-5. Goodness of Fit of Beta-Binomial Model
Study Data Log10 of Maximum Likelihood of Fit Percentile of the Null Distribution a
DePaola et al., 2002 (Hood Canal) -34.35 0.44
Kaufman et al., 2003 (0 hr) -24.50 0.27
Kaufman et al., 2003 (24 hr) -25.63 0.59
Kaufman et al., 2003(all data) -55.77 0.38
DePaola et al., 2003a -62.94 0.13

athe null distribution of maximum likelihood values was estimated by 1,000 Monte Carlo samples of hypothetical (bootstrap) outcomes based on the MLE of the parameters to the observed data

Of the three data sets analyzed, the Hood Canal study (DePaola et al., 2002) was the only study with samples taken from the Pacific Northwest. The results of model fits to this data set were taken to be representative of the Pacific Northwest region. Estimates based on the pooled 0 and 24-hour time points of the Kaufman et al. (2003) study were used to model all other areas of the country. Collectively, the estimate of the mean based on combined 0 and 24-hour time points of the Kaufman et al. study (2003) was lower than that of the DePaola et al. (2003a) study, but not significantly so. Furthermore, although the 0 versus 24-hour time points suggested differences in the percentage pathogenic distribution, an estimated mean of 0.18% based on pooling of the data was used for the purpose of risk assessment. Previous studies of the percentage of V. parahaemolyticus isolates that are pathogenic in retail samples (FDA/ISSC, 2000) and studies of the growth rate of pathogenic versus nonpathogenic strains (Cook, 2002a) do not support the hypothesis that there is any appreciable difference in percentage pathogenic at retail versus at harvest levels. The inferred uncertainty distributions of mean percentage pathogenic and the underlying bootstrap uncertainty distributions for a and b are shown in Figures A5-6 and A5-7, respectively.

Image of Figure A5-6. Histograms of Bootstrap Uncertainty Distributions of Mean Percentage Pathogenic Based on the Beta Distribution Model of Sample-to-Sample Variation of Percentage Pathogenic Vibrio parahaemolyticus (Pacific Northwest and Non-Pacific Coastal regions).

Figure A5-6. Histograms of Bootstrap Uncertainty Distributions of Mean Percentage Pathogenic Based on the Beta Distribution Model of Sample-to-Sample Variation of Percentage Pathogenic Vibrio parahaemolyticus (Pacific Northwest and Non-Pacific Coastal regions).


Image of Figure A5-7. Bootstrap Samples of Uncertainty of Mean Percentage Pathogenic and Dispersion Parameter of the Beta Distribution Model of Sample-to-Sample Variation of Percentage Pathogenic Vibrio parahaemolyticus (Pacific Northwest and non-Pacific Coastal regions).

Figure A5-7. Bootstrap Samples of Uncertainty of Mean Percentage Pathogenic and Dispersion Parameter of the Beta Distribution Model of Sample-to-Sample Variation of Percentage Pathogenic Vibrio parahaemolyticus (Pacific Northwest and non-Pacific Coastal regions).

Post Harvest Module

Based on the distribution of V. parahaemolyticus/g at harvest, the Post-Harvest Module predicts the distribution of V. parahaemolyticus/g at the time of consumption by modeling the influence of various factors on the outgrowth and die-off likely to occur during storage. Principally, this is determined by distributions of time that oysters are exposed to ambient temperature during harvest before first refrigeration, various times of transport and storage together with distributions of temperature of exposure and rates of growth or decline versus temperature.

Distribution of Air Temperature at Oyster Harvest
Oysters are typically harvested from shallow water estuaries (e.g. 1-3 m depth). Consequent to the action of wind and tides, one would expect a correlation between the day-to-day variations of the water temperature and that of the air temperatures. This was confirmed from examination of the data from the NBDC buoys selected as being representative of the four harvest regions.

To facilitate the Monte Carlo simulation of the risk assessment model, the correlation between the air and water temperatures that oysters are exposed to was implemented by using the distribution of the difference between air and water temperature as a model parameter. This difference and water temperature were then used as inputs to determine the distribution of air temperatures that oysters are exposed to during harvesting. This approach was adopted in order to assure proper correlation between water and air temperatures.

Analysis of the NBDC (1997) data revealed that the relationship between air and water temperature changes during the course of the day. This is due to the fact that the temperature of air is more variable over a 24-hour period than that of the water. Figure A5-8 shows the mean difference between air and water temperature as a function of the time of day at the Dauphin Island Buoy (1997 data). The relationship between air and water temperatures in the Gulf are different from more northern areas of the country in that the mean water temperature is always warmer than that of the air. In northern areas of the country, i.e., the other 3 harvest regions in the assessment, the mean water temperature is cooler than that of the air during the spring and summer, and the reverse is true during the fall and winter.

Image of Figure A5-8. Variation of mean hourly air-water temperature differences for Dauphin Island, Alabama Buoy (NBDC, 1997).

Figure A5-8. Variation of mean hourly air-water temperature differences for Dauphin Island, Alabama Buoy (NBDC, 1997).


Image of Figure A5-9. Correlation of Daily Midday Air and Water Temperature Measurements for Dauphin Island, AL (NBDC, 1997).

Figure A5-9. Correlation of Daily Midday Air and Water Temperature Measurements for Dauphin Island, AL (NBDC, 1997).

To minimize complexity, the relationship of air versus water temperature at 12:00 pm (midday) was assumed to be a reasonable average for the purpose of estimating a distribution of temperatures that oysters are typically exposed to during the course of oyster harvesting (i.e., harvesting starts in the early to mid morning and may last through mid-afternoon). Figure A5-9 shows the correlation of air versus water temperatures measured at midday at the Dauphin Island Buoy in 1997. Figure A5-10 shows the distribution of the difference between the same air and water temperature data as a function of water temperature. Clearly, the difference between air and water temperature is more variable when the water temperature is lower.

There were no noticeable differences in this relationship between air versus water temperatures within the same season across different years, either for Dauphin Island or any of the other NBDC sites utilized for the assessment. Consequently, seasonal distributions of the difference between air and water temperature were obtained by pooling all available years of data. The distribution of the difference in air temperature versus water temperature was then approximated as being Gaussian within each region/season classification, with mean and standard deviation estimated by the method of moments. As was the case with water temperature, this summary distribution is only an approximation since the air temperatures in the NBDC data exhibit the same degree of skewness as was discussed above with regard to water temperatures (Figure A5-3). The air-water temperature difference is also slightly skewed but less than that of either air or water temperature alone.

Distribution of Time to 1st Refrigeration

The distribution of time that oysters are exposed to ambient air temperatures during harvest (i.e., prior to refrigeration) was derived based on the distribution of duration of harvest and a distribution for the time when individual oyster lots are collected during the harvest. For the first distribution, the only information identified was minimum, maximum and most likely durations obtained by interviews with harvesters in several Gulf Coast states (GCSL, 1997). Based on this information estimated distributions of duration of harvest were taken to be Beta-PERT distributions with specified minimums, maximums and modes for each region and season combination.

The relative proportion of the harvest caught during the course of harvesting operations may vary somewhat from one harvest area to the next. However, with the exception of time required to return to dock from the harvest area, a constant harvesting operation was assumed to be typical of the majority of harvest areas. A time of 1 hour was considered typical of time to return to dock from the harvest areas. Constant harvesting operation implies that the distribution of the catch within the harvest period is uniform. Thus the time of collection of oyster lots, relative to time of 1st refrigeration, was taken to be a continuous uniform random variable with minimum time equal to 1 hour and maximum equal to the duration of harvesting operation. Distribution of time to 1st refrigeration for individual lots is the mixture of the distribution functions for duration of harvesting operation and the distribution of time of collection for a given duration of harvesting.

Image of Figure A5-10. Differences Between Midday Air and Water Temperatures as a Function of Water Temperature for Dauphin Island, AL (NBDC, 1997).

Figure A5-10. Differences Between Midday Air and Water Temperatures as a Function of Water Temperature for Dauphin Island, AL (NBDC, 1997).

Growth Rates

The growth and die-off rates were based on estimates obtained from the published literature (Miles et al., 1997, Gooch et al., 2002, FDA/ISSC, 2000). Statistical criteria used to obtain these estimates are described in the respective references. The growth rate study in oysters (Gooch et al., 2002) was conducted at only one temperature (26 °C). Therefore, a study of growth rate in broth (axenic) culture (Miles et al., 1997) across a range of temperatures (and water activities) was used as a basis to extrapolate growth rate in oysters to temperatures higher and lower than 26 °C. The uncertainty associated with this prediction was addressed in the assessment by incorporating an uncertainty distribution for the relative growth rate in broth culture versus in oysters. The uncertainty distribution selected for this ratio was taken to be a triangle distribution with minimum of 3, mode of 4 and maximum of 5. The mode of 4 corresponds to the best estimate of the ratio of the predicted growth rate in broth culture at 26 °C, using the Miles et al. (1997) model, compared to the growth actually observed in oysters held at 26 °C (Gooch et al., 2002).

More specifically, with respect to the growth rate, Miles et al. (1997) obtained worse case estimates based on the fastest growing of four strains that were studied. For each combination of temperature and water activity, the extent of bacterial growth observed was modeled using the Gompertz function and an estimate of the maximal rate of growth was obtained. A secondary model was then used to estimate the effect of environmental parameters (temperature and water activity) on the maximal growth rate. The model that was assumed by Miles et al. (1997) was of the square root type:

Image of the square root type equation

where

mm = maximal growth rate (log10 per minute)
aw = water activity
T = temperature (in degree Kelvin)

The estimates of the parameters that were obtained are:

b = 0.0356
c = 0.34
Tmin = 278.5
Tmax= = 319.6
aw,min = 0.921
aw,max = 0.998
d = 263.64

The parameters Tmin, Tmax, aw,min, and aw,max denote the range of temperatures and water activity over which growth can occur. The authors validated their model by comparison of model predictions with observed rates in eight other studies of growth in broth model systems obtained from the literature.

To use the (1997) et al. equation as a prediction of growth rate in oysters it was assumed that water activity of oysters does not vary substantially with a nominal value equal to the optimal value of 0.985 predicted to occur under broth culture conditions. At this water activity, the predicted growth rate in broth at 26 oC (78.8 °F) is 0.84 log10 per hour, which is approximately a 7-fold increase in density per hour. This is four times greater than the rate of growth observed for V. parahaemolyticus in oysters held at 26 oC (78.8 °F) (Gooch et al., 1999). Therefore, based on this observation, prediction of the growth rate in oysters at temperatures other than 26 °C (78.8 °F) was obtained by dividing the predicted rate for broth culture by a factor of four. This assumes that the growth rate in oysters is a constant fraction of the growth rate in broth at all temperatures. The influence of this assumption in the risk assessment was evaluated by considering this factor as an uncertainty parameter varying according to a triangle distribution in the range of 3 to 5 with a mean of 4. This gives an indication of the sensitivity of our conclusions to the magnitude of the relative growth rate in oysters versus broth culture but does not fully address the uncertainty in so far as it is conceivable that the relative growth rate could be temperature dependent. Although the appropriateness of the assumption has not been fully validated, the ambient temperature of Gulf Coast is close to 26 °C (78.8 °F) from April through October and this is a region and season for which the largest number of V. parahaemolyticus cases is associated.

A plot of the resulting model prediction for mm as a function of either temperature or water activity is a unimodal function with a maximum value and zero growth rate outside of the predicted range of temperatures and water activity favorable for growth. To use this equation as a prediction of growth rate in oysters it was assumed that water activity of oysters does not vary substantially with a nominal value equal to the optimal value of 0.985 predicted to occur under broth culture conditions. At this water activity, the predicted growth rate in broth at 26 oC (78.8 °F) is 0.84 log10 per hour, which is approximately a 7-fold increase in density per hour. This is four times greater than the rate of growth observed for V. parahaemolyticus in oysters held at 26 oC (78.8 °F) (Gooch et al., 1999; Gooch et al., 2002).

Oyster meat temperature
Air temperature was used as a surrogate for oyster meat temperature for oysters harvested by dredging and intertidal. For oysters harvested in intertidal areas, additional growth of V. parahaemolyticus was considered as described below.

Effect of Intertidal Exposure in the Pacific Northwest
Unlike other areas of the country, a significant fraction of oysters harvested in the Pacific Northwest are collected when oyster reefs are exposed during the course of the tide cycle. Exposure to the air and consequent radiative heating of oysters in bright sunlight can elevate oyster temperatures substantially above that of the water temperature. To model the effect of intertidal harvesting on V. parahaemolyticus densities in the Pacific Northwest, a distribution of oyster temperature during intertidal exposure was developed based on the observational data available and an estimate of the proportion of days that are subject to cloudy, partly cloudy or sunny conditions. Radiative heating, leading to oyster temperatures well above ambient air temperature was considered likely for sunny days and unlikely for cloudy days.

Estimates of the fraction of the Pacific Northwest catch that are harvested during intertidal cycles were obtained based on data for harvest volume from selected areas of Washington State. This information was combined with expert opinion concerning the fraction of harvest from each area that is collected intertidally rather than dredged (i.e., from submerged reefs) (Kaysner, 2002). The harvest volumes of selected areas of Washington State are shown in Table A5-6.

Table A5-6. Average Total (Dredged and Intertidally Picked) Oyster Shellfish Harvest (in pounds) in Selected Areas of Washington State by Season (Yearly Averages Based on 1990-2001 data)
Area Winter Spring Summer Fall
Hood Canal 389,000 480,000 378,000 416,000
North Sound 328,000 308,000 254,000 245,000
South Sound 844,000 574,000 437,000 595,000
Willapa Bay/Grays Harbor 324,000 259,000 205,000 198,000
Total 1,886,000 1,620,000 1,274,000 1,454,000

Note: Harvest intended for the shucked market was excluded since this market is typically not intended for raw consumption


Table A5-7. Average Area-Specific Oyster Shellfish Production Expressed as a Percentage of the Total Oyster Shellfish Harvest (in pounds) for Selected Areas of Washington State by Season (yearly averages based on 1990-2001 data)
Area Winter Spring Summer Fall
Hood Canal 20.6% 29.6% 29.7% 28.6%
North Sound 17.4% 19.0% 19.9% 16.8%
South Sound 44.8% 35.4% 34.3% 40.9%
Willapa Bay/Grays Harbor 17.2% 16.0% 16.1% 13.6%

Table A5-8. Percentage of Area-Specific Oyster Shellfish Harvest Collected During Intertidal Exposure
Area Percentage Intertidal Percentage Dredged
North Sound 75% 25%
Hood Canal 95% 5%
South Sound 90% 10%
Willapa Bay/Grays Harbor 10% 90%

An estimate of the overall percentage of the harvest collected intertidally in Washington State for each season was obtained by weighting the expert opinion on area-specific percentages for intertidal collection (Table A5-8) by the percentages that each area contributes to the total shellfish harvest (Table A5-7). The shellfish harvest excludes the portion of the total harvest intended for the shucked market. The harvest statistics indicated that virtually all oysters intended for the shucked market were harvested by dredging. In these calculations it was assumed that the area-specific percentages for intertidal versus dredged harvest (Table A5-8) do not vary by season. Thus, the seasonal fraction of the total oyster harvest collected intertidally was calculated as:

% intertidal 1 = ∑ihi*pi

where hi is the percentage of total harvest collected in the ith area and pi is the percentage of the harvest in that area collected during intertidal exposure. The results of these calculations give an estimate of 75% of the total harvest being collected intertidally. There was very little variation in this estimated percentage across different seasons; therefore the aggregate average was used as an estimate for all seasons in the risk assessment model. The remaining 25% of the harvest, being dredged, was not subject to predictions of growth that may occur due to elevated oyster temperatures. The Post-Harvest growth of V. parahaemolyticus in this latter portion of the harvest was treated in a manner similar to that applied to the other 3 regions of the country.

Studies of V. parahaemolyticus densities in Washington State (DePaola et al., 2002; Herwig and Cheney, 2001) provide some observational data bearing on the extent to which oyster temperatures are elevated during intertidal exposure. In the DePaola et al. (2002) study, a total of 17 temperature measurements were taken over a period of a week in conjunction with the microbiological analysis of oyster samples collected at the end of the exposure cycle (i.e., after full or "maximum" exposure to ambient air temperatures and sunlight). Across this set of measurements, the minimum, maximum, and mean oyster temperatures were 23.3, 32.6, and 27.5 °C, respectively. The distribution of temperatures was almost uniform over the range from minimum to maximum with the 25%-tile and 75%-tile of the distribution being 25 and 30 °C, respectively. Compared to air temperatures, oyster temperatures after maximum exposure were an average of 5 °C (9 °F) greater than air temperature. The maximum difference was 8 °C (14.4 °F). The mean oyster temperature was equal to the median and thus was not a consequence of any extreme observations. At the time of the study the water temperature in the Hood Canal estuary was slightly less than 17 °C. Elevated oyster temperatures, relative to that of the air were also reported by Herwig and Cheney (2001).

National Weather Service (NWS) historical data indicate that during the summer, in the Pacific Northwest, meteorological conditions are evenly divided between cloudy, partly cloudy and sunny conditions. A higher proportion of cloudy days occur during the winter but, given that summer is the higher risk season, the proportion of sunny versus cloudy days during the summer was considered to be more pertinent. Based on this information and the range of oyster temperatures observed in the study by DePaola et al. (2002), average difference between oyster temperatures versus air during intertidal exposure was modeled as being uniform in the range of 0 to 10 °C. The duration of exposure to ambient air and radiative heat was assumed to be uniform in the range of 4 to 8 hours; e.g., in consideration of the likely variation in the depth of oyster beds in relation to tide height and flows. The duration of oyster harvesting for intertidal was assumed to commence at the start of oyster collection. Given the estimate of a minimum of 2, mean of 8 and maximum of 11 hours for duration of harvest for the Pacific Northwest (i.e., in regard to ISSC time-to-refrigeration guidelines, ISSC&FDA, 1997), it was assumed that oysters harvested intertidally would reach refrigeration in a maximum of 11 hours from the start of collection. The duration of transport time (in hours) after intertidal exposure was therefore taken as:

Max(Beta-PERT(2,8,11) - Intertidal Exposure Time, 1)

Growth during the period of transport was assumed to occur at a rate commensurate with air temperature, as oysters are typically collected by boat or barge after being briefly cooled by water when the tide comes in and the oysters are retrieved for transport to processing facilities. A minimum transport time of 1 hour was assumed.

Distribution of Time to Reach No-growth Temperatures and Duration of Cold Storage
There is little data available to precisely quantify the distribution of the length of time required for oyster lots to reach no-growth temperatures after being placed in cold storage after transport from the harvest areas. Therefore it has been assumed that the distribution is uniform between 1 and 10 hours. This distribution was chosen to represent a mixture of both variability and uncertainty. A maximum time of 10 hours was selected based on literature of cooling studies with other food products, primarily meat. Studies of oyster equilibration to air temperature were undertaken by GCSL, in part to validate the reasonableness of this distribution (Cook, 2002b). The component of these temperature studies pertinent to the distribution of time to reach no growth temperature during 1st refrigeration was conducted with initial oyster temperatures of 25 °C and cooler temperature of 4 °C. In this study the temperatures of selected oysters within a "sack" were continuously recorded during cooldown. Oysters located toward the center of the sack/lot cooled more slowly than those on the outside. Temperatures of individual oysters decreased exponentially, reaching the cooler temperature at times ranging from 7 to 8.5 hours. This was considered confirmatory of the maximum of 10 hours assumed for the distribution given that the loading of commercial coolers is likely to be heavier and more variable than that typified by conditions in the experimental study. Furthermore, all commercial coolers may not be consistently at or below 4 °C.

Consumption Module

The distribution of the dose of pathogenic V. parahaemolyticus ingested per serving was estimated based on combining distributions of (a) the number of oysters consumed; (b) the weight of oysters consumed; and (c) the density of total V. parahaemolyticus per g and (d) the percentage of total V. parahaemolyticus per g that are pathogenic. Estimated distributions of the number of oysters consumed and the weight of oysters consumed is addressed here.

Number of Oysters per Serving
The modeled distribution of the number of oysters per serving was taken to be equal to the empirical distribution observed in response to a consumer survey conducted in 1994 by the Florida Agricultural Market Research Center (Degner and Petrone, 1994). In this study, the average number of oysters eaten per occasion was reported to be 13.8 with a range of 1 to 60. Given the relatively large number of respondents in the survey (n=319) and the evident multimodal characteristics of this distribution (6, 12, and 24 oysters/serving being the most probable), the empirical distribution was taken as an estimate rather than attempting to summarize the data by fit of a parametric distribution.

The survey was conducted in a coastal area, where consumption of oysters per eating occasion may be expected to be higher than in inland areas of the country. However, no other suitable sources of data were identified with respect to consumption patterns nationwide and it was judged that the potential bias in using the distribution as a nationwide estimate was minimal in comparison to other modeling uncertainties impacting estimated dose per serving.

Distribution of Meat Weight per Oyster
Given a distribution of the number of oysters per serving, an estimate of meat weight per oyster is needed to determine a distribution of the meat weight consumed per serving. The most relevant data identified to estimate the gram weight of oysters was the ISSC/FDA retail data (FDA/ISSC, 2000; DePaola, 2002). In this study, 339 of the 370 oyster samples collected from wholesale and retail locations were weighed prior to microbiological analysis. Samples generally consisted of composites of 12 oysters (range, 4-15) and this included both the oyster meat and the mantle fluid. The average oyster (i.e., meat and mantle fluid) weight per sample was calculated by dividing the total gram weight of the composite sample by the number of oysters in the sample. The resulting distribution of average oyster weight per sample was found to be positively skewed. The distribution is shown in Figure A5-11.

Image of Figure A5-11. Distribution of Average Oyster (Meat and Mantle Fluid) Weight Over Samples of Composites of 4-15 Oysters Collected From Retail Establishments (FDA/ISSC, 2000; DePaola, 2002).

Figure A5-11. Distribution of Average Oyster (Meat and Mantle Fluid) Weight Over Samples of Composites of 4-15 Oysters Collected From Retail Establishments (FDA/ISSC, 2000; DePaola, 2002).

Although there were some apparent differences in the mean oyster weight distribution by region and season of harvest, the differences were not large. A single estimate of the distribution of average gram weight per oyster based on pooling all of the data was considered appropriate and this estimate was assumed to apply to oysters harvested from all regions and seasons. A lognormal distribution was fit to the observed average oyster weight data in order to obtain a smooth estimate of the average oyster weight, rather than using the empirical distribution of the data. The maximum likelihood estimates obtained corresponded to a geometric mean average oyster weight per sample of 15.2 grams and a geometric standard deviation of 1.4 grams (Figure A5-11).

Samples in the retail study consisted of composites of both oyster meat and mantle fluid. Accordingly, a correction was applied to infer the average meat weight per oyster consumed. Oyster mantle fluid is typically not consumed with the oyster meat. The distribution of the ratio of meat weight to total (meat and mantle fluid) oyster weight based on measurements of individual Gulf Coast oysters collected during the Kaufman et al. (2003) study is shown in Figure A5-12. Although there is a distribution of percentage meat weight per oyster the coefficient of variation is very small. The mean of the distribution is 90%. Given the relatively small coefficient of variation, an average percentage was used, rather than the distribution, to determine a distribution of oyster meat weight consumed from the distribution of oyster weights shown in Figure A5-11.

Image of Figure A5-12. Distribution of Oyster Meat Weight as the Percentage of Total Oyster (Meat and Mantle Fluid) Weight Over Samples of Individual Oysters (Kaufman et al., 2003).

Figure A5-12. Distribution of Oyster Meat Weight as the Percentage of Total Oyster (Meat and Mantle Fluid) Weight Over Samples of Individual Oysters (Kaufman et al., 2003).


Table of Contents

horizontal rule
horizontal rule