Responding to Vaccine Safety Signals during Pandemic Influenza: A Modeling Study

Background Managing emerging vaccine safety signals during an influenza pandemic is challenging. Federal regulators must balance vaccine risks against benefits while maintaining public confidence in the public health system. Methods We developed a multi-criteria decision analysis model to explore regulatory decision-making in the context of emerging vaccine safety signals during a pandemic. We simulated vaccine safety surveillance system capabilities and used an age-structured compartmental model to develop potential pandemic scenarios. We used an expert-derived multi-attribute utility function to evaluate potential regulatory responses by combining four outcome measures into a single measure of interest: 1) expected vaccination benefit from averted influenza; 2) expected vaccination risk from vaccine-associated febrile seizures; 3) expected vaccination risk from vaccine-associated Guillain-Barre Syndrome; and 4) expected change in vaccine-seeking behavior in future influenza seasons. Results Over multiple scenarios, risk communication, with or without suspension of vaccination of high-risk persons, were the consistently preferred regulatory responses over no action or general suspension when safety signals were detected during a pandemic influenza. On average, the expert panel valued near-term vaccine-related outcomes relative to long-term projected outcomes by 3∶1. However, when decision-makers had minimal ability to influence near-term outcomes, the response was selected primarily by projected impacts on future vaccine-seeking behavior. Conclusions The selected regulatory response depends on how quickly a vaccine safety signal is identified relative to the peak of the pandemic and the initiation of vaccination. Our analysis suggested two areas for future investment: efforts to improve the size and timeliness of the surveillance system and behavioral research to understand changes in vaccine-seeking behavior.


Supporting Information 1 Age-Structured Disease Transmission Model
The regulatory decision model on vaccination policy required inputs on expected vaccination, vaccination-associated benefits, and vaccination-associated risks. To model these inputs, we adapted a published, age-structured disease transmission model [1]. Specifically, we added an influenza vaccination adoption function modeled as a Bass diffusion process [2]. We used the outputs of the adapted transmission model to seed the regulatory decision model.

Transmission Model Flow and Variable Definition
We used an age-structured disease transmission model to simulate the spread of pandemic influenza. A graphical display of the flow of persons is shown in Figure 1: Solid lines are transition paths between compartments. Dashed lines represent modes of influenza transmission. The top "track" represents unvaccinated individuals and the bottom "track" represents vaccinated individuals. The variables displayed above are written in vector notation, as each variable is a vector whose elements are age-specific.
The model variables were as follows, where the subscript i denotes the ith age group and the time unit t is in days: : the frequency-dependent influenza incidence rates for individuals in the ith age group at time t, also known as the force of infection; ! ( ) : the conditional probability that an individual in age group i will vaccinate at exactly time t since introduction of the vaccination, given that the individual has not vaccinated before that time; also known as the hazard function;

Additional Variables
!" : the average daily contact rate of an individual in the ith age group with an individual in the jth age group; !" ! : the per contact transmission probability of a susceptible individual in the ith age group by an unvaccinated infectious individual in the jth age group; !" ! : the per contact transmission probability of a susceptible individual in the ith age group by a vaccinated infectious individual in the jth age group; ! : the age-specific vaccine effectiveness against influenza infection for an individual in the ith age group; ! ! : the age-specific vaccine effectiveness against influenza-associated hospitalization for an individual in the ith age group; ! ! : the age-specific vaccine effectiveness against influenza-associated death for an individual in the ith age group; ! : 1/the average incubation time for an individual in the ith age group; ! ! : 1/the average duration of infectiousness for an unvaccinated individual in the ith age group; also known as the recovery rate; ! : the proportion of infected persons in the ith age group, also known as the attack rate; ! : the probability of vaccination in the ith age group.

