Measuring unequal distribution of pandemic severity across census years, variants of concern and interventions

Abstract

Background: The COVID-19 pandemic stressed public health systems worldwide due to emergence of several highly transmissible variants of concern. Diverse and complex intervention policies deployed over the last years have shown varied effectiveness in controlling the pandemic. However, a systematic analysis and modelling of the combined effects of different viral lineages and complex intervention policies remains a challenge due to the lack of suitable measures of pandemic inequality and nonlinear effects.

Methods: Using large-scale agent-based modelling and a high-resolution computational simulation matching census-based demographics of Australia, we carried out a systematic comparative analysis of several COVID-19 pandemic scenarios. The scenarios covered two most recent Australian census years (2016 and 2021), three variants of concern (ancestral, Delta and Omicron), and five representative intervention policies. We introduced pandemic Lorenz curves measuring an unequal distribution of the pandemic severity across local areas. We also quantified pandemic biomodality, distinguishing between urban and regional waves, and measured bifurcations in the effectiveness of interventions.

Results: We quantified nonlinear effects of population heterogeneity on the pandemic severity, highlighting that (i) the population growth amplifies pandemic peaks, (ii) the changes in population size amplify the peak incidence more than the changes in density, and (iii) the pandemic severity is distributed unequally across local areas. We also examined and delineated the effects of urbanisation on the incidence bimodality, distinguishing between urban and regional pandemic waves. Finally, we quantified and examined the impact of school closures, complemented by partial interventions, and identified the conditions when inclusion of school closures may decisively control the transmission.

Conclusions: Public health response to long-lasting pandemics must be frequently reviewed and adapted to demographic changes. To control recurrent waves, mass-vaccination rollouts need to be complemented by partial NPIs. Healthcare and vaccination resources need to be prioritised towards the localities and regions with high population growth and/or high density.

Keywords: Agent-based modelling; COVID-19; Pandemic inequality; SARS-CoV-2; Urbanisation effects.

Fig. 1

Model of the natural history of three COVID-19 variants: ancestral (blue), Delta (green), and Omicron (red). The illustrated profiles are sampled from 2 random agents. Each profile rises exponentially until reaching the infectivity peak, followed by a linear decrease until full recovery. Vertical lines mark the mean incubation period for the three considered variants (ancestral: blue, Delta: green, and Omicron: red), with the means following a log normal distribution. The mean incubation period and recovery period for each of the variants are reported in Appendix Table 6. The inset shows $R_0$ of the three considered variants

Fig. 2

Five simulated intervention policy scenarios. PRE-VAC: preemptive vaccination prior to the pandemic. NPIs: non-pharmaceutical interventions. Policies are considered to be more stringent moving from left to right. The macro- and micro-parameters for NPI-related policies are summarised in Appendix Table 7. Parameters relating to the vaccination coverage and vaccine efficacy are summarised in Appendix Table 8

Fig. 3

Pandemic Lorenz curves measuring inequality in distribution of the pandemic severity. The black line represents the line of equality where each SA2 contributes equally to the cumulative incidence. A curve closer to the line of equality (i.e., Lorenz Curve A, shown in red) indicates that the contributions of SA2 residential areas towards the aggregate cumulative incidence in response to a specific intervention policy A are more equally distributed than the contributions of these areas under policy B which are traced by the curve shaped further away from the line of equality (i.e., Lorenz Curve B, shown in blue)

Fig. 4

Impact of different intervention policies on pandemic severity for three considered variants simulated for two census years (top row: 2016; middle row: 2021; bottom row: relative change between years). Each column compares the impact of five intervention policies for one variant of concern: a ancestral; b Delta; c Omicron. See Fig. 2 for a detailed description of the considered intervention policies. Coloured shaded areas around solid lines show standard deviation. Each profile corresponds to one intervention policy and is computed as the average over 100 runs

Fig. 5

A comparison of pandemic severity for different policies across three considered variants (ancestral: blue; Delta: green; Omicron: red) and two census years (solid line: 2021; dashed line: 2016). The severity of each variant is measured by cases per million (top row). The change in incidence (bottom row) is calculated as the difference of incidence cases per million between two census years. Each column compares the impact of three variants for one intervention policy: a Policy 1; b Policy 4; and c Policy 5. See Fig. 2 for a detailed description of the considered intervention policies. Coloured shaded areas around solid lines (in bottom row) show standard deviation. Each profile (solid and dashed lines) corresponds to one intervention policy and is computed as the average over 100 runs

Fig. 6

Positive correlation between the usual residential population difference and the peak incidence difference between 2016 and 2021 at SA2 resolution for three considered variants: Ancestral (blue), Delta (green), and Omicron (red). Dashed lines represent linear fitting for each of the profiles (see Appendix Table 10 for statistical analysis). Data points corresponding to each SA2 are computed as the average over 100 runs. Total number of overlapping SA2 between 2016 and 2021 census years: 2147. Pearson correlation coefficients: $r_{\mathit{Ancestral}}=0.7717$, $r_{\mathit{Delta}}=0.6447$, $r_{\mathit{Omicron}}=0.9002$

