Assessing the benefits of early pandemic influenza vaccine availability: a case study for Ontario, Canada

New vaccine production technologies can significantly shorten the timelines for availability of a strain-specific vaccine in the event of an influenza pandemic. We sought to evaluate the potential benefits of early vaccination in reducing the clinical attack rate (CAR), taking into account the timing and speed of vaccination roll-out. Various scenarios corresponding to the transmissibility of a pandemic strain and vaccine prioritization strategies were simulated using an agent-based model of disease spread in Ontario, the largest Canadian province. We found that the relative reduction of the CAR reached 60% (90%CI: 44–100%) in a best-case scenario, in which the pandemic strain was moderately transmissible, vaccination started 4 weeks before the first imported case, the vaccine administration rate was 4 times higher than its average for seasonal influenza, and the vaccine efficacy was up to 90%. But the relative reductions in the CAR decreased significantly when the vaccination campaign was delayed or the administration rate reduced. In urban settings with similar characteristics to our population study, early availability and high rates of vaccine administration has the potential to substantially reduce the number of influenza cases. Low rates of vaccine administration or uptake can potentially offset the benefits of early vaccination.

-Comparing vaccination strategies impact on the relative reduction of overall deaths, hospitalizations and clinical attack rates.
Points represent the mean relative reduction of the relevant outcome when compared with the same scenario without vaccination. Each column represents a combination of vaccine efficacy and pandemic transmissibility. Each row represents a different outcome (first row deaths, second row hospitalizations and last row clinical attack rates). Solid lines represent the strategy PVS1 and the dashed line RVS. In all panels, none of the paired curves difference is statistically significant (the credible intervals [not shown] overlap). All panels show results with age-specific contact data calibrated on a Canadian retrospective study [32].  Figure S4 -Relative reduction of the age-specific clinical attack rates.
All panels show the mean relative reduction in age-specific clinical attack rates for vaccination strategies RVS (top half) and PVS1 (bottom half). We observed, as expected, a higher reduction among the under 5 group than the 5-18 one for PVS1 strategy.  Vaccintion time lag (days)

