Estimation of time-dependent arbovirus infection risk in blood and tissue donations

Brad J. Biggerstaff

doi:10.1515/scid-2020-0001

Article Open Access

Estimation of time-dependent arbovirus infection risk in blood and tissue donations

Brad J. Biggerstaff

Published/Copyright: November 13, 2020

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Statistical Communications in Infectious Diseases Volume 12 Issue 1

Abstract

West Nile virus (WNV) outbreaks raise the concern of WNV infection in donated blood and blood products destined for transfusion. We describe methods we developed to estimate time-dependent risk of WNV infection in donated blood, including improvements not previously detailed. The methods are then extended for use in estimation of the risk of WNV infection in donated cadaveric tissues by introducing stratification and stratum-specific weighting to address novel aspects of this application. Data from the WNV outbreak in Colorado in 2003 are used to estimate risk for donated cardiac tissue.

Keywords: blood safety; Poisson process; public health; tissue safety

Introduction

Relatively soon after the documentation in 1999 of West Nile virus (WNV) circulation in North America, it was recognized that during a WNV outbreak the probability or “risk” of WNV infection in blood donors may reach sufficiently high levels to pose a threat to the blood supply. In two papers, Biggerstaff and Petersen (2002, 2003) used surveillance data for WNV neuroinvasive disease (WNND; e.g., meningitis, encephalitis, acute flaccid paralysis) collected by local and state public health agencies and reported to the Centers for Disease Control and Prevention (CDC) along with exogenous information on WNV viral kinetics in blood, the inapparent-to-apparent infection ratio, and the overall rate of asymptomatic infection to estimate the time-dependent risk of donors being WNV-infected at donation in regions and during timeframes of WNV epidemic activity. In the flavivirus family, WNV is a mosquito-borne virus, infection with which may result in no illness (asymptomatic infection), a relatively mild febrile illness, and severe illness involving neurological processes that may lead to long-term incapacitation and death. More severe outcomes are less common, and severity tends to increase with age. By far the most common route of WNV transmission to humans is from the bite of an infected mosquito. Other modes of transmission have been documented, however, including laboratory exposure; from mother to baby during pregnancy, delivery, or breast feeding; and via blood transfusion or organ donation. Coincidentally, the first documented case of WNV transmission from donated, infected blood occurred in 2002, just weeks after the first estimates of such risk were published for the 1999 outbreak in Queens, New York City (Biggerstaff and Petersen 2002; Centers for Disease Control and Prevention 2002; Iwamoto et al. 2003).

Technical details of the initial statistical methods described here were used in the original paper for the Queens outbreak were presented in its appendix, with less formal description provided in the body of the text. Subsequent applications to other geographic areas, time scales, and viruses (Biggerstaff and Petersen 2003; Petersen, Tomashek, and Biggerstaff 2012; Shang et al. 2007, 2016) introduced improvements, extensions, and simplifications to the methods, though none of these papers contained detailed information on these changes from the original development. A more recent application (Biggerstaff and Petersen 2020) of the general methodological ideas to estimation of the risk of WNV infection in cadaveric tissue donations necessitated further methodological enhancements, and the present paper collects all of these updates into a single, detailed presentation.

We begin by providing background on the estimation goal and the information needed to reach it in Background & information required. Methodological and computational development then recapitulates and expands the appendix in the original paper for the Queens outbreak (Biggerstaff and Petersen 2002), with extensions to incorporate the asymptomatic proportion and uncertainty in this, as well as uncertainty in the inapparent-to-apparent infection ratio (Biggerstaff and Petersen 2003); and subsequently an extension to incorporate multi-year estimation (Petersen, Tomashek, and Biggerstaff 2012). Refinements for estimation of WNV in cadaveric tissue donations introduces new refinements to these methods to accommodate stratification and differential weighting by stratum to address specific characteristics of estimation for WNV infection in cadaveric tissue donations (Biggerstaff and Petersen 2020).

In Development using a Poisson process model for the underlying case onset times we detail how our resampling implementation can be viewed as assuming an inhomogeneous Poisson process model for the underlying WNND case onset times and estimating the associated cumulative rate function using nonparametric Gaussian kernel smoothing, which we indicated but did not detail in the Discussion of Biggerstaff and Petersen (2002). In Application: cardiac tissue, Colorado, 2003 we illustrate the updated methods on a real data set, the full analysis for which is reported in Biggerstaff and Petersen (2020). The paper closes with a summary of the computational resources used followed by a discussion section.

Background & information required

The recognition that the probability or risk of WNV transmission from donated blood may be unacceptably high during an active, regional outbreak of WNV was the motivation for the original applications discussed here. Because WNV activity is necessarily seasonal in North America, being a mosquito-borne virus, estimation methods were developed to reflect this time-dependence. In the United States, WNND cases are reported through established public health surveillance to the CDC with accompanying information, including the date that the patient first exhibited symptoms. Because WNND is a severe illness, surveillance systems capturing its occurrence in the US are thought to be essentially complete, and so epidemic curves developed using these data can be used to characterize the population-level transmission dynamics in the general population. These WNND case onset data, therefore, form the basis for the estimation methods used here. For the example from the WNV oubreak in Colorado, 2003, the WNND case onset data can be summarized in a epidemic curve, shown in Figure 1 for cases aged 13–55 for the application we present.

Figure 1:

Number of WNND cases aged 13–55 years by date of symptom onset, Colorado, 2003.