System of Differential Equations
The model is a system of 15×n ordinary differential equations, with n age-groups and i=1, 2,…, n. The time unit t is in days. The model did not incorporate demographic processes because we considered the vaccination policy for a single pandemic influenza period. Additionally, the model did not incorporate concepts associated with pre-existing immunity or waning immunity. Finally, the model did not consider the effects of the use of antiviral medications or other non-pharmaceutical interventions in mitigating the influenza epidemic.
The model equations are as follows:

Calculation of Death and Hospitalization Rates
To calculate the influenza-associated death and hospitalization rates, we relied on estimation using information on the per capita death and hospitalization rates and concepts from survivability analysis [3].
where M i is the familiar case fatality ratio, given as influenza cases resulting in death divided by influenza cases.

Vaccination Adoption Function
The influenza vaccination adoption function was modeled as a Bass diffusion process [2] given by the following equation: where p i and q i are the Bass coefficients of innovation and imitation, respectively. Ω is the total population of those expected to vaccinate. In Bass diffusion models, all individuals that are expected to vaccinate will do so.
N is the total population at the start of the influenza period. θ i is the proportion of the total population in the ith age group.

Reproduction Number
We used the next-generation matrix technique [4,5] to derive ! , the basic reproduction number in our model. ! is the average number of new (secondary) influenza cases that an index influenza case generates over the course of its infectious period. The reproduction number is the spectral radius, or the largest eigenvalue, of the "nextgeneration" matrix FW -1 . That is, The next-generation matrix technique requires identification of m infected compartments, which contain x k infected individuals where k =1, 2,…, m. W k (x) : the net transfer rate of individuals in the kth compartment for reasons other than the incidence of infections; We employed matrix notation to describe these variables. ⊗ indicates element-byelement matrix multiplication.
In our model, all individuals who were not infected at t=0 were susceptible. ! was calculated at time t=0. Therefore in the next-generation matrix FW -1 , S R , S U , and S V were all calculated at time t=0. Because FW -1 has two rows of zero submatrices, ! was the spectral radius, or largest eigenvalue of the reduced FW -1 , shown below.
Algorithmically, in our model, we set ! at the desired level and then found β U and β V to achieve that level. For simplicity, we set β U = β V for all age groups, which implied that the probability of transmission of influenza was independent of both age and vaccination status.

Epidemiologic Parameters of Influenza
We ran two scenarios, utilizing the framework set forth by Reed et al. [6] that characterized influenza pandemics by both their transmissibility and severity. The first scenario had lower transmissibility and severity, and we referred to this scenario as the "mild" scenario. The second scenario was modeled after the 1957 pandemic, which had higher transmissibility and severity than the mild scenario but was not quite as transmissible and severe as the 1918 pandemic. We referred to this second scenario as the "severe" scenario.
For all scenarios, we assumed that the duration of infectiousness, influenza hospitalization rate and influenza death rate were independent of vaccination status (i.e., an ineffective vaccination did not impart partial protection against hospitalization or death). That is,

Severe Scenario
We set R 0 =2.0 in the severe scenario [19].

Vaccination Parameters
We estimated Bass coefficients p i and q i from data from the Post-Licensure Immunization Safety and Monitoring (PRISM) system [23]. See Section 3. We assumed these coefficients were reflective of the entire US population.
We assumed that vaccine effectiveness against influenza-related hospitalization and influenza-related death to be the same as vaccine effectiveness against influenza. That is, Notes: a The Bass p i and q i values were fit when t is measured in months. See Section 3 for more details.

