Fast and accurate dynamic estimation of field effectiveness of meningococcal vaccines

Background Estimating the effectiveness of meningococcal vaccines with high accuracy and precision can be challenging due to the low incidence of the invasive disease, which ranges between 0.5 and 1 cases per 100,000 in Europe and North America. Vaccine effectiveness (VE) is usually estimated with a screening method that combines in one formula the proportion of meningococcal disease cases that have been vaccinated and the proportion of vaccinated in the overall population. Due to the small number of cases, initial point estimates are affected by large uncertainties and several years may be required to estimate VE with a small confidence interval. Methods We used a Monte Carlo maximum likelihood (MCML) approach to estimate the effectiveness of meningococcal vaccines, based on stochastic simulations of a dynamic model for meningococcal transmission and vaccination. We calibrated the model to describe two immunization campaigns: the campaign against MenC in England and the Bexsero campaign that started in the UK in September 2015. First, the MCML method provided estimates for both the direct and indirect effects of the MenC vaccine that were validated against results published in the literature. Then, we assessed the performance of the MCML method in terms of time gain with respect to the screening method under different assumptions of VE for Bexsero. Results MCML estimates of VE for the MenC immunization campaign are in good agreement with results based on the screening method and carriage studies, yet characterized by smaller confidence intervals and obtained using only incidence data collected within 2 years of scheduled vaccination. Also, we show that the MCML method could provide a fast and accurate estimate of the effectiveness of Bexsero, with a time gain, with respect to the screening method, that could range from 2 to 15 years, depending on the value of VE measured from field data. Conclusions Results indicate that inference methods based on dynamic computational models can be successfully used to quantify in near real time the effectiveness of immunization campaigns against Neisseria meningitidis. Such an approach could represent an important tool to complement and support traditional observational studies, in the initial phase of a campaign. Electronic supplementary material The online version of this article (doi:10.1186/s12916-016-0642-2) contains supplementary material, which is available to authorized users.


S1 Estimating the vaccine effectiveness with the screening method
Vaccines are usually evaluated by means of a number, called vaccine effectiveness (VE), calculated as one minus a measure of relative risk (RR) [1]: The relative risk is a ratio of two risks: the risk that a certain event occurs in a group of people divided by the risk that the same event happens in a second distinct group of people. In this context, the two groups are composed by vaccinated and unvaccinated individuals respectively, while the event of interest is some observable related to a particular infectious disease, for instance clinical observation of illness, possibly supported by biological tests to confirm the infectious agent, or sampling of pathogens in asymptomatic individuals. The risks are probabilities, often approximated by rates.
In case-control studies [2], cases are ascertained, and their vaccination statuses are retrospectively observed. Then, a random sample of the same population to which the cases belong (the controls) is drawn to give information about the proportion of vaccinated persons. Thus, cases in case-control studies represent a sampling fraction of all cases and controls represent a different sampling fraction of the healthy susceptible population. In order to avoid possible potential biases, the populations should share epidemiologically relevant characteristics, such as sex, age or geographical area, depending on the disease 1 under study. The ratio of vaccinated over unvaccinated cases is divided by the ratio of vaccinated over unvaccinated controls. This is not a ratio of two risks, but rather an odds ratio (OR). However, for rare disease, OR RR: where p v and p u are the probabilities of disease (or risks) in vaccinated and unvaccinated groups, d v and d u are the respective numbers of observed cases. n v and n u are the population in the two groups, so n = u + h with h v and h u are the controls, healthy individuals that did not experienced the disease or any other event of interest. The odds are defined as A variant of the case-control studies is the screening method [3]. As in case-control studies, it requires to compare the proportion of vaccinated among the cases and the control population. However it differs from other retrospective methods since control is achieved through external data: the vaccine coverage is estimated from a source that is external to the study, that has to be highly accurate.
Let's imagine to estimate the vaccine effectiveness by directly comparing the proportion of cases occurred in the vaccinated group (PCV) with the proportion of the total population that has been vaccinated (PPV). Then, Equation S2 can be rewritten as [3,4]: Instead of drawing a sample of the population and estimating PPV from it, as in casecontrol studies, the screening method employs the assumed true value of PPV from monitoring systems. Employing disease occurrence data, it is not possible to detect indirect effects of vaccines. An indirect effect would reduce the number of disease cases in all the population, independently on they vaccinal status. Therefore, PCV, that is a number of cases in vaccinated divided by the overall number of cases, would remain unchanged, and so would the VE. Moreover, as the indirect effectiveness increases, the samples of observed cases decrease, making the VE estimation always more difficult and even impossible, in the extreme case of perfect immunization against transmission that would imply null samples sizes.