Mostashari et al. (2001) reported an estimate of the ratio of the number of WNV infections to the number of WNND cases in a population subject to WNV transmission. Since we have available the WNND case onset dates from surveillance, we can use this inapparent-to-apparent infection ratio to estimate the number of infections in the population. This ratio is the key parameter in relating the observed surveillance data to population-level infections. Original applications used a single inapparent-to-apparent infection ratio for the whole population; more recent work (Carson et al. 2012; Williamson et al. 2016) has refined these estimates by age group and sex, which we are, therefore, now able to incorporate into our analyses, as we show below.

When focusing specifically on estimation for WNV infection in blood and tissue donations, we restrict estimates of the number infections in the population in two ways. Because transmission of WNV occurs only when a donation contains or is infected with viable virus (termed viremic in blood), and because the duration of such infection is transient, estimates of when a potential donor was infected are necessary to characterize the risk. Using the WNND case onset times as anchor times, we use estimates of the duration from infection to symptom onset—the incubation period—and then estimates of the duration of infection to do this. Finally, because donors are necessarily asymptomatic—that is, they do not exhibit any symptoms of illness—and are screened for symptoms that would require deferral from either blood or tissue donation, we account for this using an estimate of the proportion of asymptomatic infections (Mostashari et al. 2001).

Methodological and computational development

We develop the methods considering specifically applications to estimation of WNV transmission risk via blood or tissue donation. The details, however, are written to be able to be applied more generally. The intent in the development, therefore, is that whenever densities or distributions or associated parameters are required to be specified, the user would specify such in the application under consideration. Whenever possible in application, we make distributional assumptions and specifications based on published data. For example, below we specify a Weibull distribution for the duration of WNV viremia based on data reported in Southam and Moore (1954), but in another application, a different distributional assumption for the duration of the pathogen’s duration in blood or tissues may be more appropriate. In such case, in the algorithm presented in Estimation algorithm below, the reader should utilize the reader’s own distributional model assumptions.

General approach, motivation, and estimation

For each i=1,2,…,n, let Y_i be observed WNND case onset times (recorded to the day) in chronological order, and set Y₁=0 and Y_n=T. Further, for each case i, let V_0i be the duration from time of infection to symptom onset, and let V_1i be the duration of infection. Assume that the timing and duration of infection are independent of symptom onset time, and that the V_0i are independent and identically distributed, as are the V_1i; finally assume that the V_0i are independent of the V_1i.

The case onset dates Y_i are recorded to the day, so to reflect the underlying continuous nature of the infection process in the environment and population we statistically smooth these observations by setting Xi=Yi+hεi, where h>0 is a smoothing parameter and the εi are independent standard normal deviates. We refer to “smoothing” here somewhat informally, though we note that in computation, generation of a value X_i by adding the scaled noise hεi for a generated Gaussian deviate ε_i is equivalent to sampling from a Gaussian kernel density smooth of the original {Yi}. Further, we note that specification of h is important and related to the problem of density estimation or Poisson intensity estimation, and we expand on this in Development using a Poisson process model for the underlying case onset times. We have used the value of Sheather and Jones (1991), though on occasions when this value for h was not computable due to insufficient data, we used the rule-of-thumb estimator given in Scott (1992). To preserve seasonality for multi-season data, we have estimated h using within-season information, for example by considering only day-of-year for multi-year observations of data with annual cycles, or by computing different h values for each year and averaging these or resampling from them.

Viewing the smoothed case onset times X_i as anchor times, the values Xi−V0i and Xi−V0i+V1i represent the initial and terminal times of infection, respectively, for case i. For a given time, t∈(0,T), then, case i is viremic at t whenever Xi−V0i<t<Xi−V0i+V1i. We count the number of cases viremic at t using the function

(1)Vt=∑i=1nIXi−V0i, Xi−V0i+V1it

where IA(x)=1 if x∈A and IA(x)=0 if x∉A.

We interpret a particular realization of this random function V(t) as representative of such realizations, where the random mechanisms are the environmental conditions and personal behaviors leading to exposure and subsequent advancement to infection in the population. The expected value EVt is the mean of such a theoretical population and is the basis for the quantity we want to estimate.

As we saw, the inapparent-to-apparent infection ratio, ρ, with estimate ρˆ, provides a means to inflate the observed counts of WNND cases to estimate the number of infections in the population. In particular, writing u₊ for the next integer greater than u, the quantity ρˆ+ V(t) estimates the number of individuals in the population who are viremic at time t. It will be convenient below to recast the quantity ρˆ+ V(t) as follows. Since V(t) counts the number of observed cases, V(t) is a non-negative integer at each t. (In addition to the desire to produce an integer result, the use of ρˆ+ is expected to provide a conservative estimate, since E[ρˆ+]≤E[ρˆ]+.) Then write ρˆ+V(t) as a simple expansion of the expression in Eq. (1),

ρˆ+ Vt=∑j=1ρˆ+∑i=1nIXi−V0i, Xi−V0i+V1it=∑j=1ρˆ+∑i=1nIXij−V0ij, Xij−V0ij+V1ijt

which is just replication of V(t) a total of ρˆ+ times, setting each Xi(j)=Xi, etc.

Because we are interested in estimating the probability of a donor being viremic, and since donors are necessarily asymptomatic, we randomly truncate infection times at symptom onset for the proportion of infections expected to develop symptoms. This gives us, finally, an estimate of the number of asymptomatic and WNV infected individuals in the population at time t by adjusting ρˆ+V(t) to