US Demographics and Contact matrix
The US population was modeled using the US Census 2012 National Population Projections, Middle Series:  [26].
We used the 2013 projection, which lists the total US population as 316,438,601 persons ( ). ! is the probability of being in the ith age group.  Like Medlock and Galvani [1], we parameterized the contact matrix based on a study of daily contacts of eight European countries [27]. We used the data collected for physical contacts only. There were 15 age groups in the European data, divided in five-year blocks, i.e., ages 0-4, 5-9,…, 65-69, and 70+.
We adjusted the 15 age groups to 17 age groups more useful to our study: year old age groups. The contact rates from the European study's age group 0-4 were assumed to apply equally to our age groups 0-1 and 2-4: !! = !! and !! = !! , where subscripts 1 and 2 indicate age groups 1 and 2. Likewise, the contact rates from the European study's age group 15-19 were assumed to apply equally to our age groups 15-17 and 18-19.
The output of the contact matrix ( !" ) was a contact rate given in average contacts per day of age group i with age group j.
The final contact matrix is shown in Table 5.

Initial Conditions
We required an initial estimate of the infected population at the time of the model start.
Our model began August 1 in the year of the pandemic.
We began with an initial number of cases: ! at time t=0. We distributed these initial cases according to age based on the age distributions cited in the supplemental technical appendix in [18]. That is, we defined φ k to be equal to the case probabilities by the kth age group shown in Table 6 below. In order to further distribute these cases among the 17 age groups in our model, we rescaled as follows: We have already defined θ i as the probability of an individual being in the ith age group where the age groups and their probabilities awere defined in Table 4.
Let ! denote the membership of the kth age group where the age groups are defined as in Table 6 and k = 1,2,…,5.
Probability that an influenza case is in the ith age group= θ i y k All other compartments were initially zero.

Mild Scenario
For simplicity, we adopted the φ k values given in Table 6 despite the fact that they are based on data from the 2009-2010 pandemic. We set ! to be equal to 0.0001*( ) or ~32,000 persons. We primarily chose this value such that the peak of the pandemic influenza followed the peak of vaccination (see Figure 3 in the main text).

Severe Scenario
For simplicity, we adopted the φ k values given in Table 6 despite the fact that they are based on data from the 2009-2010 pandemic. We set ! at 100 cases so that the peak of vaccination aligned with the peak of influenza.

Delayed Vaccination Campaign Start
We employed a delayed start to the vaccination campaign such that the vaccination became available at some time after August 1, when t = t V , which we set to 30 days. At this time, we initiated the Bass diffusion process. Thus,

H1N1 Face Validation of Age-Structured Disease Transmission Model
We performed a simple face validation of our age-structured disease transmission model using similar assumptions to those used by [28].

Epidemiologic Parameters
We set = 1.55 in the H1N1 Face Validation scenario.

Demographics
The

Initial Conditions
We required an initial estimate of the infected population at the time of the model start.
Our model began August 1 in the year of the pandemic. As before, we set an initial number of cases at time t=0, and distributed these initial cases according to age based on the age distributions cited in the supplemental technical appendix in [18].
Shrestha et al. provided national estimates for the number of H1N1 cases in the entire month of August [16]: 1,605,760 cases. Truly, we would have liked to know the number of infected persons on August 1. We set ! to be equal to 650,000 persons. We chose this value such that the peak of the pandemic influenza with a population that achieved vaccination coverage as described in Table 9 occurred around November 1.
We employed a delayed start to the vaccination campaign such that the vaccination became available at some time after August 1, when t = t V , which we set to 60 days. This is consistent with the availability of H1N1 vaccine.   figure, we are happy to be this close. We were heartened to be so close on hospitalizations, which is the measure that CDC deems to be most reliable. On deaths, we are using the AHDRA reporting figures, which is an underascertainment, but CDC did not generate a multiplier to be applied to the AHDRA death totals as it did for hospitalizations [18]. Shrestha et al. apply a multiplier to the their hospitalizations figure [16]. If we do the same to calculate deaths, then we calculate 13,273 deaths, which is unsurprising given the correspondence among hospitalizations.

Vaccine Safety Surveillance System Model
We simulated sequential database surveillance for a) influenza vaccination-associated febrile seizures and b) influenza vaccination-associated Guillain-Barré Syndrome (GBS). We simulated surveillance by modeling an "enhanced" Post-licensure Rapid Immunization Safety Monitoring (PRISM) system, which is the supporting infrastructure for seasonal influenza surveillance [23]. Broadly speaking, we performed the following steps: 1 -Calculated database-specific chronological exposure estimates for each database using an age-structured delay differential equation model; 2 -Adjusted the chronological exposure estimates to reflect database-specific processing delays; 3 -Used the adjusted exposure estimates, along with other adjustments, to generate information on outcomes of interest; 4 -Aggregated database-specific information and performed sequential statistical analyses on simulated data.
A more detailed discussion of the methodology of these simulations is discussed elsewhere [30]. All analyses were completed using MATLAB® and R.

