Evaluation of the effect of chickenpox vaccination on shingles epidemiology using agent-based modeling

Model structure and agent-characteristics

Statecharts related to within-host dynamics of the model are shown in Fig. 1A. A disease statechart where agents are in protected, susceptible or disease states determined an agent’s probability of contracting chickenpox or shingles. The chickenpox and shingles vaccination schedules’ statecharts represented which vaccination each agent received. The chickenpox vaccination schedule was modeled based on the Alberta VZV vaccination policy (Government of Alberta, 2015). Agents become due for chickenpox vaccines at the ages of 12 months for the first dose and at 4–6 years—for the second dose (Fig. 1B). A representation of one-dose shingles vaccination given to individuals 50 years or older is depicted in the statechart shown in Fig. 1C; however, currently shingles is not part of the publicly funded schedule in Alberta and therefore was not included as part of this analysis. An agent receives a vaccine dose with a probability based on its vaccine attitude as described below. If an agent receives the second dose but has not yet received the first, it may also receive a catch-up of the first dose with a certain probability. These probabilities are specified by parameters, as shown in Table 1. A demographic statechart represented Alberta births, ageing and mortality characteristics (Fig. 1D).

An agent’s chance of being infected with chickenpox was dependent on whether they came into contact with someone with infection and the risk of transmission on such a contact. In comparison, an agent’s likelihood of getting shingles was determined by their individual waning immunity timer (i.e., after chickenpox infection, a countdown on “Immunity Waning Time” initiates; when it completes, the agents become susceptible to shingles) as well as the number and degree of boosting they received. The immunity waning time was represented by a formula derived from Ogunjimi et al. (2015) and further by calibration (Fig. S1; Table 2). This equation inherits assumptions whereby force of reactivation is represented by a gamma distribution, initial CMI is represented by a normal distribution and the waning of immunity rate for shingles is a fixed rate that can be altered using a coefficient—waning of immunity coefficient (WoI) (Ogunjimi et al., 2015). Furthermore, to account for the small but sharp increase in childhood shingles cases as observed in Russell et al. (2014), we included a small proportion of the population between the ages of 0 and 19 (5%) who had a short waning of immunity timer, allowing the model to account for individuals who may have a weaker force of boosting, lower initial CMI (e.g., immunocompromised) or a very quick waning of the immunity rate.

An agent’s chance of getting shingles was determined by the number of years they were protected through boosting. The model assumed progressive boosting, as postulated by Guzzetta et al. (2016) and therefore the number of years of boosting protection was calculated by multiplying the number of times an agent comes into significant contact with the VZV by the duration in years of each boosting event. The number of boosting events was determined through the distance-based contact network and the duration of each boost was equal to the number of years of added protection from shingles derived through each boost. Altering the quantitative values of DoB and WoI of shingles resulted in significant changes in the incidence of shingles in a population.

Contacts, network and spatial context

The VZV model represented agents in a stylized geographic area, where agents are connected to other agents based on their proximity to one another (a distance-based network) (Fig. S2). When agents are randomly placed in the model environment they are connected to all other agents within their connection range (connectionRange_Norm). We included a low-density periphery and a higher density central region to better represent how disease spreads in a typical public health district in Alberta (or comparable jurisdiction) spanning an urban center and rural regions. Approximately 20% of our population were part of low density regions and the remainder, 80% resided in a high-density region.