S1.1 Screening method's sample size analysis
Estimates of vaccine effectiveness must be accompanied by a measure of their precision (or error), usually expressed as a 95% confidence interval (CI). The precision increases as the study sample increases. When designing effectiveness studies, it is crucial to calculate the sample size of cases needed to achieve a desired precision. For the screening method, Farrington [4] derived a formula to calculate the number of cases x to be observed to estimate a hypothetical true value of VE true with a lower 95% confidence limit above a certain threshold VE low , depending on the proportion of the population that has been vaccinated (PPV): Table S1 shows the sample sizes x calculated for different values of PPV, VE true and the distance (VE true -VE low ) of lower 95% confidence limit from the true effectiveness.

S1.2 The case of Bexsero
It is convenient to convert sample sizes into observational times, by employing the disease incidence observed in the past in the cohorts of interest.  In order to convert sample sizes into observational times, we employed recent reported incidences, assuming that future incidence will be the same if no vaccination is implemented. The purpose is to provide a conservative estimate of the observational time as the MenB disease cases have occurred with a decreasing trend during the last years. The total number of MenB disease cases reported during the last three epidemiological years for which these numbers are publicly available, i.e. 2011/12, 2012/13 and 2013/14, were 587, 595 and 424, respectively [5], corresponding to an average 535.3 cases per year. The effectiveness study will take into account only infants of one year of age, belonging to the fully vaccinated cohort. Public Health England does not provide data with such a fine age-stratification. Thus, we estimated the number of cases z past a=1 occurred in the age class of interest by multiplying the mid 2011 to mid 2014 average of the yearly total number of cases, k, times the normalised distributions of MenB cases by age observed in 2007 and in 2008, here denoted with y a=1 (data available from MATS databases [6] and stratified by single year): It resulted that z past a=1 = 70.2 cases per year.
Let's now assume that, if no vaccination is implemented, during the 2015/16 and the following epidemiological years, the expected number of MenB disease cases in infants of 1 year of age will be: In this case, it is possible to calculate also an expected yearly number of cases z a=1 that will be reported in case of vaccination, providing hypothetical values of VE true and PPV: Diving the sample size x calculated by fixing VE true , VE low and PPV (Equation S5) by the appropriate guess of number of cases z a=1 expected during one year, it is finally possible to translate sample sizes in observational times T :

S2 Dynamic model for meningococcal carriage and disease
The dynamic model is structured into two parts: 1. a transmission model ( Figure 1A in the main text), that, depending on the vaccine effectiveness values given as inputs, simulates the asymptomatic spreading of meningococcal carriage and the immunization campaign implementation in a specific population. The model output is the number of asymptomatic infection events.
2. An observational model ( Figure 1B in the main text), that allows to build likelihoods given time series of IMD cases notified by public health authorities and the transmission model outcomes.

S2.1 Transmission model
The following notation is used to the define the compartments of the transmission model: • S a = susceptible unvaccinated of age class a, • C a = carriers unvaccinated of age class a, • SVI a = susceptible vaccinated and immune of age class a, • CVI a = carriers vaccinated and immune of age class a, • SV a = susceptible vaccinated but not immune of age class a, • CV a = carriers vaccinated but not immune of age class a.

Transition
Description Parameter waning of immunity ω a and τ imm CVI a → CV Table S3: Allowed transitions between compartments.
The population is divided into 160 age groups of three months size, plus one group that includes all individuals ≥ 40 years old. Thus, the index a ranges from 0 to 160.
Individuals are moved between compartments through three broad types of transitions: demographic, vaccinal and epidemiological. All the possible transitions are summarized in Table S3.