(2)ρˆ+ V*t=∑j=1ρˆ+∑i=1nIXij−V0ij, Xij−BijV0ij−V1ijt

where Bi(j)∼iidBernoulli(1−pa), and p_a is the probability that an infection is asymptomatic. This expression ρˆ+V*(t) is the statistic we base inference upon, and its expected value is the focus of our inference.

Our computation of an estimate of E[ρˆ+ V*(t)] proceeds using statistical simulation, using both the bootstrap and Monte Carlo integration. (We previously described the Monte Carlo integration as akin to imputation, as this approach samples from or imputes values from specified densities in computing integrals against these distributions.) The motivation for this approach is seen in the following. For each i=1,2,…,n, assume that Xi∼FiX, for cumulative distribution functions (CDF) FiX. Write F_i for the CDF of Y_i. Then, assuming ρˆ is a consistent estimator of ρ and is independent of all other random variables (which will be the case when estimated from an exogeneous source), with E[ρˆ+]=ρ+,

(3)E[ρˆ+ V*(t)]=∑i=1nE[ρˆ+]P[Xi−V0i<t<Xi−Bi(V0i−V1i)]=∑i=1nρ+ P[t+Bi(V0i−V1i)<Xi<t+V0i]=∑i=1nρ+(P[Xi<t+V0i]−P[Xi<t+Bi(V0i−V1i)])=∑i=1nρ+{∫0∞FiX(t+v0)fV0(v0)dv0−∫0∞∫0∞∑b=01FiX(t+b[v0−v1])pab(1−pa)1−bfV0(v0)fV1(v1)dv0 dv1}=ρ+{∫0∞GX(t+v0)fV0(v0)dv0−∫0∞∫0∞∑b=01GX(t+b[v0−v1])pab(1−pa)1−bfV0(v0)fV1(v1)dv0 dv1}=ρ+{∫0∞GX(t+v0)fV0(v0)dv0−pa∫0∞∫0∞GX(t+v0−v1)fV0(v0)fV1(v1)dv0 dv1−(1−pa)GX(t)}

where GX(u)=∑i=1nFiX(u) and noting that ρ₊ is computable from the distribution of ρˆ.

Computation of the integrals in last expression may proceed directly, either analytically or numerically, if we assume particular distributions for the FiX,fV0, and fV1, and values for p_a and ρ. A sensible, computationally accessible alternative is to couple the bootstrap with Monte Carlo integration. We enumerate the computational steps involved after listing our distributional assumptions. The estimates for the F_i, and so for G, follow from the use of smoothed bootstrap resampling of the original onset time observations Y_i and subsequent smoothing to X_i; a specific motivation for this is detailed in Development using a Poisson process model for the underlying case onset times below. In cases when estimates for the mean incubation period are unavailable, we assume a quartic density for V₀, so that

fV0(v0)=158(d−c)(1−[2v0−(c+d)d−c]2)2I[c,d](v0)

where d>c>0 are chosen to correspond to the distribution for duration from infection to symptom onset. We adopted the quartic density when the only information we had on this variable was that the incubation period is roughly c to d days. The density f_V0 is symmetric with mean (c+d)/2 and variance (d−c)/20. Sampling from the quartic may be done using, for example, acceptance-rejection sampling. When more information regarding the distribution of V0 is available, for example if there are published estimates of the mean and standard deviation of onset duration, we use a gamma density for fV0(v0).

In the original application for WNV transfusion risk via blood transfusion (Biggerstaff and Petersen 2002), when considering how to utilize available information for the distribution of V₁ to inform specification of fV1(v1), the density representing the duration of WNV infection, we resampled from the observed values directly reported in Southam and Moore (1954) as representative of the distribution (so we employed bootstrapping). To be able to expand our methodology, we subsequently adopted the Weibull density for this duration, and we estimated the shape and scale parameters by matching moments. We use the Weibull below, noting little difference in our experience when the parameters for the Weilbull are specified to match the sample mean and sample variance of the data in Southam and Moore (1954).

Because we use reported estimates pˆa and ρˆ, we incorporate uncertainty in these estimates directly into the Monte Carlo integration as well. In particular, we assume a beta density for pˆa and a gamma density for ρˆ.

Estimation algorithm

With the distributional assumptions and approximations given, we compute the estimate of E[ρˆ+ V*(t)] as follows. For each of j=1,2,…,J simulations:

Resample with replacement n values from the originally observed case onset dates Y1,Y2,…,Yn: write these Yk(j), for k∈1,2,…,n; throughout this algorithm, k indexes the n resampled cases.
Compute the smoothing parameter h_j using the n values Yk(j).
Generate n independent standard normal deviates, εk(j), and then compute Xk(j)=Yk(j)+hj εk(j).
Generate n independent values V0k(j) from an assumed distribution fV0(v0).

In Biggerstaff and Petersen (2002) we used a quartic density with a=1 and b=5 and in subsequent applications used a gamma density with mean incubation period, μˆ0, and variance, σ02, with these moments the same as for the stated quartic based on references given in (Biggerstaff and Petersen 2002). To be conservative, we began to use the a gamma density to permit longer incubation periods than permitted by the quartic, which is bounded.

Generate n independent values V1k(j) from a specified distribution fV1(v1).

We use a Weibull density with mean μˆ1 and variance σˆ12, where these are specified from available estimates.

Generate a value ρˆ (j) from an assumed distribution fρˆ(r); round this value up to the next integer and write as ρˆ+ (j).

