Here we provide a brief description of BARD. A thorough discussion of this system can be found in Hogan et al. [3]. To simplify our exposition, we follow the terminology in Lawson and Kleinman [1] and generically refer to the basic administrative spatial units employed for aggregating the disease data (e.g., block groups or zip codes) as tracts. The main symbols used in this paper and their respective meanings are listed in Table 1.

Our disease data--counts of ED visits with RCC--correspond to a region surrounding the city of Pittsburgh. This region consists of seven counties: Allegheny, Armstrong, Beaver, Butler, Lawrence, Washington, and Westmoreland. The historical ED data for our experiments were provided by 10 EDs operated by one health system. These data represented nearly 30% of all ED visits in the surveillance region. We divided the total period spanned by the available ED data into a training period (1 January 1999 to 31 December 2004), which was used to train BARD's detection algorithm, and a test period (1 January 2005 to 31 December 2005), used in our evaluation experiments. Finally, the weather data were provided from the National Weather Service, while the populations and central zip codes came from the ESRI ArcGIS Desktop product.

The data that describe the commuting patterns were provided by the U.S. Census Bureau. These data were collected during the 2000 Census, have national coverage and are provided at the census tract level: each commuting flow denotes the average daily number of commuters between a home census tract and a work census tract. The commuting flows can be naturally modelled by a weighted, directed graph in which nodes denote tracts, arcs represent commuting flows, and the weight of an arc denotes the number of commuters in the corresponding flow. We extracted from the nationwide commuting dataset the commuting flows for which both the residence tract and the work tract belong to our surveillance region. The total number of commuters in this intra-region subset of flows was 1,005,566.

As a pre-processing step, we needed to transform the commuting flows from the census tract level provided by the U. S. Census Bureau to one of the two levels supported by the BARD simulator: zip code, or block group. In this paper we performed the latter transformation mainly because the block-group level is finer-grained than the zip-code level and thus has the potential to lead to a more accurate estimation of the release location by BARD. Note that block group is a sub-unit of census tract (i.e., each census tract is partitioned into two or more block groups). To perform the desired transformation, we split each commuting flow between a pair of census tracts T₁, T₂ into several smaller-sized flows according to the following two intuitive rules: (i) the number of workers coming from each constituent block group of T₁ was assumed to be proportional to the block group population; (ii) the number of workers going to T₂ was assumed to be evenly divided among its constituent block groups. The final commuting graph for our region consisted of 1991 nodes (block groups) and 324,402 arcs (commuting flows).

We refer to the version of BARD simulator that takes commuting into account as the BARD-C simulator. BARD-C simulator takes as input a representation of the commuting graph G, in addition to the surveillance region and the anthrax release parameters. The key observation employed in the incorporation of commuting in simulation is that if a release occurs during working hours, a worker would be exposed at the dosage level corresponding to his/her work tract rather than the dosage level of his/her home tract. This idea is implemented in BARD-C by modifying the Step 3 of Algorithm 1, namely the computation of the ED-visit count for each tract. To realize this modification BARD-C first computes the ED-visit count corresponding to each commuting flow separately and then re-aggregates the workflow-level counts according to the home tract. More precisely, the number of cases for each tract i is computed as

(5)

where

(6)

In eqs. (5)-(6), c_i denotes the total number of cases who live in tract i, Out(i) ={i} ∪ {j|(i, j) is an arc of G} denotes the union of set {i} with the set of out-neighbours of tract i in the commuting graph G (i.e., the tracts where people living in tract i go to work), n_ij, i ≠ j, denotes the weight of the arc (i, j) in the graph G, i.e., the number of people living in tract i and working in tract j, and n_ii denotes number of people who either don't work, or who both live and work in tract i. Finally, c_ij, i ≠ j, denotes the number of cases who live in tract i and work in tract j, c_ii denotes number of cases who either don't work, or who both live and work in tract i and θ_j denotes the probability that a person who is exposed in tract j will visit an ED, given that a release has occurred.

