Study area

The whole of Mainland China was selected as the study area. Mainland China is defined as all regions under the direct jurisdiction of the People's Republic of China (PRC), excluding Hong Kong and Taiwan. All autochthonous H7N9 cases have so far been isolated from Mainland China, whereas Hong Kong and Taiwan have not yet reported autochthonous human or animal H7N9 cases.

The study area is approximately 9,596,960 km², and is extremely diverse in terms of physical and human geography. Generally, the eastern and southern half of China consists of fertile, low-lying land and is more amenable to agricultural activity and human habitation, whereas the western and northern half of China largely consists of uninhabitable desserts, mountains and high plateaus, which are much less amenable to any anthropogenic activity.

All analyses were conducted in an unprojected, geographic coordinate system WGS84 (World Geodetic Survey 1984) using a spatial resolution of approximately 1km² (0.0083 decimal degrees)–the study area is made up of 40 046 238 of these individual cells.

Data

Disease occurrence data was collected from the EMPES-i georeferenced disease data repository compiled by the Food and Agricultural Organization (FAO) (http://empres-i.fao.org/eipws3g/, last accessed 21 September 2015). The Chinese Ministry of Agriculture avian influenza surveillance reports (www.syj.moa.gov.cn, last accessed September 2015) and official H7N9 World Organization of Animal Health (OIE) reports (www.oie.int, last accessed September 2015) were also used to extract the names of locations (e.g. LBMs, poultry farms, parks or wetlands) reported to have a positive sample of H5N1 or H7N9.

Maxent is designed to use only coordinate data which accurately reflect the true location of the case or outbreak. Hence, to maintain the integrity of our Maxent SDMs, we separated our disease occurrence data into two data sets: case records where the coordinates provided were associated with an acceptable level of accuracy and precision of at least 2 decimal places (to reflect a spatial resolution of approximately 1km²) were assigned to an “exact” data set; whereas case records with low or unknown location quality and low precision were assigned to an “unexact” data set.

EMPRES-i records included information on the locality, and an evaluation of the quality of location coordinates. All records where the precision of coordinates were labelled either administrative region centroids, or unknown were classified as “unexact”–only records where the precision of coordinates were specifically assigned “exact” were included in our “exact” data set.

For ministry and OIE records, latitude and longitude coordinate data were retrieved through searching location names in Chinese and English (using Baidu and Google mapping services, respectively). Coordinate data was assigned “exact” only when location names matched registered and georeferenced market (or farm/park) locations on Baidu and Google. Locations were matched by name, and address details when available (street name, district name, county name, city and province). Some location names provided by reports did not specify a distinctive name (only details such as district, county and province were provided)–these locations were not included in our exact data set.

Overall, we obtained coordinates for 267 H5N1 records, and 1289 H7N9 records. We obtained “exact” coordinates for 52 H5N1 records and 69 H7N9 records, and we have provided these data sets in a downloadable excel format (S1 and S2 Files). We summarised our H5N1 and H7N9 exact data sets by year of outbreak or human case onset, type of host the virus was isolated from (e.g. human, animal or environmental sample), and locality type (e.g. market, farm, village). These are listed in S1 and S2 Tables, and plotted in S1 and S2 Figs. We refer to each record as a ‘case’ throughout this paper.

Environmental variables for SDM construction were obtained from the WorldClim database (www.wordlclim.org, last accessed 12 October 2015) and are listed in Table 1. Each layer was available at a resolution of 30 arc-seconds (0.008333° x 0.008333°). Elevation data was downloaded from the Shuttle Radar Topography Mission (SRTM) 90m Digital Elevation Database v4.1 (www.cgiar-csi.org, last accessed 12 October 2015). This data was initially available at a resolution of 1km². The data was used to derive the slope, aspect and composite topographic index using the Spatial Analyst extension of ArcMap 10.2, as described in Minnesota Department of Agriculture [41]. Topographical layers were resampled to align to the WorldClim layers.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Environmental layers for Species Distribution Model (SDM) construction.

https://doi.org/10.1371/journal.pone.0174980.t001

We obtained human population data from the LandScan (2014)^™ High Resolution global Population Data Set (copyrighted by UT-Battelle, LLC, operator of Oak Ridge National Laboratory under Contract No. DE-AC05-00OR22725 with the United States Department of Energy), available from (http://web.ornl.gov/sci/landscan/), last accessed October 2015. This data was originally available in a resolution of 30 arc-seconds and was resampled to align to the extent of WorldClim layers. Values at each cell represent the average or ambient human population distribution.

Domestic chicken population data were obtained from Livestock Geo-Wiki (http://www.livestock.geo-wiki.org/, last accessed October 2015), and details on data set construction were provided in Robinson et al. [42]. The authors provided separate domestic chicken population rasters: (i) for chickens occupying intensive production systems (typically high density poultry holdings located nearby urban areas) and (ii) for chickens in extensive production systems (typically low density poultry farms which have locations amenable to prediction based on agricultural/land properties). Data was provided at a resolution of 1km² and was resampled to align to the extent of WorldClim layers. Values for each cell represented the number of birds per km². A summary of the data are provided in S3 Table.

LBM data were requested from authors of previously published SDMs [27,30]. Base maps of Chinese administrative regions (primary and secondary) were obtained from the GADM database of Global Administrative Areas (http://www.gadm.org/, last accessed October 2015).

Species distribution model construction

The Maxent software package (version 3.3.3k, available from https://www.cs.princeton.edu/~schapire/maxent/) [43] was used to construct 8 SDMs. Maxent software [43] is primarily used in ecological sciences to model species distributions although the algorithm has previously been used in regional and global studies of avian influenza [17,18,20,22,23,25]. Maxent uses a machine learning algorithm to produce an estimate of the ecological suitability of a species using a set of environmental layers and a set of accurate case coordinates [43]. It does this by fitting a probability distribution to a set of cells in a study region—it finds the probability distribution of maximum entropy (i.e. the most spread out distribution) subject to a set of environmental constraints that represent the species distribution [43]. For each cell in the study area, a relative environmental suitability for species presence is calculated (referred to as suitability values). Maxent consistently produces relatively accurate results, even when only small sample sizes are available (generally considered as n <100) [44–46]. A major advantage of Maxent is that it only requires presence data, making this method more advantageous over other genetic algorithms and regression methods. Instead of prediction of presence or absence, Maxent estimates the relative environmental suitability for presence of a species, which also allow for finer estimates. Additionally, most publically available epidemiological data consist of records of human infections and animal outbreaks ascertained from passive surveillance, whilst systematic active surveillance using formal random sampling methodologies are not routinely performed with regards to AIVs—and details of negative cases or outbreaks (i.e. ‘controls’) are not made publically available.

All environmental layers for model construction are listed in Table 1. The environmental data was included because climate factors are known to influence AIV occurrence. Examples include modulating survival of AIV in the environment, and modifying animal host behaviour (e.g. seasonal migration of ducks and geese) or human behaviour, or host susceptibility to infection. Observational studies and spatial-temporal analyses have shown AIV infections display seasonality patterns in humans and animals [47,48], solar cycles for example have been associated with pandemic influenza activity [49] and animal transmission studies of human influenza viruses have shown that changes in temperature and humidity affect transmission potential and virus binding [50]. Previous studies have also shown decreasing or increasing relative humidity influences virus stability (discussed in [50]). Ambient temperature can impact virus transmission due to physiological changes in the hosts which increase susceptibility to infections [51], for example, in wild bird species which carry LPAI AIVs transmission can increase dramatically during migratory seasons, as birds congregate in large numbers in water bodies, and their immunity levels are likely to be low due to the physical exertion required during migration. In addition, during breeding seasons the number of immunological naïve birds increases substantially around this time, thus increasing the amounts of AIV shed in the population [52].

Temperature and precipitation can influence survivability of AIVs in the environment: H5N1 can remain infectious for 1 day at 37°C, 5 days at 24°C and 8 weeks at 4°C in dry and wet faeces, and laboratory experiments have shown H7N9 virus stability changes with varying temperatures [53]. Also, weather variations may influence human activities—for example more people go to the markets on a sunny day, which may increase virus transmission [54].

Topographical variables (elevation, slope, aspect and composite topographic index) were also used for model construction. Topographical variables, particularly elevation and slope, are commonly used in ENMs, and have been included in previous AIV models [55]. The reasons for their inclusion were (i) they are widely available for most countries, and (ii) in some studies, they have been found to contribute significantly to AIV occurrence, possibly as a surrogate indicator for some unknown variable related to AIV risk [56]. Elevation is thought to relate to land coverage—for example, areas of high elevation are dominated by forests, whilst flat areas (such as plains, deltas and coastal regions) are most likely be used for agriculture or urban development [56]. Slope and aspect may similarly influence AIV distribution through vegetation and land coverage. CTI is a measure for the tendency of water to pool, and may indicate areas suitable for rice cropping or duck farming which have been associated with H5N1 occurrence in Asia [57]. A discussion of CTI’s influence on H5N1 transmission in Vietnam is provided in Saksena et al. [58]. Aspect values indicate the compass direction which slope faces, and can be derived from slope values [21,59].

A variety of landcover and satellite data sets exist, and have been used in previous SDM models [20,22,23,25,60]. The reasoning behind their use is similar to the reasoning behind use of climate variables (they are able to capture seasonality and thus waterfowl migration patterns). For our model, landcover variables and satellite data were not used as our SDMs should already account for landcover and seasonality as they were created using precipitation, temperature and topological variables.

Maxent was run using Auto features. The test:training ratio was set to 40:60 as recommended by Phillips and Dudík [43]. Logistic output formats were selected as recommended in Phillips and Dudík [43]. In logistic output formats, each cell is assigned a probability (a value between 0 and 1) which can be interpreted as the estimate of the probability of species presence or as relative suitability. No thresholds were selected. Prevalence was set to 0.5 (default value) as we were unable to obtain reliable H5N1 and H7N9 prevalence or detection rates for all of mainland China at high resolution. For H7N9, there is more uncertainty around prevalence rates as this virus is asymptomatic in birds [3]. For each study, 100 models were created and relevant outputs were averaged.

For each virus subtype we created 4 different Maxent models which varied by number of variables and study area (model parameter settings remained the same): Models 1–2 were built using the full set of environmental variables and the full study area; models 3–4 were built using a subset of variables and the full study area; models 5–6 were built using the full set of variables and a subset of the study area; models 7–8 were built using a subset of variables and a subset of the study area.

For models 3,4,7 and 8, which included a subset of the environmental variables, 11 of the 23 environmental variables were finally chosen to be included in the subset (similar to Sarkar et al. [16]): Annual mean temperature, mean diurnal range, maximum temperature of warmest month, minimum temperature of coldest month, annual precipitation, precipitation of wettest month, precipitation of driest month, plus the four topological variables (see Table 1).

For models 5–8, which were created using only a subset of the original geographic study area, we selected only primary administrative regions (provinces, municipalities, autonomous regions) which had at least one H5N1 or H7N9 case (see S3 Fig), with the following exceptions: (i) the provinces of Shaanxi, Sichuan and the Chongqing municipality don’t contain any H5N1 and H7N9 cases from our exact data set however we include these provinces as they are geographically contiguous with selected provinces; (ii) Xinjiang Uygur autonomous region covers an expansive area, however contained only one H7N9 case and was hence excluded from the analysis. In total, 22 of 31 primary administrative regions were selected as the subset study area.

We evaluated each of the 8 models based on the receiver operating characteristic (ROC) and area under the curve (AUC) values. The ROC plot is a plot of sensitivity and 1–specificity [61]. In Maxent, sensitivity refers to how correctly the model has predicted presence, and specificity corresponds to a measure of correctly predicted absences. A model with good fit will produce a ROC curve which maximises sensitivity for low values for 1 –specificity or false positive rate. The AUC quantifies the significance of the ROC, with AUC values of 0.5 indicating the model is no better than random, and a value of 1.0 indicating perfect fit. For each model, we averaged the testing and training AUC, and standard deviation of testing AUC, over 100 replicates.

Overfitting can be an issue with maximum entropy methods, particularly when data is sparse, and can result in distributions clustering around the case data. The Maxent software has a relaxation component (regularization) which is designed to counteract correlations. Overfitting is more likely to occur when large numbers of environmental layers are used due to correlations between variables. A sign of overfitting is when a model performs very well (i.e. has high AUC results on the training data set), however performs poorly for unseen data (i.e. for the test data set). We adopt the same methodology as Sarkar et al. [16] to test whether models built with the full set of environmental variables (n = 23), were suffering from overfitting. First, the difference between test and training AUC scores was calculated for each set of 100 replicate models. Normality was tested using the Shapiro test and as not all data were normally distributed, the non-parametric Mann-Whitney-Wilcoxon test was used to decide whether the distribution of AUC differences of models with using all variables (n = 23), compared to models with fewer variables (n = 11), were identical. All statistical computations were performed using R.

The way environmental factors drive AIV occurrence differs according to region [22,23]. Having such a large study area (i.e. the whole of mainland China) means there is a large selection of background points, and as there are only a sparse number of disease cases this may influence reliability of predictions [62–64]. If case data only fall into a subset of the study area, Phillips [63] recommends a solution is to only draw background points from this subset to improve model performance. Models 5–8 were developed using a subset of the original study area of mainland China, and compared to models 1–4 in terms of differences in distributions, AUC and overfitting.

Risk model construction

SDMs alone can only provide estimates of the probability of virus presence. In order to estimate the likelihood of zoonotic transmission, a measure of transmission efficiency must be accounted for. As described in Hill et al. [33], for zoonotic transmission to take place, a susceptible human must be within range of an infected animal (or an animal environment contaminated by high viral shedding). The level of opportunity for human exposure from the virus is proportional to the product of the number of infected animals and the number of susceptible humans [65]. To estimate the risk of circulating viruses to cause human infection, the SDM outputs for H5N1 and H7N9 are combined with human and animal population density. We use domestic chickens as the representative animal host, as these animals make up the highest proportion of China’s poultry sector [66], are the most commonly identified animal host of H5N1 and H7N9 [67], and virus shedding occurs at a higher rate in chickens compared to other avian species [5].

We modify a formal risk assessment methodology described in Sarkar et al. [16]. Risk models were constructed for each virus subtype independently, whereby a value between 0 and 1 was computed for each cell representing the relative risk of a human infection of H5N1 or H7N9 compared to other cells. These models combine ecological factors with demographic and agricultural factors known to modulate AIV transmission, and disregard variations of risk from human interventions such as those described in the introduction (animal biosecurity, vaccination, LBM closures etc).

For each virus subtype, we use a simple multiplicative model for computing the risk, r_k, for each cell k in the study area as shown in Eq (1).

(1)

This equation provides a measure of the risk posed to humans in cell k in regards to the likelihood of becoming infected with H5N1 or H7N9. The variable, p_k, is the estimated prevalence of H5N1 or H7N9 in cell k (relative to other cells in the landscape) from the Maxent output. The variables, h_k, and c_k, represent human and chicken population, respectively, in cell k. Computed r_k values were normalised, as shown in Eq (2), by dividing by the highest value computed over all cells, k, in the study area. The final result, R_k, is the relative risk of a human infection of H5N1 or H7N9 posed to each cell: (2)

A log transformation was used for the human population density and intensive and extensive chicken population density variables, due to their highly skewed nature. Intensive and extensive chicken population densities were summed, and treated as a single variable, c_k, in the model. Both human and chicken population densities were normalized to the highest value over all cells, k, in the study area. A frequency histogram showing the distribution of each variable for the set of exact case locations for each virus subtype is provided in S4 Fig.

Based on the exact case data (S4 Fig) the majority of cases for both H5N1 and H7N9 fell in regions with medium human density and medium poultry densities, and less cases fell in low and high density regions. For chicken density this is possibly because H5N1 is more likely to occur in chicken farms with poor hygiene and biosecurity practices, regardless of the size or density of the farm. There is extreme variation in terms of biosecurity, hygiene and farming practices among chicken farms in Asia—and such variation is difficult to account for using the chicken density data sets that are currently available [55]. Stevens et al. [31] addressed this limitation by assuming a quadratic relationship between H5N1 risk and chicken density, with the highest risk areas being those with a medium density of chickens (between 500–5000 heads/km²). Their reasoning was based on a review of the literature [55], which found that most large commercial chicken farms (where chicken densities can be very high) often implement appropriate biosecurity regulations and are expected to have little risk of H5N1 occurrence. In contrast, medium-sized poultry farms (representing semi-commercial and backyard producers) are often the ones that are unregulated and have less stringent infection control. Subsistence farmers with very small flock densities are considered to have negligible risk due to insufficient hosts to sustain an outbreak. For human density the explanation is less obvious. One possibility is to do with the changing food consumption practices, which are more evident in urban areas—a literature review of food consumption practices in China found that people are moving away from the traditional practice of purchasing fresh market produce and preparing meals at home, to now eating outside of the home (e.g. in restaurants) and eating processed foods [68]. Hence, we assume that transmission risk is lower in low human density areas and extremely high human density areas. There are however, no large scale surveys of consumer preference specifically for frozen or live poultry in China.

Given these observed relationships between the exact case locations and respective spatial variables, the human and chicken density variables were reclassified to a scale of 0 to 1 based on a simple Gaussian fuzzy membership function prior to being used in the risk function. This transformation is illustrated in Eq 3. (3) Let x be the transformed (log and normalized) human or chicken population variable in cell k, and let the parameter, c, represent the ideal membership for the set and any values where x = c is assigned the value of 1. As x values move away from the midpoint, their membership to the set gradually decreases. When x values are too distant from the ideal definition, they are no longer considered to be in the set and are assigned zeroes. The parameter,σ, represents the spread or width of the Gaussian function (Eq 3). Values of c and σ were selected for each variable based on the distribution (i.e. mean and standard deviation) of each variables in the set of locations where exact H5N1 and H7N9 cases lie. For H5N1 we used c = 0.532 and σ = 0.183 for chicken density, and c = 0.619 and σ = 0.198 for human density. For H7N9 we used c = 0.601 and σ = 0.103 for chicken density, and c = 0.689 and σ = 0.146 for human density. All spatial analyses were performed in ArcMap 10.2 [41].

Species distribution models

Based on model evaluation results (summarised in Table 2 and Fig 1), we chose SDMs 3 and 4 (see Fig 2) to be included in our risk model. Final SDMs 3 and 4 are shown in Fig 2, and the remaining models are shown in S5–S7 Figs. Prediction capacity was high for 6 of 8 models (AUC > 0.90). H7N9 SDM models (SDM 2,4,6 and 8) consistently performed better than H5N1 in terms of AUC. SDM models 1–2, suffered from overfitting based on Mann-Whitney-Wilcoxon tests (P-values turned out to be less than the significance level, hence the null hypothesis is rejected and AUC differences are non-identical). However overfitting was not an issue for SDMs 5–6. Reducing the study area (SDMs 5–8) did not make considerable differences in terms of suitability distribution (see S6 and S7 Figs), or AUC results (see Table 2). The term “suitability distribution” refers to the distribution of these suitability values in geographic space. For each subtype, the suitability for occurrence is displayed on a continuous scale from 0 (least suitable) to 1 (highly suitable). On visual appraisal, all 4 H5N1 suitability distributions were similar to each other; likewise, all 4 H7N9 suitability distributions were similar. H5N1 suitability distribution was distinctly more geographically diffuse across south-eastern China, whereas H7N9 suitability distribution was clustered densely around distinct areas (Zhejiang and Guangdong provinces). For H5N1, 5.8% (3 of 52) exact cases fell into low suitability areas (p_k <0.25), and for H7N9, 8.7% (6 of 69) exact cases fell into the low suitability areas (p_k<0.25).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Species Distribution Model (SDMs) evaluation results.

Top left panel shows boxplots of the area under the curve (AUC) for test data. The top right panel shows boxplots of AUC for training data. The bottom left panel shows boxplots of standard deviation of test data. The bottom right panel shows boxplots of test and training AUC differences.

https://doi.org/10.1371/journal.pone.0174980.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Species Distribution Models (SDMs) built using Maxent.

The first panel shows H5N1 (SDM 3) and the second panel shows H7N9 (SDM 4). Suitability values for each cell (approximately 1km²) are represented on a continuous scale of low (light grey) to high (dark grey). SDMs were built using Maxent software version 3.3.3k (available from https://www.cs.princeton.edu/~schapire/maxent/). SDMs were developed using environmental variables, created using data from: the WorldClim database (www.wordlclim.org), the Shuttle Radar Topography Mission (SRTM) 90m Digital Elevation Database v4.1 (www.cgiar-csi.org). Data sources used to obtain the case locations to build SDMs include: the Food and Agricultural Organization (FAO) (http://empres-i.fao.org/eipws3g/), the Chinese Ministry of Agriculture Avian Influenza Surveillance Reports (www.syj.moa.gov.cn), the World Organization of Animal Health (OIE) reports (www.oie.int). Base maps were obtained from the GADM database of Global Administrative Areas (http://www.gadm.org/). Maps were built using ArcMap 10.2.

https://doi.org/10.1371/journal.pone.0174980.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Species Distribution Model (SDM) evaluation results.

https://doi.org/10.1371/journal.pone.0174980.t002

Risk analysis

Maps representing the results from this risk analysis are presented in Fig 3. For each subtype, the estimated relative risk of human infection, for each cell (approximately 1km²) is represented on a continuous scale from 0 (least risk) to 1 (highest risk) indicated by level of grey shading. H5N1 high risk areas were more diffusely spread around the south-eastern quarter of China, primarily in areas surrounding the Yangtze River delta. H5N1 high risk areas extended towards the Sichuan Basin, Jianghan Plain, river regions of Tibet, and Qinghai lake, as well as plains in Inner Mongolia, Shaanxi, Hebei and Henan. H7N9 high risk areas were more concentrated to the south-eastern coast line starting from Shandong and extended south towards Guangxi, encompassing nearly all of Jiangsu and Shanghai, and plains in Zhejiang, Anhui, Jiangxi, Hunan and Hubei, and Guangdong provinces.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Risk models.

The first panel shows H5N1 and the second panel shows H7N9. Final relative risk values for each cell (approximately 1km²) are represented on a continuous scale of low (light grey) to high (dark grey). Data sources used to develop risk models include: species distribution models (SDMs) 3 and 4 which were produced in this study, domestic chicken population data obtained from Livestock Geo-Wiki (http://www.livestock.geo-wiki.org/), human population data from the LandScan (2014)^™ High Resolution global Population Data Set (http://web.ornl.gov/sci/landscan/). Base maps were obtained from the GADM database of Global Administrative Areas (http://www.gadm.org/). Maps were built using ArcMap 10.2.

https://doi.org/10.1371/journal.pone.0174980.g003

For each subtype, to validate the performance of our risk model, we extracted normalised relative risk values R_kat exact case locations and observed where these points fell on a risk scale of low (0.0 > R_k ≥ 0.25), low-medium (0.25 > R_k ≥ 0.50), medium-high (0.5 > R_k ≥ 0.75) and high (0.75 > R_k ≥ 1.0). In Fig 4, for each subtype, we display the number of exact cases which fell into each of the four categories. The left chart shows categorisation of H5N1 exact case locations corresponding to the H5N1 risk model. The right chart shows categorisation of H7N9 exact case locations corresponding to the H7N9 risk model. For exact H5N1 cases, 7.7% (4/52) of points fell in the low risk category, 30.8% (16/52) of points fell in the low-medium risk category, 36.5% (19/52) of points fell in the medium-high risk category and 25.0% (13/52) of points fell in the high risk category. For exact H7N9 cases, 8.82% (6/68) of points fell in the low risk category, 13.2% (9/68) of points fell in the low-medium risk category, 44.1% (30/68) of points fell in the medium-high risk category and 33.8% (23/68) of points fell in the high risk category. Based on these results, the H7N9 risk model fitted more to our data set than H5N1. Human population data from LandScan (2014)^™ for China had the smallest raster extent, hence our final risk analysis models were clipped to match this extent. One exact H7N9 case did not fall within this extent and hence was excluded from validation tests.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Risk model validation.

Charts show the number of exact case locations which fell into a low risk area (0.0 > R_k ≥ 0.25), low-medium risk area (0.25 > R_k ≥ 0.50), medium-high risk area (0.5 > R_k ≥ 0.75) and high risk area (0.75 > R_k ≥ 1.0). The left chart shows categorisation of H5N1 exact case locations corresponding to the H5N1 risk model. The right chart shows categorisation of H7N9 exact case locations corresponding to the H7N9 risk model.

https://doi.org/10.1371/journal.pone.0174980.g004

We additionally validated our risk model using unexact locations of H5N1 and H7N9 cases, however as we don't know the exact location for the unexact cases, we recorded the maximum R_k within a 5km radius of each unexact case (see S8 Fig). For H5N1, 20.9% (45/215) of points fell in the low risk category, 13.4% (29/215) of points fell in the low-medium risk category, 37.2% (80/215) of points fell in the medium-high risk category and 28.3% (61/215) of points fell in the high risk category. For H7N9 15.3% (187/1221) of points fell in the low risk category, 21.5% (262/1221) of points fell in the low-medium risk category, 39.7% (485/1221) of points fell in the medium-high risk category and 23.5% (287/1221) of points fell in the high risk category.

Final relative risk values, R_k’were aggregated across all secondary administrative areas (prefectures, municipalities, cities) to quantify the relative risk posed to a particular area. The mean value per area was taken as the aggregated relative risk value. High risk areas were considered as those with a mean > 0.5 and are highlighted in Fig 5. For H5N1, these areas included: Foshan, Zhongshan, Dongguan, Zhanjiang (Guangdong); Beihai (Guangxi); Ezhou, Wuhan (Hubei); Nanchang (Jiangxi); Nanjing (Jiangsu); Hefei (Anhui); Xiangtan (Hunan); and Haikou (Hainan). For H7N9, these areas included: Shanghai municipality, Jiaxing, Zhoushan (Zhejiang); Foshan, Zhongshan, Shantou, Dongguan (Guangdong); Changzhou, Nantong, Suzhou, Wuxi, Taizhou, Zenjiang, Yangzhou, Nanjing (Jiangsu), Ma’anshan, Wuhu, Tongling, Chaohu (Anhui); Beihai (Guangxi); Ezhou (Hubei); and Nanchang (Jiangxi).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. High risk areas.

The first panel shows H5N1 and the second panel shows H7N9. Final relative risk values for each cell (approximately 1km²) are represented on a continuous scale of low (light grey) to high (dark grey). High risk areas (1) represent secondary administrative areas (prefectures, municipalities, cities) areas with a mean-aggregated relative risk value (> 0.5), and high risk areas (2) are those areas which have not yet reported a case. Data sources used to develop risk models include: species distribution models (SDMs) 3 and 4 which were produced in this study, domestic chicken population data obtained from Livestock Geo-Wiki (http://www.livestock.geo-wiki.org/), human population data from the LandScan (2014)^™ High Resolution global Population Data Set (http://web.ornl.gov/sci/landscan/). Base maps of Chinese administrative regions were obtained from the GADM database of Global Administrative Areas (http://www.gadm.org/). Maps were built using ArcMap 10.2.

https://doi.org/10.1371/journal.pone.0174980.g005

We also highlighted high risk areas which have not yet reported a H5N1 or H7N9 case (High risk areas (1) in Fig 5). An area is considered to have a H5N1 or H7N9 case if an exact or unexact case fell within the area. For H5N1 these included: Dongguang, Foshan and Zhongshan (Guangdong); Beihai (Guangxi); and Haikou (Hainan). For H7N9 these included: Beihai (Guangxi) and Zhoushan (Zhejiang). The computed risk map can also be downloaded from S3 and S4 Files in ASCII format.

In previous SDMs, LBM density showed strong associations with H7N9 [27,28,30] and H5N1 [30]. LBMs are locations where AIV transmission can be especially amplified due to the mixing of chickens from different farms, overcrowding, and stressful (and thus immunosuppressive) poultry housing conditions. Inclusion of this variable into the risk analysis would account for the high transmission intensity seen in LBM environments.

We were able to collect LBM data (from [27,30], previously purchased from a commercial mapping provider www.autonavi.com), however data collection methods were not available for scrutiny. Due to uncertainties regarding reliability and validity of this data set, we did not incorporate a LBM variable into the formal risk model. Instead, we presented our risk models overlayed with LBM data aggregated to secondary administrative areas (see Fig 6). LBM data was only available for 43 (of 344) prefectures, municipalities and cities within 10 (of 31) provinces in China. Only 10 (of 52) H5N1 exact points and 45 (of 69) H7N9 exact points fell into these regions. Secondary administrative areas with highest LBM counts included Jiaxing (Zhejiang), Beijing, Guangzhou (Guangdong) and Shanghai.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Risk models overlayed with live bird market density.

The first panel shows H5N1 and the second panel shows H7N9. Final relative risk values for each cell (approximately 1km²) are represented on a continuous scale of low (light grey) to high (dark grey). Different sized circles represent the live bird market density per secondary administrative areas (prefectures, municipalities, cities), however data were only available for 43 (of 344) of these areas. Data sources used to develop risk models include: species distribution models (SDMs) 3 and 4 which were produced in this study, domestic chicken population data obtained from Livestock Geo-Wiki (http://www.livestock.geo-wiki.org/), human population data from the LandScan (2014)^™ High Resolution global Population Data Set (http://web.ornl.gov/sci/landscan/). Live bird market data were requested from authors of previously published SDMs [27,30]. Base maps of Chinese administrative regions were obtained from the GADM database of Global Administrative Areas (http://www.gadm.org/). Maps were built using ArcMap 10.2.

https://doi.org/10.1371/journal.pone.0174980.g006