We use a Weibull density with mean μˆρ and variance σˆρ2, where these are specified from available estimates.

Replicate each the resampled case infection time, (Xk(j)−V0k(j),Xk(j)−V0k(j)+V1k(j)),ρˆ+ (j) times; identify these replicates using index r with superscripts (j,r) rather than (j).
Generate a value pˆa(j) from an assumed distribution fpˆ(p).

We use a beta distribution with mean μˆpa and variance σˆpa2, where these are specified from available estimates.

Generate n×ρˆ+ independent values Bk(j,r) from the Bernoulli distribution with probability 1−pˆa(j). This will allow independent truncation of 100(1−pˆa(j))% of the total number of generated “population” infections’ WNV infection times in simulation j.

With these values in-hand for each simulation j=1,2,…,J, compute on a grid of values 0=t1<t2<⋯<tg<⋯<tG=T, for g=1,2,…,G,

Eρ+ V*tgˆ=1J∑j=1J∑k=1n∑r=1ρˆ+IXkj,r−V0kj,r, Xkj,r−Bkj,rV0kj,r−V1kj,rtg=1J∑j=1JWjtg

This curve is an average of J simulations of step functions, Wj(tg), where the heights of the steps are the estimated number of asymptomatic, viremic indivdiuals in the population at each time point, where we “connect” or linearly interpolate values between grid points. For simplicity, write E(t)=E[ρ+ V*(t)]ˆ.

Confidence bands

We compute confidence bands based on the simultaneous bootstrap-t approach used in intensity estimation for Poisson processes, as derived in Cowling, Hall, and Phillips (1996). We use the Poisson approximating standard deviation and compute

Tj(t)=Wj(t)−E(t)Wj(t).

for j=1,2,…,J. To compute 100(1−α)% simultaneous confidence bands C1−α over t∈(0,T), find q₁ and q₂ such that

P[q1≤Wj(t)≤q2,∀t∈(0,T)|X]=α

and

PWjt≤q1,∀t∈0,T|X=PWjt≥q2,∀t∈0,T|X,

for equal-tailed bands, where X represents all the data and exogeneous information available. Given these constraints, set

C1−α={(t,y):t∈(0,T),max[0,E(t)−q2 σˆ]<y<E(t)−q1 σˆ},

where σˆ=E(t).

Finally, the sought-after, time-dependent risk of a WNV-infected donation is computed directly as the estimated number of such infections, E(t), divided by the population size, N, and we write P(t)=E(t)/N. Similarly, the confidence bands for the risk are computed as C1−α/N.

Summary measures

As part of the simulation estimation, we also estimate the mean risk over the course of the outbreak by computing the average height of the risk curve. We do this by computing, for each simulation j, the area under the curve Wj(t) using Simpson’s Rule, then dividing this by the duration of the outbreak, T. The mean of the simulated areas gives the final estimate of the average risk, and we use the 2.5th and 97.5th quantiles of these simulated areas as 95% confidence intervals for the average risk. Lastly, we estimate the maximum risk over the course of the outbreak by reading off the maximum height of the estimated curve, along with a 95% confidence interval for this maximum from the confidence bands for the curve; the location of the maximum gives an estimate of the date of maximum risk.

Refinements for estimation of WNV in cadaveric tissue donations

In certain applications, as with the motivating example here of estimation of WNV infection in cadaveric tissues, stratification by variables for which key parameters vary should improve estimation. In the present example, the inapparent-to-apparent infection ratio is known to vary by age, with symptomatic infection (and severity of illness) more likely with increasing age. A straightforward refinement of the methods in Methodological and computational development to incorporate this is to use stratification as follows. Assume that the population is stratified into S strata of sizes N_s, for s=1,2,…,S, and let N=∑sNs be the population size. Indicate case onset times, etc., in stratum s by adding this index to the subscript throughout, so write, for example, Y_si, X_si, Vs(t),Vs*(t),ρs,Ws j(t),Es(t),Ps(t). Stratum-specific estimates Ps(t)=Es(t)/Ns are readily computed by restricting the methods of Methodological and computational development to strata of interest.

Combining the strata for a population-level estimate of the risk may proceed naturally by weighting the stratum-specific estimates by relative population size, Ns/N. This weighting provides an overall estimate of the number of asymptomatic and infected donors. In some applications, however, further weighting to refine these estimates to individuals who either are more likely to donate, if this information is available, or who are actually eligible to donate may be warranted. This is the case in the cadaveric tissue donation application, where the proportions of asymptomatic and infected individuals are not expected to be the same as or even equally proportional among age groups, since donors must be recently deceased. Because death rates differ by age and there are available age group-specific inapparent-to-apparent infection ratio estimates, we are able to include weighting by stratum-specific death rates to account for the different probabilities of asymptomatic and infected individuals becoming deceased (and so eligible tissue donors). For s=1,2,…,S, let δ_s be the stratum-specific (age group-specific) annual death rates, and let Δs=T365δs be the death rates adjusted to the observed outbreak time duration of interest (where we have assumed that T is in days, and note that T365 will be greater than 1 for multi-year estimates, and we can naturally adjust for leap years). For each stratum, then, ΔsEs(t) is an estimate of the number of asymptomatic, WNV-infected potential cadaveric donors at time t, where we do note that the the use of death rates this way is intended to reflect the cumulative death rates over the whole of the observed outbreak (0,T), apportioned uniformly over the time of interest; time-dependent deaths rates [δs(t)] could be incorporated directly if available. The corresponding stratum estimates for the risk of an asymptomatic, viremic donor are therefore ΔsEs(t)/Ns, so that these are combined to the population level, incorporating relative population size, as