S2.1.1 Demographic transitions
The demographic transitions are births, ageing and deaths. Births consist in introducing new susceptible individuals into the first age class. Deaths remove individuals from all compartments, with a mortality rate µ a that depends on the age class a. Ageing consists in moving all the subjects belonging to a compartment of a given age class to the same kind of compartment of the following age class, except for the last age class, from which subjects can only exit when they die.

S2.1.2 Vaccinal transitions
Vaccinal transitions (dashed arrows in Figure 1A in the main text) act according to vaccine schedules. Individuals are removed from compartments of unvaccinated according to the reported fractions of individuals covered by the campaign, namely γ a . The precise time points t vac at which individuals are moved depend on the vaccine and the vaccination strategy. Individual can be moved only if they have received all the recommended doses and only after a time needed to enhance an immune response. Vaccinated individuals are divided into two classes: vaccinated immune and vaccinated non-immune. Here the direct vaccine effectiveness VE dir comes into play. Only a fraction VE dir of the vaccinated gains immunity to IMD, the remaining 1-VE dir goes to non-immune compartments. However, after a certain time period, vaccines cease to be efficacious, so that vaccinated subjects lose their acquired immunity.
The waning of immunity is modelled by means of transitions SVI→SV and CVI→CV. Such transitions are activated after a time span τ imm following the precise time of vaccination t vac . The duration of immunity distributions are usually right-skewed, i.e. a fraction of subjects remains immune for a longer time if compared to the average. To take into account this feature, we model a fraction ω a of immune individuals that lose their immunity after a time τ imm , so that the remaining 1 − ω fraction of the immune individuals have a long-lasting immunity.

S2.1.3 Epidemiological transitions
The epidemiological transitions (infection and recovery) are responsible for the transmission of meningococcal carriage in the population. These transitions are stochastic in order to incorporate intrinsic noise. Therefore, the corresponding parameters can be interpreted as probabilities. They are obtained from rates, but the model's time scale has been chosen so that rates are a good approximation of probabilities [7,1]. Specifically, all the stochastic processes are modelled as binomial processes.
All the carriers have the same chance ρ of recovering, regardless of their age and vaccinal state. If not immune, susceptibles can become carriers with a probability λ a . Individuals in SVI compartments are instead subjected to a probability of becoming carriers reduced by a factor VE ind , if compared to individuals in S and SV of the same age. The force of infection λ a is defined as: where N a (t) = C a (t) + CVI a (t) + CV a (t) + S a (t) + SVI a (t) + SV a (t) is the total population. The matrix elements m a,a represent the number of contacts a subject of age a has with individuals of age a during a time-step. Multiplied by the prevalence of carriage by age, it provides the pattern of potentially infectious contacts. This is then rescaled by β a , that incorporates an age-dependent risk of meningococcal transmission and nasopharingeal colonization given a contact with an infected person. At each time step, the number of individuals at risk of developing invasive disease, namely subjected to transitions S a → C a and SV a →CV a , are counted. These counts are depicted in Figure 1A of the main text by the red swollen arrows. They are denoted with J a for non-vaccinated, and JV a for vaccinated not immune. The Js and the JV s are the transmission model's outcomes, to be connected to the time series of laboratory confirmed reported IMD cases, that are the system's observable.
The progression to invasive meningococcal disease is not part of the transmission model. In fact, its contribution to the spreading of meningococcal carriage is negligible, since it affects a very small fraction of the susceptible S and SV that have been infected. The risk of IMD given infection will be instead used in the observational model.

S2.2 Observational model
The observational model described in the main text ( Figure 1, panel B) highlights the link between the transmission model outcomes J a and JV a and the time series of IMD cases by age and vaccination status to be collected from the monitoring systems, such as Public Health England. We explicitly modelled an age-dependent probability θ a of observing an IMD case given an underlying, not directly observable (but reproducible through the transmission model) infection event. We assumed no under-reporting, since IMD cases are so severe that is quite unlikely that they go unnoticed. For what concerns the laboratory confirmation (included serogroup ascertainment), we assumed that there is no bias across ages in this process.
The progression-to-disease process is modelled through a binomial distribution with parameters J (JV for the vaccinated non-immune) and θ. Therefore the number of IMD cases of age a observed during a certain time τ are distributed as : DV obs a (τ ) Bin(JV a (τ ), θ a ).
(S12) θ a is an age-dependet risk of disease given infection multiplied by a probability of laboratory confirmation and by a probability of reporting. Since we assume a perfect reporting, the latter equals 1. The probability of laboratory confirmation is a constant, whose precise value is not relevant to our aim. Therefore the θs can be called risks of disease given infection, but implying that they include also a probability of laboratory confirmation.