The populationDensityUrban and the populationDensityRural determine the number of agents in a respective area per unit of distance in the model and along with the connectionRange determine the number of connections between agents. The contact rate (baseContactRate_Norm) determines the number of contacts (e.g., messages) an agent makes to a connected agent per day. It is through these messages, sent through connections and contacts, that individuals can “transmit” chickenpox infection and provide boosting of VZV. To represent the increased contact time and range of connectivity among day-care and early school-aged children, we incorporated different connection ranges and contact rates based on age (Mossong et al., 2008). Those aged one to nine years were considered “preferential contacts” in our model, and therefore have a higher contact rate (baseContactRate_Pref) and connection range (connectionRange_Pref) when interacting with individuals within the same age range. We chose this group to represent “Preferential Contacts” as these are the age groups with the highest rates of chickenpox infection as described by Kwong et al. (2008). These network adjustments ensured more realistic contact network assumptions than those of random-mixing and compartmental models, such that not only age-group preferences of contacts were captured, but also increased global connectivity due to bridging effects of younger age groups was considered, both in a spatial context. To characterize the transmissibility of VZV, we included parameters whereby a susceptible agent will have a certain probability of infection per contact message received from a normal chickenpox case (probCPDiseaseOnContact), a breakthrough case of chickenpox (probCPDiseaseOnContactWithBreakthrough) and a shingles infection (probCPDiseaseOnContactWithShingles). The probability of infection in our model is dependent on whether a susceptible agent would come in contact with and receive a message from an infected case of chickenpox and shingles. For example, an infected agent in the preferred age group (i.e., a child) may send 20 messages per day, each to a random connected agent. Considering that the mean number of connections per agent in this age group is 10.87, the number of messages per connection-day equals to 20/10.87 = 1.8 ≈ 1–2. Outside of the preferred age group, this equates to 3/7.42 = 0.40 < 1 message per connection-day (Table 1; https://figshare.com/articles/Chickenpox_and_shingles_ABM/5294647/1). Some of these messages do not lead to a receiving agent becoming infected because the receiving agent is not susceptible, or an agent is not selected by a model to be connected to an agent sending a message (not being part of contact pool). A complete age-specific connection and contact frequency profile is available at https://figshare.com/articles/Chickenpox_and_shingles_ABM_connection_and_contact_frequency_profile/5532763/1.

Parameterization

To parameterize the model, we conducted a comprehensive literature review of chickenpox and shingles disease, including modeling, epidemiological, and immunology studies and drew evidence from Alberta’s Interactive Health Data Application (IHDA) (Alberta Health, 2017). The main parameters are listed in Table 1. We drew demographic data, including Alberta life tables, and population and age distributions, and Canadian fertility rates from Statistics Canada (2008, 2016a, 2016b). Chickenpox vaccination parameters, such as those associated with primary vaccine failure and waning of vaccination immunity were derived from literature (Gershon, Takahashi & Seward, 2012; Bonanni et al., 2013; Duncan et al., 2017). We built the mechanism whereby vaccine coverage was generated by our model based on distribution of vaccination probabilities and population vaccination attitudes as described by Doroshenko, Qian & Osgood (2016). We classified all individuals into three groups: those who accept, reject and are hesitant to receive vaccination, and we assigned vaccination probabilities for each of these groups. We used calibration to ensure that model-generated vaccine coverage rates were comparable to those reported in Alberta (Alberta Health, 2017). For the baseline scenario in the main experiment, chickenpox vaccination maintained an average coverage for the first dose of 85.58% (95% CI [85.54–85.62]) and second dose coverage of 80.28% (95% CI [80.24–80.32]) across all years and all model runs.

Calibration and validation of the model

We calibrated our model using a step-wise approach.

First, contact and mixing patterns, as well as the probability of infection per contact message, were determined using a pattern-oriented modeling approach. Pattern-oriented modeling involves identifying multiple patterns of behavior in a complex system and constructing model structures and parameters that replicate those patterns (Grimm, 2005). The rationale for using this approach is that transmission of infection in our model (as in real life) is driven by combination of the location of agents, pathogen transmissibility per discordant contact and intrinsic characteristics of agents reflecting their susceptibility (or lack of) to infection. Therefore, to represent this multi-faceted transmission in our model, we varied three model parameters by hand (connection range, contact rate, and probability of infection on contact message), until the model outputs approximated empirical data, specifically age-specific incidence rates for chickenpox (Kwong et al., 2008). Calibration for this variable was completed once our model generated patterns of disease very similar to what was observed in empirical data (Fig. 2). The best values for probability of infection on contact message, contact rates and connection ranges derived from this calibration process are shown in Table 1. The probability of transmitting infection from a chickenpox breakthrough case and from shingles case were assumed to be 0.3 of the probability of transmitting infection from a normal chickenpox case (Gershon, Takahashi & Seward, 2012). As we only had empirical data for medically-attended chickenpox (Kwong et al., 2008), we tested different underreporting factors and found that reported chickenpox cases constituted approximately 40% of all cases in our model (an underreporting factor of 2.5, consistent with previous studies that suggest the degree of underreporting can range between 2.5 and 7.7 in a country where chickenpox is notifiable) (Ciofi Degli Atti et al., 2002). To further validate our model, we tested model outputs against secondary attack rates for chickenpox reported in the literature. Our model generated a secondary attack rate of 41% in comparison to attack rates observed in the household settings, ranging mostly between 60% and 78% among unvaccinated contacts (but reported as low as 31% among individuals > 15 years old) (Gershon, Takahashi & Seward, 2012; Ceyhan, Tezer & Yildirim, 2009; Seward et al., 2004). We considered this a realistic estimate as contacts in our population intuitively should be less intensive than the household contacts represented in the literature. Furthermore, the overall chickenpox incidence and the age-specific proportion of the population with varicella antibody when vaccination was disabled in the model was consistent with Canadian data prior to vaccination (Fig. S3) (Public Health Agency of Canada, 2012; Ratnam, 2000).

Second, we ran an automated AnyLogic calibration experiment to determine the values for four other unknown parameters in our model, DoB, WoI (which was used to calculate the waning of immunity rate), exogenous infection rate (i.e., the number of chickenpox cases brought in from outside the simulated population) and shingles connection range (i.e., modification factor for the connection range for people with shingles, to account for the closer contact required to spread VZV through shingles in comparison to chickenpox). The calibration experiment automatically varied these parameters over continuous ranges in order to minimize the objective function quantifying the difference (square root of the average of square differences) between the model and actual data (i.e., age-specific incidence of shingles as described in Russell et al. (2014)) (Fig. 2). A full description of the data used in model fitting is shown in Table S1. This automated calibration experiment predicted an exogenous infection rate of 17.8 per year and a shingles connection range modifier of 0.124. These values are shown in Table 1. Furthermore, based on the automated experiment values for DoB of 0.42 years, and a WoI of 0.45/year, we produced age-specific incidence rates for shingles and chickenpox consistent with empirical data.

Third, using the automated calibration values described above as a starting point, we then further investigated whether different combinations of DoB and WoI could also reproduce empirical data. We compared the incident rates resulting from a range of plausible values for DoB and WoI to the empirical incidence rates for chickenpox and shingles (Russell et al., 2014; Kwong et al., 2008; Topping, Høye & Olesen, 2010; Railsback & Grimm, 2012). At this stage, we considered paired values of DoB and WoI to have satisfied calibration if the resulting overall rate of chickenpox and shingles, over a model run of 100 years on 50,000 population, was within 10% of the empirical values (Topping, Høye & Olesen, 2010; Railsback & Grimm, 2012). Based on these guidelines, we found nine other values for a combination of DoB and WoI that met this step of calibration.

Finally, we then ran the 10 combinations that met the criteria noted above for a period of 100 years on a 500,000-population sample (10-fold larger than in the previous step). Output from each run was then compared to empirical data on age-specific incidence rates for chickenpox and shingles. We calculated the absolute median difference between model incidence rates of shingles and empirical data for each age group to determine which combination had the smallest difference as determined by lowest sum of residuals squared. The combination with the lowest sum of residuals squared (27.95) for shingles was then considered our baseline scenario in our main experiment (which represented best fit) (https://figshare.com/articles/Sum_of_residuals_squared_-_VZV_model/5466556/1). The sum of residuals squared from each combination was then compared to the baseline scenario and tested for statistically significant difference using the Mann–Whitney U-test at the 5% level of significance. Based on the findings of the statistical tests and as a form of scenario analysis, we included both the baseline scenario along with any other DoB and WoI combinations that were statistically not different (p-value > 0.05) as scenarios in our main experiment.

Main experiment

For each of the six scenarios chosen for the main experiment, we conducted at least 30 paired runs, with and without chickenpox vaccination. Running paired runs for each scenario allowed the measurement of differences in shingles rates with and without vaccination for each DoB and WoI combination. Paired model runs were given the same random seed to ensure consistency between pairs. Chickenpox vaccination was represented as part of a two-dose schedule, as described above. In our model, chickenpox vaccination started at 25 years after the initialization of the model and continued for 75 years. The count and cumulative incidence rate of shingles cases were compared at the 5% level of significance using the Mann–Whitney U-statistical test between the runs with and without chickenpox vaccination at four different time periods, specifically at 10, 25, 50, and 75 years, following introduction of vaccination.

Sensitivity analysis

We ran several sensitivity analyses that tested how varying vaccination parameters may impact the count of shingles cases following chickenpox vaccination. Specifically, we compared a one-dose to a two-dose chickenpox vaccine schedule. In a separate set of analyses, we compared coverage rates by moving 10% of hesitant individuals to vaccine acceptors (higher vaccine coverage) in one sensitivity analysis and 10% of hesitant individuals to vaccine rejecters (lower vaccine coverage) in the second sensitivity analysis. Furthermore, we tested the impact of removing the boosting of shingles immunity (i.e., positing no added years of protection on contact with a chickenpox/shingles case) to see the overall impact of removing this biological effect from the shingles incidence estimates both before and after chickenpox vaccination.

Input calibration

Based on the last step of our calibration we found that the age-specific shingles residuals for DoB five years were not significantly different to five of the other combinations (Table 2). These combinations were then included in our main experiment as six scenarios (Scenario 2–7, sum of residuals squared for shingles ranged from 27.95 to 41.51, p-value > 0.05) (https://figshare.com/articles/Sum_of_residuals_squared_-_VZV_model/5466556/1). All sum of residuals squared for chickenpox by age group were not statistically different for each calibration scenarios (sum of residuals squared for chickenpox ranged from 249.20 to 251.00, p-value > 0.05). As one of the calibration validity tests, we tried to fit the model to empirical data with no boosting, however when we disabled boosting we could not replicate the age-specific incidence rate observed in Alberta prior to vaccination (Russell et al., 2014).

Main experiment

Chickenpox vaccination led to a large drop in chickenpox cases across all six scenarios. In the baseline scenario, the cumulative incidence of chickenpox had dropped from 1,254 cases per 100,000 person-years pre-chickenpox vaccination to 193 cases per 100,000 person-years 10 years after the vaccine implementation. The cumulative incidence of chickenpox was further reduced to 49 cases per 100,000 person-years 75 years following vaccination. In comparison, all scenarios from 10-years and 25-years post-vaccination showed significantly greater shingles incidence with vaccination compared to the no-vaccination (Table 3). However, the degree of this increase and its subsequent decline was markedly different between experiments (Fig. 3; Table 3). For instance, in Scenario 2 there was an increase of approximately 22.47 cases per 100,000 person-years after 10 years; by comparison, in Scenario 7 the magnitude of the increase was greater at 99.71 cases per 100,000 person-years over that decade (p < 0.001). At 75 years post-vaccination, cumulative incidence ranged from a decline of 69.6 to an increase of 71.3 per 100,000 person-years for two and seven years of boosting respectively (p < 0.001). The increase in the incidence of shingles was positively associated with the DoB and showed a positive relationship with cumulative shingles incidence rising with each incremental increase in DoB. By the 75-year interval, the shingles incidence in all experiments was below the rate with no-vaccination (Fig. 3); however, the cumulative incidence in Scenarios 5–7 was still higher with vaccination than without. In all our scenarios, vaccination did eventually lower the overall rate of shingles; however, the amount of time was dependent on the DoB and the WoI.

At 10-year interval, all age groups greater than 10 years old showed small but equal increase in the number of shingles cases. However, at each subsequent time point (as a greater number of persons of younger age are protected with chickenpox vaccination from initial infection with the virus), cases of shingles were progressively more concentrated in the older age-groups. In contrast, in the younger age group we observed a decrease in shingles almost immediately following chickenpox vaccination (Fig. 4).

We noted several interesting observations in our model. In the no-vaccination model runs, there was an overall increase in shingles cases even without vaccination; in the baseline scenario, the incidence of shingles was 418 per 100,000 person-years after 10 years and increased to 425 per 100,000 person-year after 75 years. Moreover, the average number of boosts per individual calculated in our model over their life-time (with no chickenpox vaccination) was 1.83 with majority of people receiving 0 or 1 boosts. We found minor differences in chickenpox and shingles incidence between urban and rural communities. The incidence of shingles was slightly higher in rural settings in the absence of vaccination but became lower with vaccination (Fig. S4). Approximately 5% of urban agents had connections with rural agents (Fig. S10).

Sensitivity analysis

None of the vaccination sensitivity analyses had a major impact on the number of shingles cases post-vaccination. As expected, we saw an increase in the number of shingles cases with a higher number of vaccine acceptors and therefore higher coverage than baseline. Conversely, we saw fewer shingles cases post chickenpox vaccination when we only gave one dose of the vaccine and when we had more vaccine rejecters, however, the effect was minimal. The difference between the low coverage (higher number of vaccine rejecters) and the high coverage (lower number of vaccine rejecters) was only 2 shingles cases per 100,000 person-years. By contrast, removing the biological effect of boosting had a substantial impact on the number of shingles cases both before and after vaccination, with cases rising to levels much higher than what was seen in Alberta (Russell et al., 2014). Without chickenpox vaccination, the shingles rate was very high at 734 cases per 100,000 person-years. Simultaneously, this is the one analysis where the rates begin to decline immediately following vaccination with a decline in cumulative incidence of 215 per 100,000 person-years 75 years after the implementation of the chickenpox vaccine (Table 4).