This refined computation of the tract-level ED counts in BARD-C could alternatively be implemented by modifying the ED-visit probability of each individual who lives in tract i (i.e., modifying the Step 2 of Algorithm 1) prior to generating the number of cases for tract i as a binomial random variate. Indeed, let us denote by the probability that an individual who lives in tract i will visit an ED, given that a release has occurred and assuming that commuting has taken place. Let the parameter θ_i denote, as before, the probability that an individual who is exposed to anthrax in tract i will visit an ED. By approximating the binomial distribution with a Poisson distribution and recalling that the sum of independent Poisson random variables is also a Poisson random variable, it is straightforward to show that the method for generating the count c_i specified by eqs. (5)-(6) is essentially the same as generating the count c_i from a binomial distribution Bin(n_i, ), where

(7)

The above two methods for incorporating commuting in simulation differ slightly from the method employed by Buckeridge [4]. In [4], other non-commuting types of mobility were also taken into account and therefore the exposure tract could be different from both the residential and the work tract. Hence, in [4] the incorporation of mobility was based on the computation of conditional probabilities P(exposed in tract i|t, live in tract j) for all pairs i, j.

After generating the number of cases for each tract, BARD-C computes the ED-visit time of each simulated case by, again, taking into account the fact that a worker would be exposed at the dosage level of his/her work tract. This idea is implemented in BARD-C by ensuring that in the Step 4(a) of Algorithm 1 the time of ED visit for a worker case is computed as a random deviate from the lognormal distribution that corresponds to the dosage level at his/her work tract.

The incorporation of commuting in simulation leads to an increase in the running time of the Step 3 of BARD's simulation algorithm by a constant factor. It is straightforward to see that this factor equals the average out-degree of the commuting graph. For the block group-level commuting graph of the Pittsburgh region that we created, the average out-degree is nearly 150. However, the factor by which the total running time of BARD simulator is increased due to the incorporation of commuting is smaller. The reason for this is that Step 3 accounts for only a fraction of the running time of BARD. In a computer with a 3 GHz processor and 2 GB of main memory, a simulation for the Pittsburgh region took nearly 2 minutes with BARD and nearly 20 minutes with BARD-C.

To verify that BARD-C simulator works as intended, we generated a random release parameter vector, consisting of the release time, location, and quantity. Then, using this fixed release parameter vector we produced 10 synthetic outbreaks with BARD-C simulator. Next, for each tract i of the surveillance region, we fit a binomial model Bin(n_i, ) to the sample of simulated counts . As before, n_i denotes the known population of tract i, and denotes the probability that an individual who lives in tract i will visit an ED, given that a release has occurred and assuming that commuting has taken place. We denote by the estimate of the binomial probability parameter learned from the count data generated with BARD-C simulator. Finally, we computed an analytical prediction for each parameter . We achieved this by first using the BARD simulator to compute the binomial probability parameters θ_j for all j ∈ Out(i) and then computing the analytical prediction of the parameter through eq. (7).

We have also developed a method for integrating commuting with BARD detector. A thorough discussion of this method is given [10]. We refer to the version of BARD detector that takes commuting into account as the BARD-C detector.

We conducted an experiment to measure whether BARD's detection and characterization performance is significantly affected by the incorporation of commuting in outbreak simulation. In the remainder of the paper, we refer to this experiment as the "simulation study" because the emphasis of this study was on the integration of commuting with BARD simulator. In the next section we discuss a second study, where the emphasis is on the integration of commuting with BARD detector. We refer to this second study as the "detection study". Figure 1 illustrates the difference between the objectives of the simulation and detection studies.

The simulation study followed a matched-pairs design. First, a set of 185 release parameter vectors r = (t, x, y, h, q) was generated. For all parameter vectors, the quantity q was fixed at 0.5 kg. The release time t was chosen uniformly at random from the test period, i.e., the year 2005. The parameters x and y where chosen uniformly at random from a circular sub-region of the surveillance region centred in the city of Pittsburgh and having a radius of 40 km; this sub-region was found to have better ED-data coverage than the whole region. Finally, the height h was chosen from a non-uniform distribution that favours lower release heights relative to the higher ones.