Relative reduction
Vaccine administration rate (per 100,000 per day) 150 300 600 1200 Figure S6 -Relative reduction of overall clinical attack rates with alternative age-specific contact rate parameterization.
Same legends as Figure 1, except that all panels show results under age-specific contact data calibrated on a French face-to face contact [33,34] (as opposed to data sourced from a Canadian retrospective study [32]  Each panel represents, for a given transmissibility and vaccine efficacy scenario, the mean relative reduction of the overall deaths (solid lines) and hospitalizations (dashed line) when compared with a scenario without vaccination. The vaccination time lag is represented on the x-axis. All panels show results under the PVS2 strategy and age-specific contact data calibrated on a Canadian retrospective study [32].
(The difference between Figure S7 and Figure 2 is that the former was run with strategy PVS2 and the latter with PVS1)  Figure S9 -Sensitivity Analysis.
A sensitivity analysis was performed on five model parameters. For each panel, the y-axis represents the range of mean relative reduction change for CAR when the associated parameter value was varied. The vaccination lag (-28 and +14 days) is represented by the x-axis. The bar fill color represents the vaccine administration rate. Each row is associated with one of the five parameters. The left column displays results for a 40% vaccine efficacy, the right column for a 90% efficacy. First row: the infectiousness ratio for asymptomatic cases was varied between 0.05 and 0.4 (baseline value was 0.1). Second row: contact age assortativity parameter was varied between 0.05 and 0.2 (baseline value was 0.1). Third row: contact rate coefficient of variation was varied between 0.3 and 1.5 (baseline value was 0.75). Fourth row: frailty standard deviation parameter was varied between 0.02 and 0.2 (baseline value was 0.1). Fifth row: maximum level of cellular immunity index was varied between 0.1 and 0.9 (baseline value was 0.45).   This supplemental material describes the details and implementation of an agent-based model used for the article "Assessing the benefits of early pandemic influenza vaccination: a case study for Ontario, Canada" published in Scientific Report (DOI: 10.1038/s41598-018-24764-7).

Overview
We used an in-silico population to simulate fairly realistically the spread of pandemic influenza in the province of Ontario, Canada, and evaluate the potential benefits of early vaccine availability. To this end, we developed an agent-based model in order to capture the heterogeneities of transmission dynamics.
The implementation of the agent-based model is not event driven. Instead, the time unit (i.e., day) is divided into relatively coarse "slices" that represent specific time-periods during the course of a day for relevant social activities. Within a given time slice, epidemiological events may occur and are computationally accounted for only once. This approximation allows significant computing time savings for an acceptable loss of realism (which depends on the size of time slices). Although each individual is uniquely identified, the information that who acquires infection from whom is omitted for performance (however, the time of transmission is always recorded to keep track of the generation interval).
The model is implemented with a coarse spatial structure. The fundamental spatial unit is called a "social place", which broadly represents any physical place that individuals may interact. Because it is preferable to have a limited number of types for social places, only the most relevant ones for which data exist are explicitly represented (e.g., households, schools, workplaces, public transportation). We also used a generic social place type that represents all other types of social places. Social places are not geolocated.
The natural history of influenza and transmission processes for each individual are determined probabilistically. Individuals are assigned age-dependent levels of immunity and frailty. These levels determine the likelihood of disease transmission, symptomatic infection, hospitalization and disease-induced death, as detailed in the following sections.
The programming language for the computational model is C ++ . The executable is wrapped as a R [26] library for output analysis. The execution time of this computer program, for a single simulation with 50,000 individuals, is about 30 seconds on a single 2.1 GHz AMD Opteron processor.

Social place
The model is implemented with a spatial structure, and its fundamental spatial unit is called a "social place". It represents real world locations where individuals may have relevant contacts for disease transmission (e.g., households, schools, workplace, public transportation). Because all types of such social places cannot be explicitly described in the model, a social place typed "other" is used for all other real-world locations. In order to have a similar hierarchical structure as real world data, social places are gathered in "area units", which constitute "regions".
The size of a social place, that is the number of linked individuals, is pre-specified with a distribution given as an input taken from provincial or national statistics, when available (see details in the following sections).

Individuals
Individuals are uniquely identified in the model with a unique identification number. The modelled characteristics of an individual include infectious status, age, level of immunity against the disease, and frailty. Immunity and frailty are modelled as real numbers between 0 and 1. They determine respectively the risk of acquiring the infection and being hospitalized as a result of infection.
Individuals are "linked" to social places. A link to a social place means that the individuals visit and have contacts in it (according to their schedule, see details later). All individuals are linked to a social place of type "household". Link to other social places depends on the characteristics of individuals: for example, all children between 5 and 18 years of age are linked to a school but none is linked to a workplace ( Figure 1).

in-silico Population construction
The population where the simulations will take place is defined iteratively. First, a large number of individuals are created and assigned to households according to the pre-specified distribution of household sizes. All households have at least one individual. Once all the households are populated, individuals are assigned an age based on pre-specified age distributions conditional on the household size. For example, a household of size one is associated with an age distribution that has a minimum value of 18 years old in order to avoid households populated with only one child. There is no individual without a link to a household.
For this study, the distributions of household sizes and population age profiles in Ontario were retrieved from Statistics Canada (2011 Census, catalogue 98-313-XCB2011022 for the household size and catalogue 98-311-XCB2011023 for ages). We did not have a direct access to the distribution of age of individuals conditional on their household size. However, we obtained the age distribution of the entire population, as well as the distribution of household sizes. Hence, it is possible to reconstruct (albeit, not in a unique fashion) the distributions of age of individuals conditional on their household size.
Let i,j be the distribution of age for the ith oldest individual in a household of size j. We chose to parametrize i,j as a Beta distribution because it allows to restrict ages between a lower-and upper-limit. This is useful, for example, when setting the age of individuals in households of size one, where the lower limit is 18 years and the upper limit is 99 years (hence no household of size one is populated with only a child). The parameters of Beta distributions are fitted such that, combined with the distribution of household sizes, they give the closest age distribution for the whole population (not conditional on household size) to the data. The fit is performed using an optimization procedure (function nlopt in R).
Construction of other social places is rather straightforward because there is no complex constraints on the age distribution of individuals linked to them. For example, individuals younger then 18 years of age are linked to schools (see section 3.2 for a detailed explanation) . The total number of a given social place and its size distribution are pre-specified as inputs. Then, individuals are assigned randomly to social places based on their schedule. Only social places tagged as "other" do not have individuals linked to them, because their purpose is to be visited randomly.
To summarize, the high level algorithm for the construction of the simulated population includes the following steps: 3. Link households with these individuals according to the distribution of household sizes.
4. Delete superfluous (not linked to any households) individuals.
5. Assign ages to individuals that are linked to a household, based on a pre-specified age distribution given a household size.
6. Create social places according to pre-specified number and size distribution.
7. Link individuals to social places according to their schedules.
The distribution of the sizes of the social places implemented are shown in Figure 2. They were informed from Statistics Canada for the households and workplaces, from the Toronto Transit Commission for the public transportation, from the Ontario Ministry of Education for the schools, and assumed to arbitrary values for the social places of type "other".

Model structure
The model implements, at the individual level, the natural history of influenza. Upon infection, a susceptible individual becomes exposed (E), i.e., infected but not yet infectious. After a period of time (drawn from the distribution of the latent period for each individual), exposed individuals become infectious and will either be symptomatic (I s ) or asymptomatic (I a ). Asymptomatic cases recover (R a ) without any disease outcomes such as hospitalization or death.  Figure 3. Compartmental representation of the epidemic model. Upon infection, individuals move from a susceptible (S) state to the exposed state (E) in which they are not infectious yet. Their infection is either symptomatic (I s ) or asymptomatic (I a ). Asymptomatic individuals will recover (R a ) after their infectious period. Symptomatic cases will either recover (R s ) or be hospitalized (H). Hospitalized patients will either recover or die from infection (D). Asymptomatic cases are not hospitalized. Depending on the intervention and its targeted population, vaccine (V ) and treatment (T ) may be administered at di↵erent stages of infection. The type of treatment for hospitalized cases is likely di↵erent (and much more intensive) than for non-hospitalized ones. So we considered that the treatment received in a hospital should not be added up with the antiviral one, which explains the absence of a treatment arrow in this diagram.
influenza vaccine, in addition to any other constraints from the chosen vaccination scenario. Finally, we do not include demographic variables of birth and natural death given the short duration of an influenza pandemic wave. Figure 3 illustrates the disease states.

Schedules and spatial movements
Every individual is assigned a daily "schedule" that determines which social places are visited at any given time of the day. There is a pre-specified probability that the visit actually occurs (based on the individual characteristics and epidemiological status, symptomatic infection may reduce this probability). For example, there is limited number of contacts if an individual's schedule allows mostly for movements between the assigned household and workplace. But social places like "other" or "public transportation" bring some additional mixing. An illustration of schedules is given in Figure 4.
We defined four types of schedules: "student", "unemployed", "worker" and "worker using public transportation". The schedule type is assigned based on the age of the individual. Individuals younger than 18 had a "student" schedule, those older than 65 had an "unemployed" one. Those of working age (between 18 and 65) are assigned either a "worker" or "unemployed" schedule. The later is chosen with a 0.10 probability. The proportion of workers using public transportation is set at 0.12 (of all the worker population) based on the num- Individuals with symptomatic infection may self-isolate, with a probability set at 0.9 [31], by staying in their household during their infectious period.

Frailty
Individuals have di↵erent levels of vulnerability to influenza infection which influences the risk of infection, severe outcomes, or even death. This vulnerability is represented in the model by a "frailty index" with a value between 0 and 1 (1 corresponds to the highest vulnerability) and depends on age. Individuals at higher risk of influenza-induced complications than the general population include those with chronic respiratory illness, cardiovascular diseases, immunosuppression, and diabetes [22]. Hence, in order to represent the correct age profile, the shape of the frailty index is fitted to survey data of chronic diseases (2014 Canadian Community Health Survey). Once the shape is fitted, the absolute level of the frailty index is fitted to the population hospitalization rate, which is specified in scenarios. We performed a segmented linear regression to fit the frailty shape f (a), as a function of age a (using the segmented package in R). Figure  5 shows the fit for frailty index.
For each individual in the simulation, the value of frailty index is drawn from a normal distribution with the mean f (a) (where a is the individual's age) and standard deviation arbitrarily set at 0.1. This deviation was chosen as the maximum di↵erence of frailty between consecutive 5-year age groups provided in data.

Humoral immunity
Humoral immunity reflects the pre-existing cross-reactive antibodies against the pathogen. The level of this immunity can determine whether or not an infection can successfully establish itself in the host. The level of humoral immunity depends on pre-existing cross-reactive antibodies as a result of prior exposures or vaccination. We label this level I h and model it as a real number ranging from 0 to 1, with 1 corresponding to full protection against influenza infection.
In the context of pandemic influenza, individuals will have little to no pre-existing cross-reactive antibodies, so the baseline protective e↵ect of humoral immunity I h is assumed to have a negligible value in our model, and we set I h = 0. However, the protective e↵ect of humoral immunity increases following vaccination. Vaccination leads to the generation of humoral immunity which reaches its maximum level over a 2-week period. This level is determined for each individual based on the frailty index and vaccine e cacy, as detailed in Section 4.3.

Cellular immunity
Cellular immunity is modelled as a determinant of symptomatic infection, and assumed to be age-dependent. Because this immunity "builds up" through life [13,12], older individuals are expected to have a higher level of cellular immunity than younger individuals. Thus, the protective e↵ect of cellular immunity for an individual of age a, I c (a), is modelled with a logistic shape curve ( Figure 6 where I ⇤ c is the maximum level of cellular immunity, and s and a pivot are shape parameters representing the slope and inflection point, respectively. Note that the values for I c are bounded by 0 and 1. The parameters s and a pivot are chosen such that the level I c is close to the immune  figure 4). Then, because our baseline pandemic scenario assumes less pre-existing cellular immunity than the 2009 H1N1 pandemic, we arbitrarily divide I ⇤ c by 2. The other parameters s and a pivot remain unchanged to keep the shape of the calibrated age-dependent cellular immunity I c .

Number of contacts
When an infectious individual is present in a social place, a pre-specified contact rate r c determines the number of contacts n c with other co-located individuals during the sojourn in that social place. For an individual of age a in a social place s, the number of contacts is drawn from a Poisson distribution with the mean drawn from a lognormal distribution.
The rate r c is log-normally distributed to allow for super-spreading events [20]. The mean of the log-normal distribution is parameterized based on the characteristics of individuals and the social places: it is modelled as the product of a baseline, constant, contact rate ⇤ and an adjustment factor (a, s). Parameter ⇤ is fitted to obtain a specific basic reproduction number R 0 for a given scenario. Values of R 0 (and hence ⇤) based on estimates from past epidemics and pandemics [1] are explored in scenarios to reflect di↵erent disease transmissibility and severity (see main text of article for scenario descriptions).
Unfortunately, there is no good agreement between empirical studies that investigated the contact rates among di↵erent age groups.
Face-to-face contact studies conducted in France [29,11,10,2] found that young children have a much higher contact rate than the rest of the population. Whereas a Canadian study [28] on the timing of laboratory-confirmed influenza infections across several seasons found that teenagers and young adults were infected earlier than other age groups (in particular young children), suggesting the age group 13-24 years old has a larger contact rate. Hence, we decided to fit the age component of (a, .) to Canadian study for our baseline scenario and as a sensitivity analysis, fit it to the French face-to-face contact study. We did not find any study that could reliably inform the relative contact rates between social places, so we set (., s) = 1 for all social places s.
The standard deviation is calculated such that the coe cient of variation is constant, arbitrarily fixed at 0.75 (di↵erent values are explored in the sensitivity analysis).

Contacts selection
Once the number of contacts from an infectious individual is determined, susceptible individuals in the same social place are selected randomly as contacts, with a constraint on age assortativity.
A predefined function '(x, y) determines the probability that an individual aged x will contact another individual aged y. In a given social place and for a given number of contacts, all susceptible individuals are assigned a "score" to be contacted based on the assortativity function '. The higher the score, the more likely the individual will be contacted. The selection is stochastic and a parameter ⇠ enables some deviation from the desired assortativity '. The scores are sorted and the susceptible contact is selected by choosing the rank of the sorted scores from an exponential distribution. If S 1 , S 2 , ... are the candidate susceptible individuals and u is the indices vector of the ranked scores, the selected susceptible will be S u[X] with X ⇠ Exp(⇠). The smaller ⇠, the more likely that the choice of the susceptible individual will comply to the predefined assortativity function '. If the selection of the ranked candidate contacts is deterministic (⇠ = 0), then age assortativity is too rigid, leading to unrealistic contact patterns.
The function ' was based on the POLYMOD study by Mossong et al [24]. This study was conducted in European countries and no equivalent study is available for Canada. Here, we implicitly assume that the age assortativity in contact patterns in Ontario is similar to that studied in European countries. The pre-specified function ' is shown in Figure 7.

Transmission probabilities
The probability that a given contact transmits the disease is reduced by a factor of ⇢ asympt if the infector is asymptomatic. We set ⇢ asympt = 0.1 [18,27] and explored other levels in the sensitivity analysis.
The probability that a susceptible contact acquires infection is a↵ected by the level of humoral immunity (I h ), which in our model is simply interpreted as the probability that infection fails to establish in the host. The probability to acquire infection, conditional upon an infectious contact, is therefore reduced by a factor

Asymptomatic and symptomatic infections
The symptomatic nature of the infection is determined probabilistically at the time of disease transmission. In the absence of vaccine, the probability of developing symptomatic infection is based on both the frailty index f and the level of cellular immunity I c : where µ is a parameter fitted to the asymptomatic fraction (i.e.., the proportion of asymptomatic infections among all infections) defined for a given scenario. Note that both the frailty and level of cellular immunity are age dependent.
Vaccination can also reduce the probability of developing symptomatic infection [3,6]. See section 4.3 for a description of the implementation of this mechanism.

Hospitalization and death
If an individual develops symptomatic infection, then there is a risk of hospitalization and subsequently, a risk of death. For any symptomatic infection, the hospitalization event is drawn from a Bernoulli distribution with the probability p hosp that depends on the individual's frailty f as follows: Parameter ↵ h a↵ects the overall number of influenza-induced hospitalizations and the frailty f a↵ects the age distribution of hospitalizations.
If hospitalization occurs, then the duration of the infectious period before hospitalization (i.e., the time to hospitalization) is sampled from a discrete distribution from Table 3 of Morrison et al. [23]. The duration of hospitalization is drawn from a lognormal distribution with the mean d hosp set at 12 days and variance 21 days 2 [16].
Death event conditional upon hospitalization is drawn from a Bernoulli trial with a probability that depends on frailty: Parameter ↵ d a↵ects the overall number of influenza-induced deaths and the frailty f a↵ects the age distribution of mortality.
The parameters ↵ h and ↵ d are obtained by fitting to the hospitalization and death probabilities specified in scenario assumptions. Ifp h is the pre-specified hospitalization probability for a given scenario, f (a) is the average frailty index for individuals aged a, and ⇡(a) is the proportion of individuals with age a in the whole population, the fitted value for ↵ h is given bŷ Similarly, for the probability of death conditional upon hospitalization, we calculatê The hospitalization and mortality probabilities chosen in our scenarios were based on the corresponding data obtained from Public Health Agency of Canada [25] Hospitalized individuals are moved to a linked social place labelled as "hospital" and stay in hospital for the duration of hospitalization (i.e., movements to other social places directed by their schedule are ignored during hospitalization).

Deployment
Antiviral treatment and vaccination are implemented as part of "interventions". There can be several simultaneous interventions of di↵erent nature (e.g., vaccination, antiviral treatment). All interventions are modelled as having the following features: • type: defines the nature of the intervention: treatment, vaccination, behavioural change (self-isolation).
• target population: defines the group of individuals that may receive the intervention. • time range: when the intervention begins and ends.
• administration rate: defines the fraction of population that receives the intervention during a given period of time.
• maximum coverage: specifies if there is a limit to the total number of individuals receiving the intervention.
Although the target groups in the population and administration rates are deterministic, the actual number of individuals receiving an intervention is stochastic. If T is the total number of individuals targeted to receive the intervention, r is the coverage, and dt is the simulation time step, then the actual number N dt of individuals receiving the intervention during the time step Moreover, the total number of individuals receiving the intervention is monitored at every time step of the simulation such that the constraint specified by the maximum coverage is fulfilled.

Antiviral treatment
Literature on the e↵ects of antiviral treatment (i.e., neuraminidase inhibitors) during the course of influenza infection indicates that it marginally reduces the duration of symptoms. It is unclear if antiviral treatment reduces the risk of severe complications, hospitalization or even death [15,14].
Hence, our model implements the e↵ect of antiviral treatment for symptomatic infection as reducing the infectious period. We considered a period of 12 hours as the minimum delay in start of the treatment following the onset of symptoms. This delay may reflect the time lag to seek care and receive antiviral treatment. Furthermore, we assume that individuals (who are not hospitalized) will not receive antiviral treatment if the time since the onset of symptoms has elapsed 2.5 days. This reflect the fact that antiviral treatment is most e↵ective if initiated within 2 days after the onset of symptoms. More specifically, if d is the infectious period drawn from the pre-specified distribution at the time of disease transmission, the reduction of infectious period ( ) due to antiviral treatment is modelled as a function of the time since symptom onset ⌧ (in days): (⌧ ) = (2.5 ⌧ )/2, 0.5 < ⌧ < 2.5 (11) This formulation indicates that if antiviral treatment starts at its earliest time (i.e., half a day) after the onset of symptoms, the infectious period reduces by a maximum of 1 day [15]. This reduction decreases with further delay in start of the treatment. In our model, antiviral drugs are considered for the treatment of symptomatic cases only, and not prophylaxis.

E↵ects of vaccine on disease transmission
Vaccination is modelled as potentially increasing the protective e↵ect of humoral immunity and decreasing the risk of developing symptomatic infection [5].
Following vaccination, each individual will have its humoral immunity level I h increased to 1.0 (its maximum value) with a probability p e↵ that represents the vaccine e cacy at the individual level. We model this e cacy as dependent on frailty: p e↵ (f ) = (1 f )" (12) where parameter " is the maximum vaccine e cacy for an individual that is fully immunecompetent. Values for " are set on a scenario basis within the estimated ranges in the published studies [7].
The increase in the level of humoral immunity following vaccination is not instantaneous, but increases linearly to its maximum level (i.e. 1.0) during a 2-week period [5]. We assumed this maximum level of protective immunity is reached with one dose of vaccine only. Influenza vaccination is assumed to have no e↵ect on individuals' frailty.
Vaccination is shown to reduce the risk of developing a symptomatic infection and severe outcomes, even if it fails to fully protect against infection [3,6]. We model this reduction based on the vaccine e cacy at the individual level [17]. Thus, the probability of symptomatic infection p sympt for a vaccinated individuals (defined in equation (5)) is multiplied by a factor 1 p e↵ (f ) if the individual is infected after the vaccine-induced immunity has reached its maximum level. See Figure 10 for an illustration of the vaccination process.

Vaccination strategy
We implemented a vaccination strategy inspired from the influenza pandemic planning strategy in Ontario (Ontario Health Plan for an Influenza Pandemic, 2013, chapter 7, immunization). This planning strategy assumes the vaccine is not available during the first months of the pandemic, whereas this study explores early (even before the first imported case) availability. The vaccination strategy implemented in the model prioritizes specific age groups and individuals in high-risk groups (i.e., patients with pre-existing chronic conditions as identified by the National Advisory Committee on Immunization). More specifically, the vaccination strategy is driven by the following algorithm: 1. The priority population group consists of individuals: • aged 5 years and younger • aged 65 years and older • with a frailty index f larger than the population average frailtyf 2. As long as the average vaccination proportion within the priority population group is below 35%, continue with prioritizing these groups 3. Once the average vaccination proportion within the priority population groups reaches 35%, other individuals can also be vaccinated. 4. The vaccination proportion is capped by age (even for the priority group) at pre-specified levels.The maximum proportion of vaccinated individuals by age were based on a study conducted in Ontario following the 2009 H1N1 pandemic influenza [9]: 45% for ages 12 and younger, 25% for ages between 12 and 30 years, 35% for ages between 30 and 55 years, and 55% for ages 55 years and older.

Vaccine administration rate
The rate of vaccination (in terms of vaccines administered per day) for seasonal influenza can be approximated by knowing the proportion of a population that is vaccinated and the vaccination roll-out period. For example, assuming that 30% of the population is vaccinated during a 3month period, the average vaccination rate is 300 per 100,000 people per day (30%/90 days). This approximation can be a starting point for the vaccination rate of a pandemic influenza. We investigated fractions and multiples of this level to reflect the shortage of vaccine supply or enhanced public health e↵orts in vaccine distribution.