(4)P(t)=∑s=1SNsNΔs Es(t)Ns=1N∑s=1SΔs Es(t)=1NJ∑s=1S∑j=1JΔsWsj(t)

where we use for convenience the same number of simulations J for each stratum. From the expression (4), we see that the order of computation—by stratum or by simulation—does not matter, and our implementation is to compute by simulation. Doing this facilitates computation of bootstrap confidence bands using P(t), which are computed as in Methodological and computational development based on the resampled estimates Wj(t)=∑sΔs Wsj(t).

Development using a Poisson process model for the underlying case onset times

In the Monte Carlo evaluation of the integrals in the expression for E[ρˆ+ V*(t)] in Methodological and computational development, we used Gaussian smoothing of the observed Y_i and bootstrap resampling (so, effectively, the smoothed bootstrap (Davison and Hinkley 1997, § 3.4)) as essentially the estimation of the cumulative distrbution functions FiX, and so of GX(u)=∑iFiX(u). On its face, this does not assume a particular distributional form for the stochastic process {Yi}.

We noted without detail in the Discussion section of Biggerstaff and Petersen (2002), however, that this approach can be seen as assuming that the process {Yi} follows an inhomogeneous Poisson process in which a nonparametric, kernel smoothing estimate for the cumulative rate function of the process is used. For completeness, we now detail this observation here.

Begin by assuming that the WNND case onset dates {Yi} (ignoring the rounding of the observations to the day) are observations from an inhomogeneous Poisson process with rate function λ(t) and cumulative rate function Λ(t)=∫0tλ(u)du. For this model, the CDF Fi(t) of Y_i is (Liao 2014)

Fi(t)=1−e−Λ(t)∑j=0i−1[Λ(t)]jj!=1−Π(i−1;Λ(t))

introducing the notation Π(x;Λ(t)) for the CDF of a Poisson random variable with rate Λ(t). From this it follows that

(5)G(t)=∑i=1nFi(t)=n[1−Π(n−1;Λ(t))+Λ(t)nΠ(n−2;Λ(t))]

Direct estimation of G(t), therefore, reduces to estimation of Λ(t). Various options are available for this, of course, including parametric regression, nonparametric spline regression, and nonparametric kernel smoothing. Employing the last of these and using a Gaussian kernel gives

Λˆ(t)=1n∑i=1n1uϕ(t−Yiu)

for smoothing parameter u>0 and standard normal density function ϕ(⋅). Use of Λˆ(t) in Eq. (5), then, provides an estimate of G(t), which is GXt above after the estimation of Λ(t), as included in Eq. (3). In doing this, we then may employ bootstrap resampling of the Yi for estimation using Λˆ(t) in the Monte Carlo computations described above, so that Λˆ(t) is computed for each such bootstrap sample. Such sampling is, however, the same as sampling from the kernel smoothing estimate Λˆ(t) and then computing the integral in Eq. (3). Further, sampling from a Gaussian-smoothed kernel estimate of the Y_i is equivalent to sampling from the original data (Y_i) and adding scaled Gaussian noise, with scale equal to the smoothing parameter in the kernel smoothing estimate; this is exactly the smoothed bootstrap (e.g., Davison and Hinkley 1997, §3.4).

Operationally, then, using bootstrap resampling and then smoothing the resampled values using Gaussian noise and then using this in the Monte Carlo integrations in Eq. (3) is equivalent to assuming an inhomogeneous Poisson Process for the {Yi} and then estimating the cumulative rate function using a Gaussian kernel smoother, then evaluating the integral in Eq. (3).

Computer implementation

The methods presented were implemented in mostly purpose-written software using the statistical package R (R Core Team 2018) https://www.R-project.org/Also used in these analyses were the boot package (Canty 2017), and the chron package (James 2018) Software is available from the author upon request. To implement the methods in other settings, data with disease case onset times and summary estimates of the modeling parameters—incubation period, inapparent-to-apparent infection ratio, viremia duration, asymptomatic rate–are needed as inputs. Stratum-specific information for the case data and modeling parameters is required to implement the weighting methods described.

Application: cardiac tissue, Colorado, 2003