Bass Adoption Model for Influenza Vaccine
We used an age-structured delay differential equation model to describe the adoption of influenza vaccine within the surveillance population. Figure 1 shows the flow of persons.

Model Variables
: the number of potential adopters in the kth database in the ith age group at time t; : the number of adopters who have not yet entered the risk window in the kth database in the ith age group at time t; : the number adopters in the risk window ("at risk") the kth database in the ith age group at time t; : the number of adopters who have exited the risk window in the kth database in the ith age group at time t; : the number of total adopters in the kth database in the ith age group at time t; : the cumulative amount of exposed-person time in the kth database in the ith age group at time t; : the conditional probability that an individual in age group i in database k will vaccinate at exactly time t since introduction of the vaccination, given that the individual has not vaccinated before that time; also known as the hazard function; : the time period after a person has been exposed to a medical product but before the person is "at risk" for experiencing the outcome of interest, also known as the induction period or the latent period.
: the time period when a person is "at risk" of experiencing an outcome of interest following some exposure of interest, also know as the risk window.
: the incidence rate of the outcome expected under the null hypothesis (i.e., that there is no excess risk in the treatment group).

Model Equations
The model is a system of 7×n×m ordinary differential equations, with n age-groups from i=1,2,…,n, and m databases from k= i=1,2,…,m. The time unit t is in months. The model equations are as follows: where and are survivability constants that result from censoring individuals who experience the outcome of interest. The occurrence of the outcome of interest was modeled as a Poisson process.
is the total population of those likely to vaccinate.

Initial Conditions for Delay Differential Equation Model
P i k (t = 0) = Ω k ω i k (t = 0) = p i k All other compartments are zero at time t=0.

Calculating Database Lag
For each database k, it was necessary to delay the chronological exposure by a processing time, allowing for "claims" to be processed [31]. This time is a random variable estimated by empirical data.
Let be the cumulative distribution function of the processing time delay θ and for z=1,2,…,x where x is the number of discrete increments in .
It is important to ensure that t is in consistent units in the processing time delay and in the exposure vectors above. If that is not the case, both vectors need to be converted to a common unit.
Let M be a temporary matrix of exposures redistributed to account for the processing delay times where M will be size (w,x). M is then horizontally concatenated to a matrix of zeros of size (w, w-1). Now, each row of this new matrix, M hat, is then shifted right as in the example below to account for delays. In the example, w=4, x=2. 1=(1,1,...,1) is a vector of length w.
From ê i k , we can calculate the adjusted vector of cumulative exposure Ê i k which has length w+x-1.

Variable Definition
the number of true positive cases among the treatment population identified in the kth database for the ith age group at time t; the number of false positive cases among the treatment population identified in the kth database for the ith age group at time t; the number of false negative cases among the treatment population identified in the kth database for the ith age group at time t; Β i k (t) : : the incidence rate of the outcome expected under the null hypothesis (i.e., that there is no excess risk in the treatment population), also known as the background rate among the comparator population

Outcomes in the Treatment Population
Using the data on exposure, we generated 10,000 simulations of outcome patterns for analysis according to the following equations. That is, we repeatedly used a Poisson random number generator for variables distributed as Poisson.

Outcomes in the Comparator Population
We also generated outcomes under the null hypothesis, but these outcomes were deterministic and do not change throughout the 10,000 simulations:

Sequential Statistical Analyses
Given the simulated data on exposure and outcomes of interest, we applied the Poisson Maximized Sequential Probability Ratio Test (MaxSPRT) [34] to the simulated data, which was part of the current analytical plan when monitoring seasonal influenza vaccination [35]. Briefly, the Poisson MaxSPRT compared current adopters of influenza vaccination to historical adopters of influenza vaccination to detect an elevated risk of designated adverse events (i.e., outcomes of interest). The null hypothesis was that the risk after influenza vaccination in the pandemic period is no greater than the risk after influenza vaccination in past seasons.
The null hypothesis of no increased risk will be rejected if the test statistic, the log likelihood ratio (LLR), reaches an upper boundary, known as the critical value. The null hypothesis will not be rejected if the total number of cases of the designated adverse event surpasses a pre-specified "upper limit" for surveillance. Alternatively, it was possible to fail to reach either boundary condition because of a lack of cases, which ultimately implies a failure to reject the null hypothesis.
For each outcome of interest, the critical value of the LLR was calculated by the userspecified upper limit of expected events and alpha level (i.e., Type I error). Upper limits in these analyses are historically selected based on the approximate number of events that would be expected under the null hypothesis in the risk interval. This value is incremented up to prevent reaching an end of surveillance before the surveillance period is over. One-tailed tests were used for alpha (i.e., Type I error).

Parameters
All surveillance parameters were modeled on seasonal influenza surveillance planned for the 2013-2014 season [35].

Surveillance Parameters
The comparator outcome rate was based on data on individuals in the PRISM system that had previously received an inactivated trivalent influenza vaccine and experienced a febrile seizure in the risk window. All other epidemiologic parameters are based on the protocol or discussion with PRISM epidemiologists.

Database Parameters
The size of the surveillance population was determined by averaging adopters of influenza vaccination in the PRISM system over three previous influenza seasons (2008-2009, 2009-2010, 2010-2011). The processing delay was an empirical distribution function based on historical claims originating from the Emergency Department.

Values References
Database Size  a Bass p and Bass q are listed by age group, but they are also calculated by database. We are not permitted to share database-specific information and so a PRISM-wide coefficient is shown. Bass p and Bass q were fit when t is measured in months.

Surveillance Parameters
The comparator outcome rate was based on data on individuals in the PRISM system that had previously received an inactivated trivalent influenza vaccine and experienced a febrile seizure in the risk window. All other epidemiologic parameters were based on the protocol or discussion with PRISM epidemiologists.

Database Parameters
The size of the surveillance population was determined by averaging adopters of influenza vaccination in the PRISM system over two previous influenza seasons (2009)(2010)(2010)(2011). The processing delay was an empirical distribution function based on historical claims originating from inpatient services at the hospital.

Values References
Database Size

Additive Multi-Attribute Utility Function
In multi-criteria decision analyses (MCDAs), outcomes of decisions are characterized by multiple, often competing criteria or attributes, some of which may be less important than others but nevertheless important enough to be considered in the analysis. Outcomes are thus multivariate. Decision makers must consider all criteria when evaluating possible decision options. A multi-attribute utility (MAU) function is a mathematical function used to characterize the overall value (or "utility") of each outcome relative to the others to represent the judgment of the decision maker [38]. The MAU is used to weight and combine multiple attributes into a single figure of merit whose expected value was maximized to select the preferred regulatory decision.
Using an expert panel preference solicitation process, we developed an additive MAU function. We also considered multiplicative and multi-linear formulations but ruled these out in favor of the simplicity of the additive model.
We use the following functional form: where x i are the attributes of interest and k i are the weights associated with those attributes.
We evaluated four attributes as follows: Every outcome of the decision model is represented as a vector of values for each of these four attributes. One example outcome follows: Attribute 1: 23,940,000 saved infections, 18,000 saved hospitalizations, 9,794 saved deaths Attribute 2: 6700 excess febrile seizures Attribute 3: 0 excess GBS cases Attribute 4: Minor Change, 10% reduction in vaccine-seeking behavior