Fig. 7

Pandemic Lorenz curves measuring distribution of pandemic effects across SA2 areas for considered variants, years and policies. Each column compares the impact of five intervention policies for one variant: a ancestral; b Delta; c Omicron. Top raw: 2016; bottom raw: 2021. Refer to Fig. 2 for a detailed description of the considered intervention policies. Each profile corresponds to one intervention policy and is computed as the average over 100 runs

Fig. 8

Comparison of pandemic waves in Greater Capital Cities (GCCs) and all other areas. Each column compares the baseline scenario (Policy 1) for 2016 and 2021 census data, for a variant of concern: a ancestral; b Delta; c Omicron. Each profile is computed as the average over 100 runs

Fig. 9

Effects of school closures combined with NPIs for three considered variants (log scale): a ancestral; b Delta; c Omicron. School closures effectively control the spread of Delta variant, producing a sharp difference in incidence. Such a bifurcation is not observed for the ancestral and Omicron variants. Each profile corresponds to one intervention policy and is computed as the average over 100 runs. Appendix Fig. 31 shows these plots on linear scale

Fig. 10

Snapshot of the Australian population: locations and geographical areas. Capital city of each state is annotated with their usual residential population (P) and geographical area (A) in 2016 (in red) and in 2021 (in blue), respectively. The state names are abbreviated as follows: Western Australia (WA), Northern Territory (NT), South Australia (SA), Queensland (QLD), New South Wales (NSW), Victoria (VIC) and Tasmania (TAS). Note that changes in the size of geographical areas are due to boundary changes

Fig. 11

Comparison of residential population density and working population density in different states and territories (ST/T) of Australia partitioned by GCCs (a), and other areas (b). Note that the population density scale in (b) is more than 50 times lower than that in (a). The Australian Capital Territory (ACT) predominantly comprises urban areas and thus does not show a comparable density for non-GCC areas

Fig. 12

High variation of the population density across three Greater Capital Cities of Australia: Greater Sydney (left), Greater Melbourne (middle), and Greater Hobart (right). Population-dense areas are highlighted in darker colours, with intensity defined by the colour bar on the right. Note that boundaries in these figures are shown at SA2 resolution

Fig. 13

Geographic representation of Greater Capital Cities (in red), international airports (having international air traffic in 2019, in green), and other regions (in blue) across Australia. The boundaries are delineated based on the 2021 SA2-resolution map of Australia. For visual clarity, some Australian islands are not shown. Hobart is not annotated as there were no incoming international passengers to Hobart during the considered period (i.e., in 2019). The 3-letter code names for the annotated international airports are described in Table 4

Fig. 14

Visualisation of the discrepancies at different resolutions observed in 2021 census commuting data. Each subplot traces combinations between UR and POW for all possible resolutions. Across the UR $\leftrightarrow$ POW flows, considered at comparable or higher resolutions (i.e., SA3–DZN, SA2–DZN, and SA1–DZN), the edge numbers are noticeably lower towards the right hand side on x-axis, with a reduced number of commuters (y-axis)

Algorithm 1

Optimised algorithm generating the 2021 travel-to-work (TTW) surrogate network

Fig. 15

Visualisation of commuting patterns between usual residences (URs) and places of work (POWs), using 2021 census data. Each commuting flow is represented by an edge connecting the centroids of the corresponding UR and POW. Red lines represent the commuting edges between URs (in SA1 level) and POWs (in DZN level) directly taken from 2021 census, while blue lines depict the reconstructed edges created by our algorithm to match the aggregate totals

Fig. 16

Results of network reconstruction using the 2021 TTW data. a Number of commuters between the reference TTW SA2–SA2 network (from ABS Data), the TTW SA1–DZN network (from ABS Data), the TTW SA1–DZN network (refined by the original algorithm of Fair et al. [37], adapted for census 2021), and the TTW SA1–DZN network (refined by the optimised algorithm). b Number of unassigned commuters, counted in SA2–SA2 network but not in SA1–DZN network, over refining iterations. c Mean squared error between the SA2–SA2 TTW network from ABS and the reconstructed SA2–SA2 TTW network from the refined SA1-DZN networks. d Pearson’s correlation coefficient between the SA2–SA2 TTW network from ABS and the reconstructed SA2–SA2 TTW network from the refined SA1–DZN networks

Fig. 17

Weighted linear regression models for the considered airports across Australia. The incoming international passengers (IIP) numbers in previous years up to 2019 are shown in blue; the projected IIP numbers in 2021 are shown in red. Solid red line represents the fitted linear regression model