We implement the updated methods with the application introduced above, estimation of the risk of WNV infected cardiac tissue during the Colorado WNV outbreak in 2003, which remains at the high end of state-level WNV outbreak experience (http://www.cdc.gov/westnile/statsmaps/). Figure 1 shows the number of WNND cases aged 13–55 years by date of symptom onset in Colorado in 2003. In this case, a total of 358 WNND cases were reported to CDC in this age range, for an outbreak lasting 99 days. Cardiac tissue donations are restricted to donors under 56 years of age, and because there were few WNND cases aged less than 13 years, these were analyzed separately in Biggerstaff and Petersen (2020). The data in Figure 1 therefore provide time-dependent information on the WNV transmission activity in Colorado throughout 2003 for potential cardiac tissue donors.

As detailed in Biggerstaff and Petersen (2020), estimates of the inapparent-to-apparent infection ratios, ρˆs, and their standard errors are available by age group stratum (Carson et al. 2012), as reported in Table 1. Note that these age grouping were selected to reflect available data and risk categories used for cardiac tissue by tissue banking industry and regulatory agencies. Also shown are the total numbers of WNND cases reported in these age groups, the annual death rates for these age groups for Colorado (http://wonder.cdc.gov), and the total population for these age groups.

Table 1:

Age group-specific WNND counts for Colorado 2003; inapparent-to-apparent infection ratios, ρˆ, and their standard errors (SE) (Biggerstaff and Petersen 2020); annual death rates, δ, per 1,000 population; and total population (Pop).

Age group	WNND	ρˆ (SE)	δ
13–20	22	742.57 (96.92)	0.58
21–40	134	526.90 (58.39)	1.12
41–55	202	286.23 (19.52)	3.86
Total	358		Pop = 2,911,979

The gamma density was assumed for the incubation period, with mean μˆ0=3 (SD σˆ0=4/7) days (Biggerstaff and Petersen 2002). The mean duration of WNV infection in cardiac tissue was estimated as μˆ1=9.39 (SD σˆ1=16.68) days based on a Weibull survival model fit to data recorded in Southam and Moore (1954) (see Biggerstaff and Petersen 2020). Additionally, assuming the asymptomatic rate of μˆpa=0.79 (SD σˆpa=0.13) as estimated in Mostashari et al. (2001), we apply the simulation algorithm described above and employ the refinements detailed in Refinements for estimation of WNV in cadaveric tissue donations using J=1000 simulations to produce the estimated risk of WNV-infected cardiac tissue donations shown in Figure 2. The solid black line is the risk curve, while the black, dashed lines are the 95% confidence bands. To illustrate the variability in the simulations, 100 of the simulated curves Wj(t) are shown in gray.

$Figure 2: Estimated risk of WNV-infected cardiac tissue among potential donations, Colorado, 2003, shown as the black, solid line, with 95% confidence bands shown as black, dashed lines. The gray lines are a sample of the simulated risk curves,Wj(t)${W}_{j}\left(t\right)$.$

Figure 2:

Estimated risk of WNV-infected cardiac tissue among potential donations, Colorado, 2003, shown as the black, solid line, with 95% confidence bands shown as black, dashed lines. The gray lines are a sample of the simulated risk curves,Wj(t).

Finally, the average risk of WNV-infected cardiac tissue donations over the duration of this outbreak was 6.3 (95% CI 3.1, 9.6) per 10,000 deaths, with a maximum risk of 12.3 (95% CI 6.8, 20.5) per 10,000 deaths on August 10.

Discussion

While the original methods described here were developed for WNV, the underlying ideas were general. Subsequent applications to severe acute respiratory syndrome (SARS) virus (Shang et al. 2007), chikungunya virus (Brouard et al. 2008; Liumbruno et al. 2008), dengue virus (Seed et al. 2009; Petersen, Tomashek, and Biggerstaff 2012), hepatitis A virus (Perevoscikovs et al. 2010), and Ross River virus (Shang et al. 2016) show that the methods may be applied in situations for which sufficient information is available on viral/pathogen kinetics, symptomatic rates, inapparent-to-apparent infection ratios, and, in some cases, surveillance coverage. In at least one application (Seed et al. 2009), the methods were applied in real time to guide decisions on suspension and continuance of blood donation collection during an outbreak. Most of these applications have been retrospective, so that case surveillance has been completed before analysis. In real-time application, time for disease progression and reporting delays of identified cases to public health surveillance systems can be expected to cause under-representation of real-time disease case counts. Using observed counts without correction for this would result in under-estimation, since too few cases would be identified in real-time because of these delays. Statistical methods of prediction of occured but not yet reported events may be considered in such cases to attempt to adjust for this expected downward bias (see, e.g., Lawless 1994).

Application of these methods yields estimates of risks of infected donations before any post-donation testing or processing. In the US, for example, blood donations are now routinely tested seasonally for WNV, using policies implemented shortly after risk estimates were published and WNV transfusion transmission was documented (U. S. Food and Drug Administration 2009). Further, donated tissues undergo post-donation processing to reduce or eliminate transmission of multiple pathogens. Risk estimates produced using the methods developed here are intended to inform applicable industries, public health officials, and regulatory agencies concerning these risks to guide requirements for intervention strategies and policies. If estimates of the impact of post-processing are available, they may be naturally applied to results of the present analyses to characterize post-processing risks.

One assumption we made was that individuals share common distributions for both the incubation period and the viral infection duration. This is likely not strictly the case in either the donor or the general population, as there may be biological factors that mediate these for different individuals or groups of individuals. For WNV there is not currently evidence that systematic differences by known factors for these quantities would be on a scale to impact the general conclusions drawn from the methods used here, which take account of known uncertainties directly in the estimation. Should such become available, the estimation framework presented allows direct incorporation of such information in an easily implemented way. We also assumed that durations of onset and of infection were the same for individuals who remain asymptomatic or who develop symptoms, as there are not available separate estimates for these durations between these two groups; should such estimates become available, they can be incorporated into the methods presented here easily using the stratification mechanism introduced. Further, if any of the parameters differ by identifiable strata (say, by sex or by age), then one may naturally incorporate any available estimates or information into the framework we introduce here for weighted, stratified analyses, as well. In the absence of such detailed information, however, we may view these assumptions as reflecting a population average distribution across any such relevant factors.

In our approach we relied on Gausian kernel smoothing to underlie the estimation of time-dependent risk curves. Naturally, other methods of smoothing, such as spline smoothing, semi- or fully parametric regression, or wavelet basis smoothing methods may have been used. We adopted the approach we present due to the the direct and ready interpretability of the conversion of the initial, discrete observations to continuous “observations” (Y_i to X_i) and to utilize the bootstrap methods in (Cowling, Hall, and Phillips 1996). Evaluation of other methods of smoothing in the risk estimation context considered here may prove an interesting future research topic.

The analysis estimating risk of WNV infection in donated cardiac tissue during the Colorado outbreak of 2003 was part of a more complete analysis reported in Biggerstaff and Petersen (2020). Interested readers may refer there for estimates for blood, and cardiac, musculoskeletal, and soft tissues for outbreaks covering the contiguous US, 1999–2014; Colorado 2003; and Dallas-Fort Worth 2012.

We have detailed and updated our approach for time-dependent risk estimation for arbovirus infection of blood donations, and to address unique aspects of risk estimation for cadaveric tissue donations, we expanded the methods to include stratification and weighting. We expect that these refinements will provide greater flexibility and increased precision for researchers employing this general approach for a general variety of pathogens.

Corresponding author: Brad J. Biggerstaff, Centers for Disease Control and Prevention, Division of Vector-Borne Diseases, Fort Collins, Colorado, USA, E-mail: bkb5@cdc.gov

Acknowledgment

The author thanks the American Association of Tissue Banks for the suggestion to pursue this work, particularly Mr. Frank Wilton, former AATB President and CEO, and also the Scientific and Technical Affairs Committee of the American Association of Tissue Banks for technical guidance, particularly the former committee chair, Dr. Jeffrey Cartmell. Further, the author thanks Dr. Lyle R. Petersen, Division Director, Division of Vector-Borne Diseases, CDC, for his collaboration on these projects and for encouragement in writing this methods paper. Thanks also to the ArboNET team at CDC for curating the data used in this study and to the SciComp team at CDC for access to their computational resources. Finally, thank you to Dr. Mark Delorey of DVBD/CDC for his careful review of the manuscript.

The findings and conclusions herein are those of the authors and do not necessarily represent the official position of the US Centers for Disease Control and Prevention.

Research funding: None declared.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Competing interests: Authors state no conflict of interest.

References

Biggerstaff, B. J., and L. R. Petersen. 2002. “Estimated Risk of West Nile Virus Transmission through Blood Transfusion during an Epidemic in Queens, New York City.” Transfusion 42: 1019–26, https://doi.org/10.1046/j.1537-2995.2002.00167.x.Search in Google Scholar

Biggerstaff, B. J., and L. R. Petersen. 2003. “Estimated Risk of Transmission of the West Nile Virus through Blood Transfusion in the US, 2002.” Transfusion 43: 1007–17, https://doi.org/10.1046/j.1537-2995.2003.00480.x.Search in Google Scholar

Biggerstaff, B. J., and L. R. Petersen. 2020. “Estimated Risk of West Nile Virus Infection in Cadaveric Tissue Donations.” Unpublished.Search in Google Scholar

Brouard, C., P. Bernillon, I. Quatresous, J. Pillonel, A. Assal, H. De Valk, and J. C. Desenclos, A Workgroup, Quantitative Estimation of the Risk of Blood Donation Contamination by Infectious. 2008. “Estimated Risk of Chikungunya Viremic Blood Donation during an Epidemic on Reunion Island in the Indian Ocean, 2005 to 2007.” Transfusion 48 (7): 1333–41, doi:https://doi.org/10.1111/j.1537-2995.2008.01646.x.Search in Google Scholar

Canty, A. and B. Ripley. 2017. boot: Bootstrap R (S-Plus) Functions. R package, version 1.3-20.Search in Google Scholar

Carson, P. J., S. M. Borchardt, B. Custer, H. E. Prince, J. Dunn-Williams, V. Winkelman, L. Tobler, B. J. Biggerstaff, R. Lanciotti, L. R. Petersen, and M. P. Busch. 2012. “Neuroinvasive Disease and West Nile Virus Infection, North Dakota, usa, 1999-2008.” Emerging Infectious Diseases 18: 684–6, https://doi.org/10.3201/eid1804.111313.Search in Google Scholar

Cowling, A., P. Hall, and M. J. Phillips. 1996. “Bootstrap Confidence Regions for the Intensity of a Poisson Point Process.” Journal of the American Statistical Association 91: 1516–24, https://doi.org/10.1080/01621459.1996.10476719.Search in Google Scholar

Davison, A., and D. Hinkley. 1997. Bootstrap Methods and Their Application, Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge: Cambridge University Press.10.1017/CBO9780511802843Search in Google Scholar

U. S. Food and Drug Administration. 2009. Guidance for Industry: Use of Nucleic Acid Tests to Reduce the Risk of Transmission of West Nile Virus from Donors of Whole Blood and Blood Components Intended for Transfusion. Also available at https://www.fda.gov/media/124300/download.Search in Google Scholar

Centers for Disease Control and Prevention. 2002. “West Nile Virus Infection in Organ Donor and Transplant Recipients–Georgia and Florida, 2002.” MMWR Morb Mortal Wkly Rep 51 (35): 790. https://www.cdc.gov/mmwr/preview/mmwrhtml/mm5135a5.htm.Search in Google Scholar

Iwamoto, M., D. B. Jernigan, A. Guasch, M. J. Trepka, C. G. Blackmore, W. C. Hellinger, S. M. Pham, S. Zaki, R. S. Lanciotti, S. E. Lance-Parker, C. A. DiazGranados, A. G. Winquist, C. A. Perlino, S. Wiersma, K. L. Hillyer, J. L. Goodman, A. A. Marfin, M. E. Chamberland, and L. R. Petersen and The West Nile Virus in Transplant Recipients Investigation Team. 2003. “Transmission of West Nile Virus from an Organ Donor to Four Transplant Recipients.” New England Journal of Medicine 348: 2196–203, https://doi.org/10.1056/nejmoa022987.Search in Google Scholar

James, D., and K. Hornik. 2018. chron: Chronological Objects which can Handle Dates and Times. R package, version 2.3-53.Search in Google Scholar

Lawless, J. 1994. “Adjustments for Reporting Delays and the Prediction of Occurred but Not Reported Events.” Canadian Journal of Statistics 22: 15–31, https://doi.org/10.2307/3315826.n1.Search in Google Scholar

Liao, M. 2014. Applied Stochastic Processes. Boca Raton, Florida: CRC Press.Search in Google Scholar

Liumbruno, G. M., D. Calteri, K. Petropulacos, A. Mattivi, C. Po, P. Macini, I. Tomasini, P. Zucchelli, A. R. Silvestri, V. Sambri, S. Pupella, L. Catalano, V. Piccinini, G. Calizzani, and G. Grazzini. 2008. “The Chikungunya Epidemic in Italy and its Repercussion on the Blood System.” Blood Transfus 6: 199–210, https://doi.org/10.2450/2008.0016-08.Search in Google Scholar

Mostashari, F., M. L. Bunning, P. T. Kitsutani, D. A. Singer, D. Nash, M. J. Cooper, N. Katz, K. A. Liljebjelke, B. J. Biggerstaff, A. D. Fine, M. C. Layton, S. M. Mullin, A. J. Johnson, D. A. Martin, E. B. Hayes, and G. L. Campbell. 2001. “Epidemic West Nile Encephalitis, New York, 1999: Results of a Household-Based Seroepidemiological Survey.” Lancet 358: 261–4, https://doi.org/10.1016/s0140-6736(01)05480-0.Search in Google Scholar

Perevoscikovs, J., A. Lenglet, I. Lucenko, A. Steinerte, L. Payne Hallstrom, and D. Coulombier. 2010. “Assessing the Risk of a Community Outbreak of Hepatitis a on Blood Safety in Latvia, 2008.” Euro Surveillance 15: 19640, https://doi.org/10.2807/ese.15.33.19640-en.Search in Google Scholar

Petersen, L. R., K. M. Tomashek, and B. J. Biggerstaff. 2012. “Estimated Prevalence of Dengue Viremia in Puerto Rican Blood Donations, 1995 through 2010.” Transfusion 52: 1647–51, https://doi.org/10.1111/j.1537-2995.2011.03529.x.Search in Google Scholar

R Core Team. 2018. R: A language and environment for statistical computing. R Foundation for Statistical Computing. Vienna, Austria. URL https://www.R-project.org/.Search in Google Scholar

Scott, D. W. 1992. Multivariate Density Estimation. Theory, Practice and Visualization. New York: Wiley.10.1002/9780470316849Search in Google Scholar

Seed, C. R., P. Kiely, C. A. Hyland, and A. J. Keller. 2009. “The Risk of Dengue Transmission by Blood during a 2004 Outbreak in Cairns, Australia.” Transfusion 49: 1482–7, https://doi.org/10.1111/j.1537-2995.2009.02159.x.Search in Google Scholar

Shang, G., B. J. Biggerstaff, A. M. Richardson, M. E. Gahan, and B. A. Lidbury. 2016. “A Simulation Model to Estimate the Risk of Transfusion-Transmitted Arboviral Infection.” Transfusion and Apheresis Science 55: 233–9, https://doi.org/10.1016/j.transci.2016.07.018.Search in Google Scholar

Shang, G., B. J. Biggerstaff, B. Yang, C. Shao, and A. Farrugia. 2007. “Theoretically Estimated Risk of Severe Acute Respiratory Syndrome Transmission through Blood Transfusion during an Epidemic in Shenzhen, Guangdong, china in 2003.” Transfusion and Apheresis Science 37: 233–40, https://doi.org/10.1016/j.transci.2007.09.004.Search in Google Scholar

Sheather, S., and C. Jones. 1991. “A Reliable Data-Based Bandwidth Selection Method for Kernel Density Estimation.” Journal of the Royal Statistical Society: Series B 53: 683–90, https://doi.org/10.1111/j.2517-6161.1991.tb01857.x.Search in Google Scholar

Southam, C. M., and A. E. Moore. 1954. “Induced Virus Infections in Man by the Egypt Isolates of West Nile Virus.” The American Journal of Tropical Medicine and Hygiene 3: 19–50, https://doi.org/10.4269/ajtmh.1954.3.19.Search in Google Scholar

Williamson, P. C., B. Custer, B. J. Biggerstaff, R. S. Lanciotti, M. H. Sayers, S. J. Eason, M. R. Dixon, V. Winkelman, M. C. Lanteri, L. R. Petersen, and M. P. Busch. 2016. “Incidence of West Nile Virus Infection in the Dallas-Fort Worth Metropolitan Area during the 2012 Epidemic.” Epidemiology and Infection 145: 2536–44, https://doi.org/10.1017/s0950268816000042.Search in Google Scholar

Received: 2020-01-28

Accepted: 2020-10-12

Published Online: 2020-11-13

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/scid-2020-0001

Keywords for this article

blood safety; Poisson process; public health; tissue safety

Creative Commons

BY 4.0