Expert Panel Preference Elicitation
We convened an expert panel of six physicians who were currently serving or had previously served on vaccination-related federal advisory boards. All physicians were pediatricians. They were asked as individuals to assume the role of regulatory decision makers, i.e., their values served as proxy for decision makers may have to actually make such decisions.
We used a fractional factorial design to select sets of outcomes (i.e., combinations of levels of the attributes listed above) for evaluation; each panelist was asked to consider each outcome as a whole and to rank the set of outcomes, then to assign the best outcome (in that panelist's assessment) the value of 100 and the worst outcome a value of 0, then to rate the remaining outcomes relative to the anchor points. This process and facilitated discussion of ratings among the panelists was to help them think deeply about how each outcome attribute weighed in the decision process and what tradeoffs were made in thinking about overall attractiveness or unattractiveness of an outcome. Later in the day, using a series of steps [38][39][40] each panelist was led through a process of developing the disaggregated parts of a MAU function which could be used to assign values to outcomes.
We performed separate elicitations conditional on the mild scenario and on the severe scenario, thereby creating separate MAU functions for the two situations. We combined each expert panelist's MAU function to derive the "average" federal decision-maker.

Attribute
Expected Notes: A scaling constant represents the relative weight given to each criterion in the utility function and together, they must sum to 1. a The mild scenario was characterized by low transmissibility, low severity, and the peak of vaccination preceded the peak of influenza transmission.
Following the initial day-long expert panel preference elicitation process, we re-contacted each expert with an analysis of the correlation between their user-valued ratings and the model of their individual preferences to assess whether we had accurately captured such preferences. We provided each expert the opportunity to change or iterate on his or her individualized model. Most experts retained their original preferences. Figure 3 shows the correspondence between the user-valued ratings and the derived model in the set of 31 mild scenario permutations. The upper panel shows the utility computed by the MAU in blue compared to the utility specified by the expert panel (green). In the lower panel, the black line is perfect correlation between the model and expert panel ratings. The actual correlation coefficient of the model and the expert panel ratings is displayed, and the blue line is the regression.  The upper panel shows the utility computed by the MAU in blue compared to the utility specified by the expert panel (green). In the lower panel, the black line is perfect correlation between the model and expert panel ratings. The actual correlation coefficient of the model and the expert panel ratings is displayed, and the blue line is the regression.

Conditional Probability Elicitation
Our regulatory decision model is concerned with the general public's reaction to federal vaccination policy changes and how that reaction translates into future vaccination seeking behavior. We lacked any public preference surveys similar to those described in [41] regarding changes in vaccination-seeking behavior as a result of government policy and therefore, we decided to elicit conditional predictions of the likelihoods of various reactions from the expert panel.
For each of the four decision options examined in our regulatory decision model, we asked the question: "We are experiencing a {Mild/Severe} pandemic and an information signal of vaccineassociated risk is received by federal decision-makers, and they select {Decision Option 1-4}, then what is the probability that: a. Public does not change their future vaccination-seeking behavior. b. Public responds by reducing their future vaccination-seeking behavior by 10%. c. Public responds by reducing their future vaccination-seeking behavior by 25%."   Notes: This is the averaged probability elicited from the expert panel. If a particular regulatory response is selected in the model, then each row represents the probability of the three levels (i.e., note that each row sums to 1.0). Therefore, regulatory responses with the highest probability of no change are associated with the highest levels of future vaccine-seeking behavior. Notes: This is the averaged probability elicited from the expert panel. If a particular regulatory response is selected in the model, then each row represents the probability of the three levels (i.e., note that each row sums to 1.0). Therefore, regulatory responses with the highest probability of no change are associated with the highest levels of future vaccine-seeking behavior.