Fig. 18

Constructing pandemic Lorenz curve for the simplified example, using the data presented in Table 9

Fig. 19

Number of SA2 areas exhibiting an incidence peak in the simulated time period, across three considered policies: Policy 1, Policy 4, and Policy 5, and three variants: a ancestral; b Delta; c Omicron. Each plot is averaged over 100 runs

Fig. 20

Correlation between the usual residential population density difference and the peak incidence difference, computed between 2016 and 2021 at the SA2 resolution for the three considered variants: ancestral (blue), Delta (green), and Omicron (red). Data points corresponding to each SA2 area are derived as the averages over 100 runs. Dashed lines represent linear fit for each of the profiles. Total number of overlapping SA2 areas between 2016 and 2021 census years: 2,147. Pearson correlation coefficients: $r_{\mathit{ancestral}}=0.4150$, $r_{\mathit{Delta}}=0.3426$, $r_{\mathit{Omicron}}=0.4509$. These are lower correlation coefficients than the ones obtained for the population size difference, as reported in main manuscript, Fig. 6

Fig. 21

Correlation between the average household size and the peak incidence between 2016 and 2021 at SA2 resolution for three considered variants: a ancestral, b Delta, and c Omicron. Data points corresponding to each SA2 are computed as the average over 100 runs. Total number of overlapping SA2 areas between 2016 and 2021 census years: 2147. Refer to table 11 for Pearson correlation coefficients

Fig. 22

There is no clear correlation between the difference in average household size and the peak incidence difference, computed between 2016 and 2021 at SA2 resolution for three considered variants: ancestral (blue), Delta (green), and Omicron (red)

Fig. 23

Strong positive correlation between the sub-population residing in large households in an SA2 area and the peak incidence, determined at SA2 resolution across two census years (2016 and 2021) and three considered variants: a ancestral, b Delta, and c Omicron. Data points corresponding to each SA2 are computed as the average over 100 runs. Pearson correlation coefficients r are summarised in Table 12

Fig. 24

Strong correlation between the difference in sub-population residing in large households and the peak incidence difference between 2016 and 2021 at SA2 resolution for three considered variants: ancestral (blue), Delta (green), and Omicron (red). Dashed lines represent linear regression for each of the profiles. Data points corresponding to each SA2 are computed as the average over 100 runs. Pearson correlation coefficients r are summarised in Table 12

Fig. 25

Strong positive correlation between the sub-population residing in small households in an SA2 area and the peak incidence, determined at SA2 resolution across two census years (2016 and 2021) and three considered variants: a ancestral, b Delta, and c Omicron. Data points corresponding to each SA2 are computed as the average over 100 runs. Pearson correlation coefficients r are summarised in Table 12

Fig. 26

Strong correlation between the difference in sub-population residing in small households and the peak incidence difference between 2016 and 2021 at SA2 resolution for three considered variants: ancestral (blue), Delta (green), and Omicron (red). Dashed lines represent linear regression for each of the profiles. Data points corresponding to each SA2 are computed as the average over 100 runs. Pearson correlation coefficients r are summarised in Table 12

Fig. 27

Pandemic Lorenz curves for the considered policies across the three variants and two census years. Each row compares the impact of a policy using 2016 census (left) and 2021 census (right) for three variants. Refer to Figure 2 in the main manuscript for a detailed description of the considered intervention policies. Each profile corresponds to one intervention policy and is computed as the average over 100 runs

Fig. 28

GCCs versus other non-urban areas. Effects of urbanisation on pandemic dynamics for the considered policies, across three variants, using 2016 and 2021 census years. We partition SA2 areas in Australia as Greater Capital Cities (GCCs), and other non-urban areas. Each column compares the impact of five intervention policies for a variant of concern: a ancestral; b Delta; c Omicron. Each plot is averaged over 100 runs

Fig. 29

GCCs and APs versus other non-urban areas. Effects of urbanisation on pandemic dynamics for considered policies, across three variants using 2016 and 2021 census years. We partition SA2 areas in Australia as Greater Capital Cities (GCCs) and those close to international airports (APs), and other non-urban areas. Each column compares the impact of five intervention policies for a variant of concern: a ancestral; b Delta; c Omicron. Each plot is averaged over 100 runs

Fig. 30

Pandemic Lorenz curves for considered variants across policies and census years. Each column compares the impact of five intervention policies for a variant of concern for GCCs and non-GCCs: a ancestral; b Delta; c Omicron. Each profile corresponds to one intervention policy for Australia (red), GCCs (green), or other non-urban areas (blue), and is computed as the average over 100 runs

Fig. 31

Effects of school closures combined with NPIs for three considered variants (on linear scale). a ancestral; b Delta; c Omicron. Figure 9 in main manuscript shows these plots on log scale. Each profile corresponds to one intervention policy and is computed as the average over 100 runs