S2.3 Transmission model's equations
The equations governing the model dynamics are many and non-linear. Here we provide the deterministic skeleton of the transmission model's system of equations.
Equation without vaccinal transitions. When a = 0 (note that there are no carrier newborns): 8 When 0 < a < 160: When a = 160: (S18) N a = S a +C a is the total population of age a, while λ a depends on two additional parameters that allow to model the meningococcal transmission patterns: β a is a rate of meningococcal transmission and colonization given contact, while m a,a is the number of contacts that an indicidual of age a has with an individual of age a during a time-step.
Complete system of equations. When a = 0 (note that there are no carrier newborns, nor vaccinated newborns): When 0 < a < 160: The complete set of deterministic equations governing the compartments of the last age class are not shown, since they are very verbose, but they can be easily derived from the above equations, as done for equation S17 and S18 . Here, N a = S a +C a , NVI a = SVI a +CVI a , and NV a = SV a +CV a . As displayed in the main text, λ a (t) is defined as: with the same β a and m a,a defined after equation S19. Note that in this deterministic formulation, here provided for illustrative purposes, the compartment populations and the numbers of cases are not integer numbers. Differently, in the stochastic model that we used, these numbers are naturally discrete when sampled from binomial distributions, or are forced to be discrete by rounding when they emerge from deterministic multiplications, as in the demographic transitions.

S2.4 Observational model's equations
It is possible to add to the previous systems of difference equations also the deterministic version of the observational model: with D and DV number of IMD cases emerging from the numbers of meningococcal transmission events J and JV counted during the time-step t.

S3 Model parametrisation
The meningococcal model is a piece of an inferential framework to estimate the effectiveness of a vaccine in the field, once an immunisation campaign is introduced. Given the quantitative nature of the objective, the particular setting of the immunisation campaign is of primary importance. The case studies on which the model has been parametrised are two: 1. The serogroup C meningococcal vaccine introduction in England, started in November 1999; 2. The serogroup B meningococcal vaccine introduction in England, started in September 2015.
Even if the two vaccines under study act in immunologically very different ways, from the point of view of the model they are analogous. Differences lie in the precise values assigned to parameters and initial populations assigned to the compartments, that can be determined by the available data. The two implementations of the model will be hereafter sometimes referred to as MenC model and MenB model.

S3.1 Initial compartments
First of all, the compartments must be filled with the appropriate numbers of individuals, i.e. initial conditions have to be set. This task was accomplished using three ingredients: • the England's population by year of age estimated at the time point of interest, immediately before the campaign start; • the average meningococcal carriage prevalence by year of age; • the fractions of carriers of serogroup B and serogroup C.
Population estimates are available from the UK Office of National Statistics [8] but they have a minimum resolution of one year. In order to fit the model's age classes, each 1-year-wide sub-population up to 39 years old is divided by four, while all the subjects older than 40 are aggregated into the last age class. The time points chosen for calibration are mid 1999 for the MenC model and mid 2014 for the MenB model.
The age pattern of carriage prevalence that we employ comes from of a recent metaanalysis, that merged several meningococcal carriage studies in Europe undertaken since the '70s [9]. As for the population, also the carriage prevalence has one year of resolution, so it is opportunely fitted to the age classes through cubic interpolation. The interpolation is done by fixing a boundary condition of null prevalence for negative ages.
The meta-analysis does not differentiate between serogroups, so it is rescaled using another source of data [10], that reported the results of three large large carriage studies, Thus, the initial populations belonging to the compartments are:  Figure S1 shows the population distributions employed for the two case studies.

S3.2 Parametrisation of the transitions
Once the initial conditions are fixed, we set the parameters that regulate the transitions between compartments. Some of these quantities are known in literature, for instance demographic rates. In other cases, some calculations or additional assumptions are required to determine unknown quantities such the βs or the θs. A summary of the quantities needed to parameterise the model is shown in the main text, Table 1.

