Meeting radiation dosimetry capacity requirements of population-scale exposures by geostatistical sampling

Overview

Nuclear detonation scenarios created by HPAC were derived for 22 North American cities and surrounding regions. Initially, locations within likely fallout areas were sampled from the potentially affected population. In each real-world scenario, a set of randomly sampled locations downwind of the event epicenter are initially selected. In each scenario, we assume the epicenter location (“ground zero”), wind direction, intensity, and weather conditions. In this paper, we refer to these locations as samples, regardless of whether the radiation is quantified from either physical emissions or by quantifying absorbed biological effects. We assume that all measurements used to derive the plume are consistently obtained by the same approach. The contours of the HPAC radiation plumes were used to specify the dose levels and random locations downwind of the epicentre. In HPAC, the contours are created from the points at locations specifying the contour threshold, which have specified exposure levels. Simulated dosimetry measurements were created from these samples and used to populate geostatistical-derived radiation maps. In the simulation, we assume that map locations between the contour boundaries exhibit the same radiation levels of the outer boundary. This introduces a source of systematic error into dose estimation, since it is likely that that the actual exposures at these locations should be interpolated between the neighboring contours that circumscribe the specific sample location. Samples were then simulated, initially at randomly selected locations, and which were then used to further refine the radiation exposure plume. After the initial set of samples were analyzed and an initial draft plume was computed, subsequent sampling locations were specified by densification. Densification is the geostatistical procedure that targets and localizes an additional small cohort of sampling locations to mitigate uncertainty in environmental measurements. These kriging and densification processes are repeated for a limited number of iterations until the coverage area and the radiation level contours of the inferred plume stabilizes (i.e. additional sampling in the affected area does not significantly alter the geographic coverage of the plume or the estimates of absorbed radiation dose). Comparisons between independent replicates of the same scenarios using different, randomized, initial sample distributions evaluated the reliability of this approach. A detailed sequential protocol is available (http://dx.doi.org/10.17504/protocols.io.ba4nigve; “Protocol for Geostatistical Determination of Radiation Dosimetry Maps of Population-Scale Exposures”).

Acquisition and Processing of United States (US) Census Data for Geostatistical Analysis

The sizes of populations impacted by simulated nuclear events in each scenario were based on 2017 US Census Data of affected counties and subdivisions. US state and sub-division boundary files (in KML [Keyhole Markup Language] format) were retrieved from the US census bureau (https://www2.census.gov/geo/tiger/GENZ2016/shp/). Population data were downloaded from the US Census (“Incorporated Places and Minor Civil Divisions Datasets: Subcounty Resident Population Estimates: April 1, 2010 to July 1, 2017”), supplemented with additional population data from the American Fact Finder and Home Town Locator (both updated July 1, 2018).

Geostatistical analysis was implemented with the ArcGIS Runtime SDK (Software Development Kit) for Java and the ArcPy package. Geostatistical mapping was performed with customized Python scripts calling the ArcGIS software toolkit (ESRI, Redlands, California, United States). For import into the ArcGIS environment, the US state sub-division boundary files were first split by sub-division name and imported to ArcMap using the ‘KMLtoLayer’ tool. Exceptions handled states with multiple sub-divisions that share the same name, so that ArcMap did not also select unintended sub-divisions found outside of the radiation plume. Sub-division names with spaces or dashes were also modified by concatenating them prior to ‘KMLtoLayer’ conversion to avoid their unintended truncation.

Derivation of ground-truth HPAC radiation plumes

Radiological release scenarios were derived for 24 different US locations (in 22 cities) by the University of Ontario Institute of Technology (UOIT) Health Physics and Environmental Safety Research Group. Radiation plumes were simulated using HPAC v4.04 (developed by the Defense Threat Reduction Agency [DTRA]), which models the dispersion of radiation (as well as chemical and biological releases) assuming only the location of the epicenter, historical prevailing wind direction and speeds, and weather conditions. A typical ground truth plume comprised 4000 sample locations, each coincident with the contour boundaries for each radiation level. The majority of the nuclear incident scenarios did not include precipitation, except for New York and Washington D.C., where rain and snow were compared with normal plume conditions. The HPAC data are represented as a series of high-density data points at geospatial coordinates which define the shape of each radiation contour level of the plume. Topological exposure contours were computed in increments of 0.5Gy, across the 0.0-7.0Gy range, with 1.0Gy intervals shown. Contours were plotted on top of a US city map layer generated by the Humanitarian OpenStreetMap project (https://export.hotosm.org/en/v3/).

Radiation Plume Reconstruction with Iterative Kriging and Densification

The HPAC-generated plume for each scenario was visualized with the ArcMap software toolkit. The plume location was used to estimate affected population size. The US census-defined sub-divisions, which overlap and/or surround the plume of interest, were determined from overlap with their respective latitude and longitude boundaries. A Python script was written to select data points (or ‘samples’) at random locations within each Census sub-division using the ArcMap tool, ‘CreateRandomPoints_management’. These random samples, which corresponded to 0.1% of the population of each sub-division overlapping and surrounding the HPAC plume, simulated the locations for dose assessment. Randomization of samples was bootstrapped a minimum of 10 times per scenario to evaluate whether this procedure impacted the final plume derived by simulated physical or cytogenetic dosimetry. Results were exported using the ArcMap ‘ExportXYv_stats’ tool, then assigned radiation level values corresponding to the adjacent outer HPAC contour by a script comparing each sample with its location within the HPAC plume. The number of random samples generated for each scenario, and how many of those overlapped the HPAC plume, is available in S1A Table. Only a subset of the random samples shown overlapped the HPAC-derived plume; all remaining samples were considered to be unirradiated. The resulting output was then re-uploaded into ArcMap to perform kriging and densification.

We used kriging, a geostatistical interpolation technique which computes the spatial autocorrelation between data points (unlike deterministic interpolation techniques) to map predicted radiation levels by sampling around the epicentre of the event and surrounding region. Various kriging methods, available through the ArcMap extension “Geostatistical Analyst”, were evaluated for this study: Ordinary, Simple, Universal, and Empirical Bayesian kriging (EBK; each kriging method described in the S1 Methods). To determine which type of kriging generated the most accurate plume, random points representing 1% of the population of the Boston (N=617,594) and Cambridge (N=105,162) subdivisions were generated by ArcMap, of which 223 points overlapped the Boston HPAC plume (predicted dose >0Gy). While no plume could be derived from Simple kriging, the methods of Ordinary, Universal and EBK were successful (S1 Figure). Ultimately, EBK was used for all further analyses featured in this manuscript, as the plume generated best represented the expected HPAC plume (with a 57% overlap to HPAC plume vs. 35% overlap for Ordinary and Universal) and has the advantage of taking uncertainty measurements into account (S1 Figure). The presence of unirradiated data points adjacent to the plume was found to be crucial for accurate kriging, since these points served as boundaries for kriging. A high number of unirradiated (0Gy) samples can depress the range of the plume; therefore, these locations were restricted to the subdivisions immediately surrounding the irradiated region. The number of unirradiated samples for each scenario replicate is listed in S1A Table. We envision that testing could be greatly reduced by initially measuring background or low level physical radiation in population scale events by aerial surveys or targeted multiplex dosimetry. Dose reconstruction has been extensively modeled over varied geographical, topological and weather scenarios (21). In an actual event, environmental physical measurements at these locations without detectable radiation could substitute for dosimetry at this boundary, in order to focus attention on sampling within irradiated regions (22).

The “Densify Sampling Network” tool of the Geostatistical toolbox indicates lower confidence regions in the kriging-derived map, i.e. regions with highest variance specifying radiation dose (17). We applied this tool to limit results to regions that would most likely exceed a pre-defined radiation level threshold. In practice, the locations selected by densification would be used to direct first responders to new locations for subsequent rounds of data acquisition in order to improve the accuracy of the kriging-derived map. Using 2Gy as the critical threshold (selection criterion QUARTILE_THRESHOLD_UPPER option), densification on one plume identified a maximum of 200 new sampling locations. We assume that locations within the 0Gy envelope surrounding the plume do not have to be sampled in subsequent kriging iterations. Densification is a compute intensive step, requiring approximately 1 hour on a desktop with an Intel i7-4770 processor [3.4Ghz] and 16GB of RAM. Note that reducing the number of requested sampling locations decreases overall processing time. The Densify Sampling Network tool would sometimes select a sample at the same latitude and longitude between iterations. Furthermore, many densification-selected samples did not overlap the HPAC plume. As a consequence, the process often did not always yield 200 unique samples with values exceeding 0Gy. New sample data were assigned radiation values based upon their locations within the HPAC-generated plume, and kriging was performed on these and the original samples to generate another iteration of the inferred plume.

Simulated analyses of population-scale, nuclear radiation scenarios

A geostatistical workflow for managing iterative computation of the inferred radiation plume in population-scale radiation events when processing new samples (Figure 1). The purpose of simulating analyses of the radiation exposures was to determine accuracy for distinguishing clinically relevant exposure levels, and the number of irradiated individuals or locations needing dose measurements. This was based on an initial random set of location-based data points representing sample measurements and locations collected by first responders, followed by measurements at additional locations which were assigned by densification and kriging.

• Download figure
• Open in new tab

Figure 1. The workflow for handling a population-scale radiation event.

First responders collect dose measurements and coordinates of the tested individuals or locations. Measurements are then collected and mapped. A dose plume is generated, and densification is used to select locations with lower confidence radiation estimates for follow-on sampling. These steps are repeated until output plume dosimetry levels converge. The resultant plume can be used to differentiate locations associated with significant exposures (≥2Gy), from those below this or other thresholds.

EBK with default values for optional parameters was selected to predict the dose value which each location received (Figure 1), which generated a dosimetry-based radiation plume establishing the spatial boundaries of graduated exposure doses. Two consecutive plumes were quantitatively compared with a heat map matrix with overlapping radiation levels. Each row or column respectively corresponded to a dose range in the current or previous plume, and the dose ranges were sorted in the ascending order downward or rightward. Each cell indicates the average overlap percentage between two dose ranges in terms of area. Therefore, two identical plumes resulted in an identical identity matrix. The dissimilarity between the two consecutive plumes was quantified by the diagonal Bray-Curtis dissimilarity (BCD) and root-mean-square deviation (RMSD) between the heat map matrix and the identity matrix. Lower BCD and RMSD values indicate better fit (23,24). For this study, BCD is computed as: where n represents the 8 topological contours of the plumes (<1Gy, 1-2Gy, 2-3Gy, 3-4Gy, 4-5Gy, 5-6Gy, 6-7Gy, >7Gy), A_computed is the percent area overlap of a particular contour of each plume, and A_ideal represents an identical overlap between the plumes (A_ideal = 1). RMSD is computed as:

The diagonal values of the heat matrix represent the overlap of each radiation dose range (the topological contours) between compared plumes (see center heat matrices in Figure 2). The converging BCD and RMSD thresholds of 1/19 and 0.1 (respectively) were computed as equivalent to heat maps with a 90% overlap across each radiation dose level. This 90% stringency of overlap between consecutive iterations of geostatistical analysis was selected as a compromise to limit the number of samples while building a strong approximation of each HPAC plume (which themselves vary in shape and size). This has been a threshold indicated in prior geostatistical studies (25,26). The iterative workflow was discontinued when either metric dropped below the threshold.

• Download figure
• Open in new tab

Figure 2. The simulated analysis of the no-precipitation scenarios of urban New York [Replicate #4] (A) and Albany [Replicate #1] (B).

Each of the three left-most maps (first row) display sample locations and EBK-derived dosimetry plumes after several rounds of densification (iterations; note gradual improvement of plume from left to right), while the right-most map displays the HPAC plume. Plume colours (from red to purple) indicate the dose ranges of the region (from <1Gy [100cGy] to ≥7Gy [700 cGy]). The total number of samples exceeding 0Gy in each iteration are indicated (in parentheses). The matrices (second row) are used to compare two plumes at each dose range, where each value indicates the area overlap (background grey level is proportional to the percent overlap) of two regions of plume pair. The right-most matrix compares the final derived plume to the HPAC plume. Below the matrices, the diagonal BCD and RMSD values between successive iterations represented by each heat matrix are indicated. The green or red colors indicate that these two metrics exceeded or did not meet the threshold for convergence of the procedure, respectively. Subsequent iterations gradually improved the computed radiation dosimetry plume, with each adding a small number of new samples. Based on the census data, the converged plumes localized 80.3% and 75% of the locations (weighted for population density) for treatment-eligible radiation exposures in these scenarios.

The geostatistical workflow was used to analyze three replicates of each scenario, with each replicate initiated with a different set of random sample locations within the plume. After each iteration of kriging and densification, the resultant plumes were compared with the preceding versions (and the HPAC-based map) and presented using heat maps indicating overlap between different radiation levels. Variation between radiation levels at the same locations in successive plumes determined whether the distribution of the initially sampled locations affected the resulting maps. BCD and RMSD were computed when comparing derived plumes to both the HPAC plume (S1A Table) and all other replicates for the same scenario (S1B Table). Since the HPAC plume cannot be directly compared with the converged dosimetry plume because their data formats are incompatible, EBK was performed on all joint vertices associated with the mapped dosimetry results. This produced a topographic map that best approximated the HPAC plume. The percentage of individuals with ≥2Gy exposures that were correctly localized, assuming a uniform population distribution within each subdivision, was determined by converting the area overlap ratio with population. Area overlap is converted to estimated population affected by multiplying the area of plume overlap by the total population contained within a subdivision, divided by the area of said subdivision.

We further evaluated the proposed geostatistical approach by testing the method under 2 suboptimal sampling conditions. We mimicked improper sampling due to an inaccurate specification of wind direction for the Albany NY scenario by providing samples which partially deviated from the affected census sub-divisions (e.g. under sampling the region overlapping the plume while oversampling a non-irradiated region; S1A Table). The angle was determined by finding the distance of mean latitude or longitude of all random points and computing an angle against the direction of the HPAC plume. To simulate variation in the precision of the radiation dose measurements, we applied random dose errors to HPAC generated radiation doses (to a maximum value of either ± 0.5Gy or ± 1.0Gy) to all initial and densification-selected samples. A representative subset of the previously analyzed radiation release scenarios was assessed for plume accuracy, and the numbers of kriging/densification iterations and samples required for convergence by iterative kriging and densification (S2 Table). These included Birmingham AL (Replicate #2), Boston MA (Replicate #1), Chicago IL (Replicate #1), Columbia SC (Replicate #1), and Columbus OH (Replicate #1 and 2).

Development of Simulated Radiation Plumes

We simulated the analysis of 28 population-scale, 10 megaton yield nuclear detonation scenarios with dosimetry data for radiation exposures corresponding to HPAC-derived dose estimates. These results were used to derive plumes of absorbed radiation exposures. The simulations included 22 in urban/high-density populated regions of the United States, and 2 rural/low-density scenarios. The process of kriging and densification was repeated until the difference between the current and previous plume reach a minimum threshold (Figure 2). The number of samples with radiation exposure used to generate the final iteration of a plume would vary between different scenarios (ranging from 58 and 347 samples within the irradiated plume; starting from an initial sample set representing 0.1% of the population; Table 1). For less densely populated regions, initial sampling was also performed at higher population densities (0.2 and 1.0%). The number of samples necessary to reach this stopping point was sometimes variable between replicates, and this variability was inversely proportional to the overall population density of the region (from US census). Cities with a high population density (>10,000 per square mile) had a 16.0% average coefficient of variation (CV) of the number of samples necessary to reach convergence, whereas the scenarios within low-density regions (<10,000 per square mile) had an average CV of 21.4% (see S3 Table for all CVs). In general, more samples would lead to a plume which better resembles HPAC (Table 1). There are exceptions where equivalent accuracy is obtained with fewer overall points (e.g. Buffalo NY; Table 1) which implies that the spatial distribution of initial set of random samples can also influence the accuracy of the converged plume.

Success in reconstructing simulated plumes from dosimetry data was based on the accuracy of predicting irradiated samples ≥ 2Gy (as defined by the HPAC plume). This radiation threshold was selected based on US government recommendations for eligibility of clinical treatment of acute radiation syndrome by cytokine therapies (27). Results are reported on the basis of the accuracy in distinguishing samples exposed to this threshold or higher from false positives or negatives. In these simulations, we have determined that false positive samples (i.e. ≤ 2Gy account for [on average] 0.5−2% of those estimated to be exposed at ≥2Gy. On average, the accuracy of the simulated plumes to predict exposures above this threshold (relative to the HPAC radiation plume) was 69.8%, ranging from 52.8% (Alexandria VA) to 90.7% (Burlington VT). The accuracy range for replicates from the same geographic region was also consistent, i.e. within 10% difference between each.

Two of the 28 scenarios exhibited outlier replicates, in which accuracy was impacted by the random sampling procedure and local population densities (Columbia SC replicate #2 [30.1%] and Columbus OH replicate #2 [44.1%]; Table 1 and S1A Table). Additional replicates were analyzed for these scenarios. While both additional replicate sets for the Columbus OH scenario closely resembled the first replicate (67-68% accuracy), one of the Columbia SC replicates continued to perform poorly (27.6% [replicate #5]). Further investigation suggested that this result might not correlate with the number of initial sampled locations within the irradiated region (12-16 samples exceeding 0Gy). Rather, it was apparent that the initial sets of sample data for these particular replicates exhibited sparse coverage over large regions of the HPAC radiation plume. Including subsequent densification steps did not always identify additional sample locations that resulted in a contiguous 2Gy plume similar to the entire HPAC plume. This issue was addressed by increasing the population fraction that was initially sampled. By combining the initial samples of the two underperforming Columbia SC replicates (#2 and #5; now representing 0.2% of the population rather than 0.1%), a plume was generated with accuracy equivalent to the best performing Columbia replicate (65.1% accuracy in 6 iterations [replicate #6]; Table 1). An independently selected random set of samples representing 0.2% of the population led to a comparable (albeit slightly lower) accuracy rate using fewer kriging and densification cycles (58.1% in 3 iterations [replicate #7]). A further increase in sampling rate to 1.0% of the population improved performance relative to the best Columbia SC replicate (+3.8% more accurate; replicate #8). Thus, in scenarios with low population densities, the population fraction sampled will need to be increased to construct a contiguous dosimetry map that resembles the actual radiation plume.

In rare cases, low population scenarios would fail to yield a plume in the initial step, regardless of kriging method used. When this occurred, the iteration methods could not progress to the following densification steps, and plume development was halted. The numbers of samples which overlapped the HPAC plume were found to be extremely low in these cases, ranging from 0-2 samples in total. There are instances in which 2 irradiated samples were adequate to progress plume development (e.g. Cincinnati urban sampling #2 [Table 1]). In these cases, a greater number of iterations were required for plume development to fulfill the established stopping criteria (N=8 iterations for the Cincinnati example). Therefore, while the minimum number of irradiated samples can be extremely low, development of accurate dosimetry plumes in such cases could require unacceptable levels of exposure for first responders.

We also compared coverage of derived plumes at a higher radiation level (≥ 3Gy), since these were expected to be more compact than the ≥ 2Gy contour, more densely populated, and more likely to result in severe clinical symptoms. Based on the degree of overlap between the derived and HPAC plumes (Table 1), 19 of 28 scenarios were more accurate at the ≥ 3Gy contour than at ≥ 2Gy level for at least one replicate (greater accuracy was seen for all replicates of 13 scenarios). In particular, the Alexandria VA, Boston MA and Buffalo NY cases respectively exhibited an average of +25.7%, +21.0% and +16.7% higher accuracies at the 3Gy threshold. In several instances, the performance of replicates of the same scenario varied (e.g Charleston SC replicates #2 (+15.2%) and #3 (-4.9%) in comparison with the 2Gy threshold). Derivation of plumes for the Charleston scenario has been challenging due to differences in population density for different initial sampling locations, for example, replicate #5. However, this issue was mitigated by sampling at higher densities (replicates #6-8), for example, the 1.0% population density replicate, in which the 3Gy contour accuracy exceeded that of the 2Gy contour (+2.9%).

Distinct Initial Sampling Produced Similar Plumes for the Same Scenarios

In each case, determination of the distribution of radiation exposure in each scenario began with an initial set of sampling locations within the region of the simulated radiation release scenario. The majority of these locations did not overlap with the HPAC plume and have therefore been modelled as unirradiated samples. As a result, the initial number of irradiated samples can vary significantly among different scenarios, based on the population of the region and overall plume size. Thus, the number of iterations of densification required to maximize the accuracy of the simulated dose plumes varied significantly between replicates (between 3 and 10 iterations per replicate; Table 1).

The initial sets of samples were randomly assigned locations within each census sub-division by ArcMap. We first determined whether this process affected accurate derivation of the plume as a consequence of the variable locations of these initial conditions. Replicate dissimilarity was quantified as BCD and RMSD (S1B Table) to determine the extent to which this source of variation could impact the accuracy of the derived radiation plumes. The BCD between plume replicates were less than 0.3 (or a >70% similarity between them [23]) for nearly all scenarios (S1B Table). New York NY and Cincinnati OH show > 90% similarity (by Bray-Curtis) across all replicates, while Albany NY, Boston MA, Birmingham AL, Camden NJ, Cleveland OH, Detroit MI, Grand Rapids MI, Philadelphia PA and Washington D.C. (urban scenario) have at least two replicates with a >90% similarity to each other.

Increasing Stringency of Overlap between Consecutive Iterations of Geostatistically-Derived Radiation Plumes

To expedite dose estimation for triage management, sampling was discontinued when the estimated populations in areas covered consecutive plume iterations were within 90% of each other by similarity. We investigated whether further cycles of geostatistical processing would stabilize or improve the derived plumes relative to the HPAC standard. We therefore mandated a more stringent threshold for consistent coverage of iterative plumes, i.e. from 90% to 99%, and determined whether the additional iterations would more closely resemble the gold standard for one replicate of each scenario (S4 Table). On average, an additional 3-4 iterations were necessary to reach this threshold, with an average of 65 additional sampling locations. In 3 scenarios, the plumes improved significantly: rural New York (+9.6%); Alexandria VA (+7.7%); and Columbus OH (+5.1%). In the rural New York scenario, the 1Gy and 2Gy contours are significantly expanded with the inclusion of additional data points (S2 Figure). However, altering the stopping criteria did not significantly improve plume accuracy in most replicates tested (+1.2% improvement, on average). Frequently, the plume would prematurely stabilize once densification failed to reveal new sampling locations. The addition of further iterations did not always improve accuracy at the 2Gy threshold. For example, addition of 43 new unique sampling locations slightly reduced the accuracy (-2.5%) of one Boston MA replicate (S4 Table). While the 2Gy contour in this instance did expand, new 1.0 and 1.5Gy samples introduced by densification made the 2Gy contour discontinuous, explaining the reduction in accuracy (S3 Figure). Contours at higher radiation levels for this plume were unaffected, and were contiguous, however.

Impact of Weather Effects on Derived Plumes

The effects of adverse and normal weather conditions (i.e. no precipitation, rain and snow) were compared for both urban and rural New York and Washington DC scenarios (Table 1 and S1A Table). The HPAC software was used to simulate the effects of rain and snow conditions in these scenarios at the same nuclear yields. The resulting HPAC plumes closely resembled the morphology, extent, and wind vector of the no precipitation scenario. At the >2Gy threshold, both of the final New York plumes with precipitation exhibited less accurate geographic coverage (72% rain and 74% snow) than under normal conditions (76-81% without weather effects; S1B Table). By contrast, the Washington D.C. urban snow scenario significantly outperformed both no precipitation and rain scenarios (82% accuracy versus 66-68% for all other weather conditions; S1B Table). These plumes were found to have a 70-76% similarity to the HPAC plume (by Bray-Curtis), similar to that of the same scenarios in the absence of weather effects (S1A Table). In these nuclear event scenarios, adverse weather events did not confound or impact the accuracy of population scale geographic dosimetry using the same approach described for normal weather conditions.

Inferred Radiation Exposures Under Suboptimal Sampling Conditions

If the information about wind direction and location of the epicenter of a nuclear event is limited and/or inaccurate, it might be expected to lead to improper sampling. Initial samples for the Albany NY scenario were used as input, under sampling the Albany sub-division (which overlaps the majority of the plume) while oversampling its neighboring sub-division, Colonie, which is due north of the epicenter of the simulated nuclear event. The bearing angle of the Albany NY scenario is N 51.9° W. For the initial sample data, the populations representing Albany and Colonie were shifted from equal 0.1% proportions of their respective populations to unequal distributions of either 0.05:0.2% or 0.01:1.0%. The 0.05:0.2% ratio simulates a wind measurement error of 29.1° north (or N 22.8° W) relative to the actual wind direction of the HPAC plume. The 0.01:1.0% ratio corresponds to a deviation of 40.9° north (or N 11.0° W). Despite this initial sampling error, inferred radiation plumes comparable to the correct plume were obtained. The accuracies of the >2Gy threshold maps were, respectively, 74.2% for the 0.05:0.2% population ratio and 74.5% for the 0.01:1.0% ratio, which compares well with the correctly sampled initial map of 75.0-78.4%. Reconstruction of the geostatistical dosimetry map required 4 iterations with sampling error, which was only one more cycle of densification-sample procurement and kriging than the map generated from the original (correctly sampled) Albany scenario. This demonstrates the ability the robustness of the densification step to select relevant sample points at subsequent iterations of kriging and densification that result in a correct dosimetry map. It appears that, in some scenarios, reconstruction of the original HPAC map can largely succeed with fewer irradiated samples at the expense of a single additional sampling cycle (e.g. iteration) for convergence.

The original set of simulations presented do not account for error in dose estimates, which can vary considerably based on the methodology used (28). The impact of random variation on measured dose was examined by deriving plumes with modified sample exposures with added or subtracted random errors. Random errors in radiation exposure (either ±0.5Gy or ±1.0Gy) were applied to both the initial samples and those obtained from subsequent densification steps. Maximum deviations were designed to represent confidence values in physical and/or biodosimetry methods, including physical dose estimation error (29), dicentric chromosome analysis (DCA; ±0.5Gy) or cytokinesis-blocked micronucleus assays (±1.0Gy), which exhibit higher variance in estimating dose compared to DCA (28, 30-31). The introduction of such measurement errors led to plumes that were generally less accurate, required a greater number of samples, and additional iterations of kriging and densification to achieve convergence (S2 Table). The degree of error in this analysis was also correlated with decreased accuracy and increased processing time (±0.5Gy error-derived plumes were more accurate than ±1.0Gy error-derived plumes in nearly all scenarios; S2 Table). We also note that dose modifications (±1.0Gy) for the Columbia SC scenario #1 completely failed to achieve complete plume coverage, and as a consequence, did not successfully derive a plume with accurate radiation exposure levels >2Gy.

Note added in proof

We recognize that the methods described here could be useful for other applications, for example as a tool for late stage epidemiology of a highly communicable infectious pathogen which disperses across a geographic region. We envision applying geostatistical inference to fill a gap in assessment after community transmission, when patient tracking is less useful as an effective mitigation strategy to identify sources of disease transmission. We suggest that densification step, which selects locations based on the variance levels of demonstrably exposed individuals might prove useful in identifying emerging hot spot locations. For example, after population isolation in a pandemic, at the point where the public health monitoring is moving into a management phase, the identification of such locations could be valuable in prevention of large waves that could overwhelm regional healthcare capacity.