Study area
The methodology scheme is summarised in Additional file 1. The study location was the city of Mashhad, capital of Razavi-Khorasan Province in north-eastern Iran. Mashhad is located between latitudes 36°10’and 36°25’N, and longitudes 59°25’and 59°46’E, covering an area of 307 km2 (Fig. 1). According to the 2016 national census statistics, the city population was almost 3.3 million then and just slightly higher at the time of the study [39]. Mashhad consists of 17 municipality regions, 175 districts and 1301 census tracts (CTs). In this study, we considered the CTs as they provide the finest resolution for accessibility analysis. The average population density of CTs were 20 052 \(\pm\) 10 983 individuals per km2. At the time of conducting this study, there were 26 active public hospitals and 271 public healthcare centres in Mashhad (Fig. 1).
Data
Data on the 26 public hospitals and 271 public healthcare centres, including capacity that depends on available equipment and staff, were obtained from Mashhad University of Medical Sciences. Demographic characteristics of all CTs (including population size stratified by age and area) were obtained as a GIS shapefile from the Statistics Centre of Iran. The addresses of all healthcare centres and hospitals, city boundary and road network were retrieved from Mashhad Municipality and updated based on the Google Maps (https://www.google.com/maps) and OpenStreetMap (https://www.openstreetmap.org/) websites. All above data are freely available for public use via the link provided in the data availability section at the end of the article.
Development of the age-integrated E2SFCA
Calculating the weighted population
To measure the PSA to vaccination centres, the age structure of each CT was integrated into the E2SFCA method as a local influential factor of COVID-19 mortality. However, we weighted each age group in the accessibility formula according to United States Centers for Disease Control and Prevention (CDC) in Atlanta, GA, USA [40]:
$$\begin{gathered} Weighted\,population = \hfill \\ Pop_{0 - 4} \times 1 + Pop_{5 - 17} \times 1 + Pop_{18 - 29} \times 10 + Pop_{30 - 39} \times 45 + Pop_{40 - 49} \times 130 + Pop_{50 - 64} \hfill \\ \times 440 + Pop_{65 - 74} \times 1300 + Pop_{75 - 84} \times 3200 + Pop_{ > = 85} \times 8700 \hfill \\ \end{gathered}$$
(1)
where Popx denotes the population for age group x. The age-related weights were applied, using ArcGIS v.10.8 software (ESRI, Redlands, CA, USA) to the CTs layer.
CA identification
The CA is the basis of the E2SFCA method [25]. According to previous studies [18, 24, 26, 41], we defined these areas by radius as 1 km (CA-1), 1.5 km (CA-2) and 2 km (CA-3) taking into account the average population density of the city (~ 20 000 people per km2). The distances chosen were the routine minimum and maximum values for defining service areas for health facilities in Iran’s major cities [42]. Moreover, the speed limit of 30 km/h (based on the average speed when driving a car in the city) was considered the basis of the travel-time analysis. The buffer analysis tool in ArcGIS software was used to calculate the CAs. The network analysis tools in QGIS v.3.20 open-access software v.3.20 (https://qgis.org/en/site/forusers/download.html) were used to identify the service areas.
Accessibility calculations
In step 1, the CAs were set at 1, 1.5, and 2-km distance to the jth healthcare location. We searched all weighted population locations (k) within the threshold travel-time zone (Dr) from healthcare centre j (CA j) using the travel-time zone (catchment) as the area formed based on the farthest accessibility boundaries (threshold) of each CT centroid to each potential vaccination centre and indicated by the distance (km) or time elapsed. This threshold was plotted based on the closest distance to vaccination centres in kilometres and we computed the vaccination capacity-to-weighted population ratio Rj within the CAs using Eq. 2 below following previous studies [24, 43, 44]:
$${R}_{j}=\frac{{S}_{j}}{{\sum }_{k\in \left\{{d}_{kj}\in {D}_{r}\right\}}{P}_{k}{W}_{r}}=\frac{{S}_{j}}{{\sum }_{k\in \left\{{d}_{kj}\in {D}_{1}\right\}}{P}_{k}{W}_{1}+{\sum }_{k\in \left\{{d}_{kj}\in {D}_{2}\right\}}{P}_{k}{W}_{2}+{\sum }_{k\in \left\{{d}_{kj}\in {D}_{3}\right\}}{P}_{k}{W}_{3}}$$
(2)
where Pk is the population of the kth CT falling within the CA j (dkj \(\in\) Dr); Sj the vaccination capacity at healthcare centre j; dkj the travel time between k and j; and Dr the rth travel time zone (r \(\in \{\mathrm{1,2},3\}\)) within the CA in question. Wr represents the distance weight for the rth travel-time within the CA calculated by the Gaussian function. The weights set (1.00, 0.68, 0.22) were used to capture the distance decay of access to the jth healthcare centre.
In step 2, we searched all locations j for people in ith CT within the 1, 1.5, and 2-km distance radius. Then, using Eq. 3 below, we summed up the vaccination capacity-to weighted population ratios Rj (calculated in step 1) for all CTs. The same distance weights derived from the Gaussian function were applied in different travel-time zones to account for the distance decay.
$$A_{i}^{F} = \sum\nolimits_{{j \in \left\{ {d_{ij} \le D_{r} } \right\}}} {R_{j} } W_{r} = \sum\nolimits_{{j \in \left\{ {d_{ij} \le D_{1} } \right\}}} {R_{j} } W_{1} + \sum\nolimits_{{j \in \left\{ {d_{ij} \le D_{2} } \right\}}} {R_{j} } W_{2} + \sum\nolimits_{{j \in \left\{ {d_{ij} \le D_{3} } \right\}}} {R_{j} } W_{3}$$
(3)
where \({A}_{i}^{F}\) represents the accessibility vaccination centre for the population at location i; Rj the vaccination capacity-to-weighted population ratio at healthcare centre j that falls within the CA of population centre i (dkj \(\in\) Dr); and dij the travel time between i and j. The E2SFCA calculations were performed under three different scenarios as follows:
Scenario 1: Accessibility to public hospitals (PHs)
In many developing countries, including Iran, PHs act as general and special care facilities and the first point of contact when the patient is referred to specialist care [26]. Hospitals are usually well-equipped and can thus be used as public vaccination sites when needed. In this scenario, the PSA to PHs as a potential vaccination centre (PVC) was calculated considering all 26 active PHs (Fig. 1). Since not all employed at the hospital were qualified and available to perform vaccinations, according to hospitals’ official, only 4% of each hospital’s staff capacity (using an average capacity of 24.85 people as vaccinators per hospital) was entered into the accessibility measurement.
Scenario 2: Accessibility to public healthcare centres (PHCs)
The PHCs, also known as primary care centres, are the second main health facilities that can be used for public vaccinations during pandemics. In this scenario, 271 PHCs with a capacity of 1 to 5 people (with an average capacity of 2.05 people as vaccinators), were included in the E2SFCA model. The number of PHCs is almost 10 times higher than that of hospitals and they are well-dispersed across the city. Therefore, the PHCs have a stronger potential to act as vaccination centres.
Scenario 3: Joint accessibility to PHs and PHCs
In this scenario, all 26 PHs and all 271 PHCs were entered into the E2SFCA model. All E2SFCA method calculations to measure the PSA index were performed using ArcGIS v.10.8.
Choosing the most appropriate centres for vaccination (the proposed model)
Many countries have used available public space for vaccination, e.g., cinemas, shopping malls, even outdoor areas such as football arenas and the like. However, we did not take this possibility into account since our primary aim was to evaluate the capability of the modified version of the E2SFCA model to measure the PSA to available medical centres. Since it is not feasible to equip and prepare any resource for public vaccination, especially in developing countries, PHs and PHCs with high scores were selected for administering vaccines based on accessibility index, service capability, distance to main roads and capacity as vaccination centre as follows:
Accessibility index
This index was derived from the calculation of PSA in scenario 3 (Fig. 2A) classified as five categories depending on the PSA for the following classes of CTs: very high (PSA = 2.50E−05 to 3.70E−05), high (PSA = 1.70E−05 to 2.50E−05), medium (PSA = 1.10E−05 to 1.70E−05), low (PSA = 6.00E−06 to 1.10E−05) and very low (PSA = 0 to 6.00E−06). The CTs falling into the very low category were considered as the highest priority for improved accessibility to vaccination. Therefore, higher weights were assigned to healthcare centres in CTs with low and very low PSAs. The weighting was employed to enhance the access to the nearest healthcare centre for people in deprived areas.
Service areas
Equitable accessibility requires optimum allocation of health facilities [45]. However, these areas do not have uniform geometric shapes, so the QGIS version of hexagonal tessellation was applied to achieve more homogeneous and geometrically defined service areas for the healthcare centres (depending on the area and population density of the city). The city was therefore divided into 50 equal 1000-ha hexagons (Fig. 2B). The radius of these theoretical service areas was set at 2 km, i.e. the maximum standard service area (coverage) for access to a health facility [46]. Then, the available PVCs closest to the geometric centres of these hexagons were given higher weights and consequently selected as centres with comparatively high suitability as vaccination sites.
Distance to main roads
Travel time and distance from the location of residents to health or medical facilities were the foundation when calculations were carried out to find locations with high accessibility [47]. In other words, medical facilities located at a convenient distance from main roads decrease travel time to the vaccination services, while larger distances to healthcare facilities (that translates into longer travel times) obstruct visits and revisits. There is no universally accepted superior relation between medical care facilities and the range of roads to reach it. According to Silalahi et al. [48], more attention should be paid to the transport network for easy access to PVCs. In this study, distances of 100 m (very desirable), 150 m (desirable) and 200 m (somewhat desirable) from the main roads were defined as one of the sets of criteria for selecting PVCs when applying buffer analysis in ArcGIS. The PVCs within the distance buffers were then assigned weights, and the PVCs outside these buffer zones with no access to main roads (especially for public transport users) were excluded from the analysis (Fig. 2C).
Vaccination centre capacity
In this study, PHs and then PHCs were ranked according to their capacity, so that centres with higher capacities were given a higher chance of selection. At the same time, the centres without proper facilities and equipment (e.g., centres with two vaccinators or less) were not given priority status (Fig. 2D).
The weighted decision-matrix method proposed by Pugh [49], practical both for simple and complex decisions [50], was used to select the most appropriate centres. First, a 297 × 4 matrix was formed. Then, each criterion was normalized through numbers between 0 and 1. According to the difference in importance of each main criterion, the weights were calculated using the analytical hierarchy process (AHP) [51] based on the opinions of 10 health professionals. Afterwards, the obtained weights from the main indices were multiplied by the corresponding numbers to arrive at the actual values of the selected criteria. For each potential vaccination centre, the obtained number from this calculation gave the final score that was finally used to rank the PVCs. The top 20% of the centres, i.e. those with the highest capability for vaccination according to the final scores, were selected. Business Performance Management, Singapore (BPMS) free web-based AHP online system [52] was used to calculate the weights of the various criteria.
Accessibility to PVCs (reallocated centres)
After selecting the PVCs, the PSAs to these facilities were measured using the E2SFCA methodology.
Evaluation of the model
Global Moran’s Index (GMI) was used to measure the spatial autocorrelation of the accessibility index in different scenarios and the proposed model. We assumed that a decrease in spatial autocorrelation of accessibility within CAs indicated improved accessibility. The GMI is defined by Eq. 4 [53]:
$$I = \frac{{n\left( {\mathop \sum \nolimits_{i = 1}^{n} \mathop \sum \nolimits_{j = 1}^{n} w_{ij} \left( {x_{i} - \overline{x}} \right)\left( {x_{j} - \overline{x}} \right)} \right)}}{{\left( {\mathop \sum \nolimits_{i = 1}^{n} \mathop \sum \nolimits_{j = 1}^{n} w_{ij} } \right)\left( {\mathop \sum \nolimits_{i = 1}^{n} \left( {x_{i} - \overline{x}} \right)^{2} } \right)}}$$
(4)
where n is the total number of spatial divisions (i.e. the CTs); xi the value of the PSA in CT i; \(\overline{x}\) the arithmetic mean for a given PSA; and wij the spatial weight matrix based on inverse distance and the Euclidean distance (distance band = 2983.7 m and the number of nearest neighbours = 4). The value of Moran’s I ranges from -1 to + 1. The further away it is from zero, the stronger (positive or negative) the autocorrelation [54]. ArcGIS v.10.8, software was used to compute the Moran’s index. A P-value < 0.01 was considered significant in the test using 99 permutations.
Visualisation of the selected PVC areas of influence
This was done by Thiessen polygons, generated around a set of points in a given space by assigning all locations in that space to the closest member of the point set, a type of spatial tessellation called Voronoi diagrams [55]. Each polygon was created by QGIS and can be seen as an area of influence of a point in the given set [55]. The CTs were introduced as the origins and the PVCs as the destinations. In Fig. 4, the hub-distance (shown as red lines) indicates the distance (in km) from the centre of each CT (origin) to the nearest PVC (destination). The boundaries (in black) of the Thiessen polygons indicate the coverage of the service area of each PVC (Fig. 4A).