Births and deaths
The number of susceptible newborns introduced in the model at each quarter of year are generated by multiplying the crude birthrate α by the total population. α was obtained by dividing by four the yearly crude birthrates estimated in the years of interest in England and Wales [8]. For the MenB model, α is fixed to the most recent report, that is 12.3 births per 1000 population in 2013.
The mortality rates µ a are derived from the notified deaths by age in England [8], that are firstly divided by the population, and then linearly interpolated to fit the model's age classes. Figure S2 shows the resulting mortality rates, plotted in logarithmic scale.

Recovery rate
The recovery rate ρ is taken as the reciprocal of the average duration of carriage τ car , so that ρ = 1/τ car .
There are only a few longitudinal studies reporting the duration of meningococcal carriage, a quantity not easily obtainable and subject to large fluctuations. Median values reported in literature range between 4.1 and 10.2 months [11,12,13,14] and no significant differences across age, serogroup and vaccination status are known. During the last years, modellers have mainly used a τ car of six months for MenB [15,16]. 6 months correspond to exactly 2 time steps of the model, thus ρ = 0.5 per time step. 6 months is here taken as the base case for duration of carriage, but other values of τ car ranging between 4 and 10 months have been tested, as a sensitivity analysis (see section S5).

Infection rates
The λs regulate the infection transitions, but they depend on two other parameters (see with carriers, β. The m aa parameters correspond to the number of persons belonging to the age class a , that one person of age a is in contact with, during three months (one time step) in England. They are rearranged from the contact rates φ ij surveyed in Great Britain during year 2005 (physical and non-physical contacts) [17]. Such contact rates are stratified to form 15 age classes (14 classes having 5 years of width, from 0-4 to 65-69 years, plus a 70+ years class), and they are daily contacts.
To get the m matrix, the φ matrix has been: (1) multiplied by the 2005 population in England arranged in the φ's age classes, to get the total number of contacts between ages; (2) grouped by ages higher or equal to 40 years; (3) divided by the 2005 population in England, also grouped for 40+ years old, to get back the number of contacts per single subject; (4) temporally rescaled to get contacts per time step of three months; (5) arranged in the a classes of the model, using a two-dimensional linear interpolation technique. The interpolation is done by fixing the boundary condition of null contacts for negative ages. The resulting contact pattern is shown in Figure S3.
The meningococcal asymptomatic transmission (or acquisition) given contact rates by age, that we called β a , are still unknown in the literature. In order to parametrize them, we imposed the time-stability of the carriage prevalence, adding the following equation: corresponding to a conservation of the prevalence c a = C a /N a for each age a (N a is the total population). Mathematically speaking, it consists in adding an equation containing no further variables that, acting as a constraint, allows to close the system of equations from S13 to S18. However, this assumption is not only a mere mathematical trick, but it is motivated by the fact that the prevalence of meningococcal carriage has a shape that is with a certain degree conserved in time [9]. Moreover, since our aim is building a tool for a fast inference of effectiveness, it seems reasonable to suppose that large fluctuations of the prevalence in short times (i.e. in less of 5 years) are unlikely, given the long duration of meningococcal carriage, and the slow dynamics of infection spreading. More detailed studies of meningococcal carriage would be precious in order to quantify transmission rates by age, or to verify and possibly modify our assumption. After constraining equations from S13 to S18 with the condition S38, it is possible to calculate β a . If a is not the last age (a < 160), for any time t, where N −1 = N newborns = α a N a and C −1 = C newborns = 0. For the last age class, a = 160, .

(S40)
The values of β calculated with the above relationships in the MenC and the MenB settings are very similar (see Figure S4).

Risks of disease given infection
The risk of disease given infection by age, θ a , is calculated as where D exp a (t) is the expected occurrence of invasive disease during a time step if no vaccination is implemented, and J exp a (t) is the number of expected infection events produced by the model during a time step if no vaccination is implemented. The calculation is based on the most recent available data before vaccination.
Specifically, for the MenC model D exp a was obtained by multiplying the normalised incidence of disease cases by age reported from January 1998 to June 1999 times the total number of cases observed in the epidemiological year 1998/1999, and then divided by 4 to get the desired time scale. However, we also tested different total numbers of cases (see section S5). For the MenB model the IMD incidence distribution employs cases reported in 2007 and 2008 in England and Wales, times the yearly average of the total number of cases reported in England from mid 2012 to mid 2015, then divided by 4 to get the desired time scale. The J exp can be calculated using the infection term of the model equations, set at the initial time step: The θ rates do not depend on time, and can make sense only to reproduce meningococcal disease in regions where it is endemic, as in the UK: their calibration is based on the additional assumption that the the risk of developing invasive disease that susceptibles of a certain age will experience once infected will not change in the future, when other susceptibles of the same age will acquire the pathogen.

Vaccinal transitions
γ a and ω a are the parameters controlling the vaccination and the waning of immunity, respectively. They are modelled as on switches that, when activated, move a certain fixed fraction of a compartmental population.
For the MenB model, γ has in fact a single value, because the modelled immunisation campaign targets a single age class of infants, that are continuously targeted by a routine vaccination. The value of γ must be assumed because coverage data are not available yet. We assume here a coverage γ = 95% for the MenB vaccine implementation. Only children that have received all the three vaccine doses are considered vaccinated, therefore, during each time step, a fraction of unvaccinated individuals in the age class a = 4, corresponding to 12-14 months of age, are moved to the respective compartments of vaccinated.
The MenC immunisation campaign covered a much larger population. It has been composed by a large catch-up campaign, targeting individuals from 4 months to 18 years of age, extended later to cover up to 24 years old, plus a routine campaign designed with a 2+4+6 months schedule. In this case, coverage data are available in literature [18].
The fraction ω of immune individuals for which, after a certain while, the immunity wanes can be obtained from serological studies. Here we adopt the simplifying assumption that after a certain time τ imm a fraction ω of the individuals belonging to a compartment of vaccinated immune is moved to a vaccinated non-immune compartment. This fraction is assumed to be ω = 100% for the MenB vaccine, while the assumed duration of immunity is τ M enB imm = 18 months after boost. For MenC, we used results of serological studies that estimated that a fraction f imm = 0.416 of those who initially gained immunity was still immune 3-4 years after the last dose [19]. Therefore ω a = 1 − 0.416 = 0.584 for MenC. τ M enC imm was fixed to 15 months for routine vaccination, 5 years for catch-up of one year old infants, 10 years for all the others [20,21].

S3.3 Sequential Monte Carlo
Both the VE are varied between 0 and 1 with a step of 0.001, so that the (VE dir , VE ind ) space of the unknown parameters is a discrete grid. For each combination of the VEs, here we explain the process of likelihood building.
1. 1000 particles are generated. Each particle represents the variables of the systems, i.e. the populations belonging to each of the 966 compartmental populations. At the beginning, all the particles are set to the same compartments' initial conditions.
2. Then all the 1000 particles evolve for the first year y (4 model's time steps) in the transmission model. At the end of the process, they will have in general different values in the compartments' space, because of the stochasticity of the model. In addition to the updated state, each particle produces also an output: the Js by age and vaccination status.
3. The J a,v (y = 1) values are combined with the risks θ a (both opportunely rescaled from age classes a to age groups b if needed) in the observational model, to calculate a binomial probabilities to observe D obs b,v (y = 1). These binomial probabilities multiplied together are the likelihood at y = 1. The particles are then pruned according to their specific likelihoods, using a weighted resampling algorithm [22], and, finally, the probabilities of the particles are averaged and a likelihood for the first year is obtained: Then, this procedure is repeated from point 2, substituting "first year" and "y = 1" with "second year" and "y = 2". The result is a second likelihood, conditional to the states and observations of the first year: Going on for the following years, using the Markov property, all the conditional likelihoods are multiplied together to eventually get the likelihood shown in equation 3 of the main text. Once the final likelihood at year Y is calculated for the particular set of VEs chosen at the beginning, then the VEs are varied, a new set of 1000 particles is created to finally produce a likelihood for another point of the VE's grid. The procedure is repeated until all the (VE dir , VE ind ) grid is explored. At the end, a surface like the one shown in Figure 1C of the main text is produced.
Technical notes: to run the SMC algorithm described above, we run multiple python scripts on computer clusters; instead of the likelihood, log-likelihood were almost always actually used, in order to transform products in sums and avoid computational underflow problems.

S3.4 Maximum Likelihood and confidence intervals
Once we generated the likelihood surface at a desired time point, it is straightforward to get the best estimates VE best dir and VE best dir from the maximum of the likelihood itself. To calculate the 95% confidence intervals (CIs) around the best estimates, it is convenient to work with the logarithm of likelihood l = log((L)). The loglikelihood l is first sliced by fixing one of the VEs on the value that maximise the l, then using the Wilks theorem [23], the 95% confidence interval can be calculated to be approximately included in the region for which: where h is a chi-squared-distributed variable with one degree of freedom, so that h = 0.381 with a confidence level of 0.05.
However, before applying the Wilks theorem, the l values must be smoothed, since they are affected by Monte Carlo variability [22]. Here this is accomplished by applying a Savitzky-Golay filter. Finally, since l depends on discrete values of VE, it is linearly interpolated to be continuous in all the VE values between 0 and 1.
Examples of sliced likelihoods produced with this method are shown in Figure 2B of the main text, where the vertical lines point out best estimate and the 95% CIs.

S4 The mathematical model realistically reproduces demography and meningococcal epidemiology in England
We developed a mathematical model to reproduce two different vaccination campaigns in England, against serogroup C and B meningococcal disease. In both cases, here we test the model running its deterministic version without vaccination campaigns for 50 years, to see if it reasonably reproduce empirical evidences. The model reproduces a realistic demographic evolution of the population in England, with fluctuations due to population differences across age groups. Figure S6 shows the evolution of the demography for the MenC setting.
The carriage prevalence remains constant for all age groups for both MenB and MenC (MenC carriage shown in Figure S7), indicating that the model has been well calibrated under the assumption of temporal stability of the asymptomatic carriage prevalence.
We calculated the serogroup-specific transmission given contact rates β a using equations S39 and S40, that basically depend on the following key quantities: the serogroup carriage prevalence by age and the duration of carriage. Direct estimates of transmission rates are never been published in literature, except for a recent work that evaluated a single carriage acquisition rate, for MenB only, in individuals of age between 10 and 25 years, finding 2.8 per 1000 person-month [24]. We tested whether our own acquisition rates are reasonably similar to the above clinical study's findings. To produce comparable numbers, we run the MenC and MenB models in absence of vaccination and then we simulated the statistical procedure adopted by the authors: we divided the overall number of transmissions J of the interested age classes by the starting susceptible population in the same age classes. We found the following acquisition rates: 1.32 per 1000 person-month in the MenC setting and 6.93 per 1000 person-moths in the MenB setting. Figure S8 shows the distribution of MenC disease cases predicted without implementing any vaccination campaign. The cases have been numerically generated combining the disease equations S33 with the deterministic version of the transmission model, defined by equations from S13 to S18. The cases occurrences by age remain roughly constant, with fluctuations only due to the demographic shifts. An analogous behaviour is observed in the MenB setting (not shown). The disease endemicity is a direct consequence of the assumption that the risk of disease given transmission remains constant over time. For the MenC setting, it is also possible to check if the mathematical model correctly reproduces meningococcal epidemiology once the immunization campaign is implemented. In Figure 3 of the main manuscript, panel A, we showed that when the VE parameters are set to the MCML best estimates, the model well reproduces the overall number of reported cases per epidemiological year. Here we show that the same argument is valid looking at different age classes. Figure S9 shows 6 panels, each one depicting time series of cases reported in England (grey bars) for 6 different age classes. The red points are the synthetic disease cases generated by our model set on our best estimates of VE dir and VE ind . In general, our simulation are in good agreement, demonstrating that the vaccination has been correctly modelled. However, the age group 20-29 yr shows a discrepancy: the real cases grew while instead of declining as predicted by our model. We note that such age class was not targeted by the catch-up campaign in 2000/01. We hypothesize that some large environmental fluctuations, not reproducible by our model, brought to the increase of disease cases. Anyway, the next year the cases drastically fell, probably mainly for the herd immunity effect. It is possible that the mixing assumptions of our model are too simplistic or not very accurate. Maybe that individuals in the 20-29 age class have more contacts within themselves and less with the remaining population, so retarding the herd immunity effect. Further studies on human contacts and their impact on the transmission of N. meningitidis would be precious to improve the understanding and modelling of meningococcal disease.

S5 Sensitivity analysis
The complexity of the inferential method and of the underlying mathematical meningococcal model required to test the sensitivity of the VE final estimates. We varied the expected number of disease cases, the duration of carriage and the duration of immunity. We tested these three key quantities because they are not accurately known and because all other quantities that are affected by uncertainties depend on them. Finally we discuss the temporal evolution of the final VE estimates. The sensitivity analyses were carried out in the MenC setting.

Expected number of disease cases
The number of disease cases by age can be subjected to unpredictable trends. We parametrized the risk of disease given infection (equation S41) using the average age-pattern of cases reported in England from January 1998 and June 1999, multiplied by total number of cases reported in England in epidemiological year 1998/99.
We tested whether different choices of the total number of cases, namely D tot , would affect the MCML inference. In particular we set D tot to be: • England cases of known age in 98/99 (base case) = 869; • England cases in 98/99 = 883; • Yearly average England cases in 97/98 1 and 98/99 = 807; • Yearly average England cases from 96/97 1 to 98/99 = 776; • Yearly average England cases from 95/96 1 to 98/99 = 728; • Low incidence scenario = 600; • High incidence scenario = 1000. Table S4 shows the final estimates of VE dir and VE ind and after two years of vaccination campaign, when D tot is changed. The direct effectiveness is almost perfectly unchanged. The indirect effectiveness gets higher as the number of expected cases is increased.     Figure 2 of the main text is shown here at the intersection of the second row from top (2 years of data) and second column (6 months of duration of carriage). The MCML method is in general sensitive to the precise value assigned to the duration of carriage τ car . However, only the indirect vaccine effectiveness is affected by changes to τ car : moving from 4 to 8 months, the 95% confidence region shifts only horizontally. While VE dir is fairly independent on τ car , VE ind is clearly affected. In the 4 months scenario, VE ind converges to less than 50%, while if τ car equals 8 months, VE ind almost reaches 100%. When τ car is set to 6 months, VE ind converges to about 70%, a value very near to the ones calculated from carriage studies [25,10], as described in main text.

Duration of immunity
We also tested the sensitivity of our results to the duration of immunity, τ imm , whose value is reported in the literature with significant uncertainty. As described in the main text, in the MenC model we have set three different values, depending from the age of the vaccinated: 15 months for routinely vaccinated infants, 5 years for infants of up to 12 months of age enrolled for the catch-up campaign, 10 years for all the other subjects included in the catch-up campaign. Varying the duration of immunity induced in subjects vaccinated in the catch-up campaign does not change the results. This is a trivial consequence of the fact that we do not use disease case notifications after 5-6 years from the beginning of the campaign.
Regarding the routine vaccination, changing the duration of immunity slightly affects our results. Figure S11 shows likelihoods of VE dir and VE ind for the MenC routine vaccination programme, calculated at different time points, by employing three different duration of immunity: 9 months, 15 months (base case [20,21]) and 24 months. Differences in the best estimates are small. However, the 24 months scenario displays a higher value of the likelihood maximum, indicating that previous estimates could underestimate the duration of immunity. We think that the immunity persistence should be further investigate to improve the modelling of this process, that in other scenarios could significantly impact the VE estimation and Public Health decisions.

Temporal evolution of the VE estimates
In addition, Figures S10 and S11 also show how the estimates evolve in time (epidemiological years, overall between mid 2000 and mid 2005). In general, the 95% confidence regions shrink as more data are used. When IMD cases notified during the first epidemiological year only are employed, the best estimates of VE ind are significantly different if compared to the best estimates obtained using more data, even if the precision after one year is very low, with a very broad 95% CI of the indirect effectiveness. On the other hand, the assessment of VE dir is extremely fast to converge and does not show systematic deviations, indicating that the direct effectiveness of MenC conjugate vaccines in England could have been estimated even after one year only.