We supplied each release vector r_i, i = 1,...,185, as input to the BARD and BARD-C simulators. By doing so, we obtained one group of 185 simulations that did not take commuting into account and another group of 185 simulations that took commuting into account. We refer to these two simulation groups as non-commuting and commuting, respectively. The synthetic cases generated in each simulation were injected into the real ED data from 2005. Next, each of the 370 semi-synthetic outbreaks was supplied as input to the BARD detector and its average performance over the non-commuting outbreak group was compared with its average performance over the commuting outbreak group.

BARD's detection performance was measured through the time-to-detection--or timeliness--and its characterization performance was measured through the t, x, y, h, and q (absolute) errors. Each of the above six metrics is a function of the alarm threshold, i.e., the threshold of the alarm statistic employed by BARD detector to discriminate between the non-outbreak and outbreak situations. For a given threshold, an alarm is considered to be false if the alarm statistic exceeds the threshold prior to the outbreak onset. The false alarm rate (for a given alarm threshold) is the number of false alarms that occur per unit of time. We measured it as the number of false alarms per year.

We employed two methods to compare BARD's performance over the non-commuting group with its performance over the commuting group. The first method was Activity Monitoring Operating Characteristic (AMOC) analysis [11]. The AMOC curve of a performance metric (e.g., timeliness) plots the false alarm rate against that metric. The initial step in constructing an AMOC curve is to compute the alarm thresholds corresponding to various false alarm rates (in our case, 1 per year to 36 per year). We computed the alarm thresholds by running BARD detector on the baseline ED data for 2005. Next, these thresholds were employed to identify the pairs of values of the form <false alarm rate, performance metric> that constitute an AMOC curve. For each simulation group, we constructed a single AMOC curve per metric by plotting the false alarm rate against the mean value of the given metric over the whole simulation group.

The second method we employed to compare BARD's performance over the non-commuting group with its performance over the commuting group was statistical testing. Our null hypothesis asserted that for a given false alarm rate and performance metric, the mean of the given metric over the non-commuting group was equal to its mean over the commuting group. Testing was carried out through the paired t test.

We performed a second experiment designed to compare the performance of BARD and BARD-C detectors over a set of simulations generated with BARD-C simulator. As mentioned earlier, we refer to this experiment as the "detection study". Here, we briefly outline the design of this study; a full discussion is given in [10]. A total of 100 semi-synthetic outbreaks were generated using the BARD-C simulator. For all simulations, the quantity q was fixed at 0.5 kg. The release time t was chosen uniformly at random from the test period, i.e., the year 2005. The parameters x, y, and h were drawn from their prior distributions. Each semi-synthetic outbreak was supplied as input to the BARD and BARD-C detectors. The detectors were executed 42 times on each simulation in increments of 4 hours: the first execution began 4 hours after the release, the second execution 8 hours after the release, and so on.

The detection performance of both detectors was measured through timeliness. AMOC analysis and statistical testing were employed to compare the performance of the two detectors. The characterization performance of both detectors was measured through t, x, y, h, and q (absolute) errors. Each characterization metric was plotted against the time interval from the release of anthrax to the beginning of the detector's execution, which we call the time-to-execution. Finally, a Wilcoxon signed ranked test was performed with respect to each characterization error.

Figure 2 shows the results of our experimental verification of the correctness of simulation with commuting. This figure plots the estimated parameter values and the analytically predicted values for the 1991 tracts of our surveillance region. The solid line corresponds to the predicted values, while the points correspond to the values estimated from the simulated data. As shown, the predicted values are in good agreement with the values estimated from the simulated data, indicating that BARD-C simulator works as specified.

Figure 3 shows the daily mean number of cases (i.e., the epidemic curves) for the non-commuting and commuting simulation groups. In this figure, the horizontal axis shows the outbreak day, while the vertical axis shows the daily mean number of cases in the surveillance region. It can be observed that the shapes of these two curves are very similar. In particular, the two curves are seemingly identical at least up to the third outbreak day--the day when BARD typically detects a semi-synthetic outbreak. These two curves differ mainly around the peak of the outbreak, where the commuting curve is slightly more elevated than the non-commuting curve. These observations imply that any differences between BARD detector's performance over one simulation group and its performance over the second group could be attributed only to the change in the spatial distribution of cases caused by commuting and not to a consistent discrepancy between the sizes of outbreaks simulated with and without commuting.

With respect to q and h errors, we found that BARD's performance over the non-commuting group was almost indistinguishable from its performance over the commuting group. Unexpectedly, the posterior mean of the release parameters q and h did not appear to converge. We leave the investigation of this unexpected outcome as a topic for future research and in the remainder of this section we focus only on the remaining four metrics.

Figure 4 shows the results of the AMOC analysis for the timeliness, t error, x error, and y error. As shown, in Figure 4(a), BARD's mean timeliness over the commuting group is nearly one hour greater than its mean timeliness over the non-commuting group, at each false alarm rate. Likewise, BARD's mean t error (Figure 4(b)) over the commuting group is nearly one hour greater than its mean t error over the non-commuting group, at each false alarm rate. BARD's mean x and y errors (Figure 4(c-d)) on the commuting group are each nearly 2 km greater than the respective mean errors on the non-commuting group.

Finally, Table 2 gives for all false alarm rates between 1 per year and 12 per year the p-values of the test for equal means of BARD's timeliness, x error, y error, and h error over the two simulation groups. From the results shown in Table 2 it can be concluded that, at each shown false alarm rate, the mean timeliness, x error, and y error over the non-commuting group are statistically different from the respective means over the commuting group, at the 0.05 significance level. The same holds for the t error with the exception of three false alarm rates: 1, 4, and 6.

Some of the results appearing in this section (Figs 5 and 6) have been previously reported in [10]. We have reproduced these results here to contrast the simulation study with the detection study. This section also presents two new plots (Figs 7, 8) that correspond to our detection study but do not appear in [10]. Figure 5 shows the AMOC analysis for the timeliness of BARD and BARD-C detectors. As seen, BARD-C's timeliness is more than one hour smaller than BARD's timeliness at every false alarm rate. A paired t test showed that, for each false alarm rate, the mean timeliness of BARD-C was statistically different from the mean timeliness of BARD, at the 0.05 significance level.

Next, we turn to the characterization performance of BARD and BARD-C detectors. As in the simulation study, we found that the posterior means of the release parameters q and h did not appear to converge as the time to execution increased.

Hence, in the remainder of this section we focus on the x, y, and t errors. Figure 6 shows box plots and plots of the medians of the x error for BARD and BARD-C. In each plot the horizontal axis shows the time-to-execution. The samples from which the box plots were constructed correspond to the different simulations. Several conclusions can be derived from Figure 6. First, the x error of both BARD and BARD-C converges to relatively small values as the time-to-execution increases. Second, the errors of BARD and BARD-C appear to converge around the same time and nearly 10-16 hours after the detection. Third, the sampling distribution of the x error is right-skewed for both detectors and in almost each execution of BARD and BARD-C there is a small number of outliers. Fourth, after the convergence, the median of the x error for BARD-C is nearly 2 km smaller than the median of the x error for BARD.

Figure 7 plots the median of the y error for BARD and BARD-C. The post-convergence median of the y error for BARD-C was nearly 0.5 km smaller than the median for BARD. Finally, Figure 8 plots the median of the t error for BARD and BARD-C. The post-convergence difference of the medians for the t error was nearly 20 minutes.

Finally, we performed a Wilcoxon signed rank test whose null hypothesis asserted that BARD and BARD-C have identical locations of the distributions of x error, y error, and t error. The results of this test showed that after the convergence of the x error, the p-values of the test remain well below the 0.05 significance level (see [10]). Similar comments can be made for the y error. For the t error the convergence of the p values is not as clear as for the x and y errors.