Data

Published phase III clinical trial data are available for 12 months¹⁷, 18 months¹⁸, 24 months¹⁹, 36 months²⁰ and 54 months²¹ after the second dose. From these publications we extracted the number of symptomatic VCD cases (\({N}_{{{\rm{symp}}}}\)) and the number of VCD cases leading to hospitalization (\({N}_{{{\rm{hosp}}}}\)) at the finest stratification available, in Excel (version 16.75).

Per-protocol outcomes and population sizes for the Qdenga phase III clinical trial are shown in Supplementary Fig. 2. \({N}_{{{\rm{symp}}}}\) was published by the trial arm \(v\) (1 for placebo, 2 for vaccinated), classified baseline serostatus \(b\) (−, seronegative; +, seropositive), and infecting serotype \(k\) (1, 2, 3, 4 for DENV1–4, respectively), within each reporting interval \(d\) (1, 1–12 months; 2, 13–18 months; 3, 19–24 months; 4, 25–36 months; 5, 37–48 months; 6, 49–54 months). The number of hospitalized cases, \({N}_{{{\rm{hosp}}}}\), was published by trial arm, baseline serostatus and infecting serotype at 1–24 months (safety set, individuals who received at least one dose of the vaccine or placebo), as well as for the last three reporting intervals. \({N}_{{symp}}\) and \({N}_{{{\rm{hosp}}}}\) were both published by trial arm, baseline serostatus and age group \(j\) (1, 4–5 years; 2, 6–11 years; 3, 12–16 years), for the first four reporting intervals \(d\) (that is, up to 36 months), and by age group and infecting serotype for reporting intervals 1–12 months and 13–24 months. Limited country-specific data were published, therefore all data and model estimates presented here are for the whole trial population across all countries.

The trial populations by baseline serostatus and trial arm were reported for each time interval, \({{N}_{{{\rm{pop}}}}}_{{bv}}(d)\), with the population further broken down by age group \(j\) up to 36 months (Supplementary Fig. 2g,h). Given that only the first VCD case is reported in the trial, we right-censor cases at the end of each reporting interval \(d\) and reduce the monitored population in the next time interval accordingly (that is, \({{N}_{\rm{{{surv}}}}}_{{bv}}(j,d)={N}_{\rm{{{po}{p}}}_{{bv}}}(j,d)-{N}_{\rm{{{{symp}}}}_{{bv}}}(j,d-1)\)). Age-specific data were not published after month 36 of the trial, therefore we estimate the age-specific trial populations at 37–48 months and 49–54 months, assuming the same age-specific probability of loss to follow-up as the first 36 months of the trial.

Cohort survival model of vaccine efficiency

To estimate VE we calibrated a Bayesian cohort survival model to the above data (Supplementary Fig. 1a). We assumed that an individual may be infected up to four times, and any infection could be symptomatic. We did not track further infections among individuals who had a detected infection. The model was solved at monthly time intervals (t = 1–54).

Fitting the initial conditions

Let \({p}_{{kj}}\) denote the probability of exposure to serotype \(k(=\mathrm{1,2,3,4})\) in age group \(j\) prior to the start of the trial \(({t=T}_{0})\). We estimate \({p}_{k3}\) (the probability of exposure in 12–16 years). For the other two age groups, we estimate the cumulative probabilities of exposure to any serotype and scale serotype-specific probabilities to match the distribution of serotype-specific exposure in age group 3. Let \({h}_{k3}\) be the serotype-specific hazard of exposure in age group 3:

The serotype-specific hazards in age groups \(j=1\) and \(j=2\) are therefore given by:

Finally, the serotype-specific probabilities of exposure in age groups \(j=1\) and \(j=2\) are:

To aid identifiability of both the FOI and probability of symptoms, we constrain the probability of past exposure \(p\) by the estimated FOI during the trial. Let \({\lambda }_{m}\) be the average monthly FOI during the trial (irrespective of serotype). The cumulative probability of exposure in the oldest age group \({p}_{m}\) is therefore:

where 14 is the mean age in the oldest age group and the factor of 12 reflects the time unit of months used in our analysis. We therefore set a beta prior for \({p}_{k3}\) with a mean \({p}_{m}\) and variance \(\mathrm{var}=0.02* {p}_{m}(1-{p}_{m})\). A value of 0.02 was chosen to ensure that the prior was informative. We reparametrized the beta distribution to obtain standard shape parameters \({\rm{{shape}}}1\) and \({\rm{{shape}}}2\).

Let \(c\) denote the number of serotypes to which an individual had been exposed (0, 1, 2, 3, 4). The probability that an individual has been exposed to serotypes \(\Omega\) within set \({h}_{c}\), defined as \({h}_{0}=\{\varnothing \}\), \({h}_{1}=\left\{\left\{1\right\},\left\{2\right\},\left\{3\right\},\{4\}\right\}\), \({h}_{2}=\left\{\left\{\mathrm{1,2}\right\},\left\{\mathrm{1,3}\right\},\left\{\mathrm{1,4}\right\},\left\{\mathrm{2,3}\right\},\left\{\mathrm{2,4}\right\},\{\mathrm{3,4}\}\right\}\), \({h}_{3}=\left\{\left\{\mathrm{1,2,3}\right\},\left\{\mathrm{1,2,4}\right\},\right.\)\(\left.\left\{\mathrm{1,3,4}\right\},\{\mathrm{2,3,4}\}\right\}\) and \({h}_{4}=\left\{\{\mathrm{1,2,3,4}\}\right\}\) is given by:

Accounting for the imperfect test sensitivity, \({{\rm{sens}}}\), and specificity, \({{\rm{spec}}}\), of the microneutralization test used to classify individuals as seronegative (−) or seropositive (+) at baseline⁴⁹, the probability of being classified as seropositive at baseline is:

We estimate \({p{\rm{SP}}}\left(j{,T}_{0}\right)\) by fitting the number of baseline seropositive individuals by age group \(j\) across both trial arms, \({N}_{{{\rm{SP}}}}\left(j,{T}_{0}\right)\) to the observed data, assuming a binomial distribution, and the observed number of participants \({N}_{{{\rm{pop}}}}\left(j,{T}_{0}\right)\) in age group \(j\) over both the vaccine and placebo arms:

Estimating vaccine efficacy

We used a modified form of the correlates of protection model first proposed by Khoury et al.²⁷ (to model VE for SARS-CoV-2 vaccines) to link the imputed neutralizing antibody titer \({n}_{{ck}}\left(t\right)\) against serotype \(k\) (Supplementary Fig. 3) with the risk ratio (comparing vaccinated and unvaccinated individuals) of symptomatic disease \({{{\rm{RR}}}_{{{\rm{symp}}}}}_{{cvkj}}\left(t\right)\):

By definition, \({{{\rm{RR}}}_{{{\rm{symp}}}}}_{{cvkj}}=1\) if \(v=1\). Here, \({L}_{{ck}}\) is the maximum potential vaccine-induced disease enhancement with serotype \(k\) associated with vaccination. This is assumed to be zero for all individuals with one or more previous DENV exposures (\(1\le c\le 4\)), such that the maximum \({{\rm{RR}}}_{{symp}}\) is 1 for all baseline seropositive vaccinees. We estimate \({L}_{{ck}}\) (≥0) for vaccinated individuals with no prior DENV exposure (\(c=0\)), such that their maximum \({{\rm{RR}}}_{{{\rm{symp}}}}\) may be >1. The parameter \({n}_{{50}_{{ck}3}}\) represents the neutralizing titer conferring 50% protection against disease (in the absence of enhancement) for the oldest age group \((j=3)\), which is scaled by \({e}^{\;{{\beta }_{{{\rm{symp}}}}}_{j}}\) for individuals in the younger age groups, \(j\in \{\mathrm{1,2}\}\). Finally, \({w}_{{ck}}\) controls the steepness of the relationship between neutralizing titer and the risk ratio.

Similarly, the risk of being hospitalized with dengue for a participant in the vaccination arm of the trial relative to one in the control arm is:

This is the same model as used for the symptomatic VCD endpoint in equation (8) above with two modifications: first, the maximum level of vaccine-induced enhancement for symptomatic VCD, \({L}_{{ck}}\), is scaled by the multiplier \({\tau }_{k}\); and second, the neutralizing titer that gives 50% protection from hospitalization, \({n}_{50}\) is assumed to be proportional to that for VCD, but scaled by the factor \({e}^{-{\alpha }_{{ck}}}\). In addition, the variation of efficacy with age (represented by the parameter \({{\beta }_{{{\rm{hosp}}}}}_{j}\)) is allowed to be different from that for the symptomatic VCD endpoint.

Note that because \({{\rm{RR}}}_{{{\rm{symp}}}}\) and \({{\rm{RR}}}_{{{\rm{hosp}}}}\) are dependent on the cumulative number of DENV exposures at month \(t\), not the number of exposures prior to vaccination, we are assuming that the order of vaccination and infection does not matter. Therefore, for example, a baseline monotypic vaccinee escaping infection up to month \(t\) has the same disease risk as a baseline seronegative vaccinee with a single breakthrough infection prior to month \(t\).

We assume that the vaccine-induced neutralizing antibody titers are characterized by an initial period of fast decay, with half-life \({{hs}}_{b}\) (decay rate \(\pi {1}_{b}=-{\mathrm{ln}}\left(2\right)/{{hs}}_{b}\)), followed by a period of slow decay with half-life \({hl}\) (decay rate \(\pi 2=-{\mathrm{ln}}\left(2\right)/{hl}\)):

Here, \({n}_{{0}_{{bk}}}\) is the fitted initial neutralizing antibody titer to serotype \(k\), \(b\) is the baseline serostatus prior to vaccination (\(b\in \{-,+\}\)), \(t{s}_{b}\) is the time at which decay switches from fast to slow, and \(t\) is the month elapsed since the second dose.

Likelihood

Within each trial reporting period \(d\), we assume a constant serotype-specific FOI over time, \({\lambda }_{k}\left(d\right)\). Given that the model is solved at monthly intervals \(t\), we set \({\lambda }_{k}\left(t\right)={\lambda }_{k}\left(d\right)\) for all months \(t\) within the reporting period \(d\).

For each month \(t\), we calculate the probability of surviving that month without infection, given as \({e}^{-{\sum }_{k\notin \Omega }{\lambda }_{k}\left(t\right)}\) (Supplementary Fig. 1a). The probability of being seronegative at the start of month \(t+1\) is therefore:

The probability of having been infected by the set of serotypes \(\Omega\) by month \(t+1\) can be represented as the sum of those who were already exposed to the serotype combination in set \(\Omega\) and who escaped infection by all serotypes not in \(\Omega\) during month \(t\), plus those with previous exposure to all serotypes in set \(\Omega\) excluding serotype \(k\) (denoted \({\Omega }_{{\rm{\backslash k}}}\)) who had an undetected infection to serotype \(k\) during month \(t-12\) (so that individuals may not be re-infected within the 12 month period of heterotypic cross-immunity, and censoring individuals with detected infections following the date of detection^20,21):

where \((1-{{\rm{e}}}^{-{\lambda }_{k}\left(t-12\right)}){{pE}}_{{\Omega }_{{\rm{\backslash k}}}{bv}}\left(j,t-12\right)\) is the incidence of infection with serotype \(k\) during month \(t-12\) in individuals previously exposed to all serotypes in set \(\Omega\) excluding \(k\). The probability of developing symptoms upon infection is denoted \({p}_{\rm{{{sym}{p}}}_{{ck}}}\). The second term in equation (12) represents the incidence of asymptomatic infection in the 12 months prior. Let \(\gamma\) = the probability of symptoms during a secondary infection, \({\rho }_{k}\) = the serotype-specific risk of symptoms during a primary infection relative to a secondary infection, and \(\varphi\) = the risk of symptomatic disease during a post-secondary infection compared with a primary infection. Thus, \({p}_{{{\rm{sym}{p}}}_{{ck}}}={\rho }_{k}\gamma \,\) for primary infections (for \(c=0\)), it equals γ for secondary infections (for \(c=1\)), and it equals \({\rho }_{k}\gamma \varphi\) for post-secondary infections (for \(2\le c\le 3\)).

The incidence of primary infection is:

The incidence of secondary infection is:

where \({{pE}}_{{\Omega }_{{\rm{\backslash k}}}{bv}}\left(j,t\right)\) is the probability of having prior exposure to serotypes in set \(\Omega \in {h}_{1}=\left\{\left\{1\right\},\left\{2\right\},\left\{3\right\},\left\{4\right\}\right\}\), not including serotype \(k\).

The incidence of tertiary infection is:

where \(\Omega \in {h}_{2}=\left\{\left\{1,2\right\},\left\{1,3\right\},\left\{1,4\right\},\left\{2,3\right\},\left\{2,4\right\},\{3,4\}\right\}\).

The incidence of quaternary infection is:

where \(\Omega \in {h}_{3}=\left\{\left\{1,2,3\right\},\left\{1,2,4\right\},\left\{1,3,4\right\},\{2,3,4\}\right\}\)

The incidence of disease and hospitalization, respectively, are:

Here, \({{p}_{{{\rm{hosp}}}}}_{{ck}}\) is the probability of hospitalization. Let \({\delta }_{k}=\) the probability that a symptomatic case due to serotype \(k\) is hospitalized and \(\epsilon\) = the risk of hospitalization in secondary cases, compared with primary or post-secondary infections. Thus, \({{p}_{{{\rm{hosp}}}}}_{{ck}}={\delta }_{k}\epsilon\) when \(c=1\), and \({\delta }_{k}\) otherwise.

Summing over the serotype combinations, the total symptomatic case and hospitalization incidence are:

The monthly symptomatic and hospitalization incidences are aggregated to match the months in each trial reporting period \(d\), and the expected number of symptomatic and hospitalized cases are calculated as:

where \({{N}_{{{\rm{surv}}}}}_{v}\left(j,d\right)\) is the observed population size in reporting period \(d\), accounting for censoring of symptomatic cases in the previous time intervals. The observed total number of symptomatic cases and hospitalizations at time interval \(d\) are assumed to follow binomial distributions:

where \({pS}\left(d\right)=\sum _{b}\sum _{v}\sum _{k}\sum _{j}{S}_{{bvk}}(j,d)/{N}_{{{\rm{surv}}}}\left(d\right)\) and \({pH}\left(d\right)=\sum _{b}\sum _{v}\sum _{k}\sum _{j}\)\({H}_{{bvk}}(j,d)/{N}_{{{\rm{surv}}}}\left(d\right)\).

For each published time interval d, the expected distribution of symptomatic cases (\({{pS}}_{{bvk}}\left(d\right)\)) or hospitalizations (\({{pH}}_{{bvk}}\left(d\right)\)) in baseline serostatus \(b\), due to serotype \(k\), in trial arm \(v\), relative to the total number of symptomatic cases (\({{\rm{cum}}S}\left(d\right)=\sum _{b}\sum _{v}\sum _{k}\sum _{j}{S}_{{bvk}}(j,d)\)) or hospitalizations (\({{\rm{cum}}H}\left(d\right)=\sum _{b}\sum _{v}\sum _{k}\sum _{j}{H}_{{bvk}}(j,d)\)), is given by:

Equally, \({{pS}}_{{bv}}\left(j,d\right)\) and \({{pH}}_{{bv}}\left(j,d\right)\), the expected proportion of symptomatic cases and hospitalizations, respectively, in baseline serostatus \(b\), trial arm \(v\), age group \(j\) are given by:

Finally, \({{pS}}_{k}\left(j,d\right)\) and \({{pH}}_{k}\left(j,d\right)\), the expected proportion of symptomatic cases or hospitalizations, respectively, due to serotype \(k\) in age group \(j\) are given by:

The expected probabilities of observing symptomatic disease and hospitalized cases given in equations (25)–(30) are assumed to follow multinomial distributions:

Model variants and sensitivity analyses

To explore whether the data supported varying parameters by baseline serostatus, serotype, outcome and age group, we considered multiple simpler nested models than those outlined above. Model variants were compared visually and quantitatively using the log-likelihood and the expected log-predictive density, estimated using the Watanabe–Akaike information criterion, to find the most parsimonious, biologically motivated model that could reproduce the trial data.

We ran a sensitivity analysis on the choice of the prior normal distribution of the enhancement parameter \({L}_{{ck}}\) centered and truncated at 0, with standard deviations {0.25, 0.5, 0.75, 1.00, 1.25, 1.50, 1.75, 2.00}. We assessed the impact of the choice of the prior distribution on the log-likelihood, posterior estimate of \({L}_{{ck}}\), and the Bayes factor estimate comparing models with and without vaccine-associated enhancement (that is, setting \({L}_{{ck}}=0\)).

We also explored the sensitivity of the VE estimates to assumptions made about: the period of heterotypic cross-immunity, which we fixed at 1 and 6 months (rather than 12 months as in our main analysis); whether multitypic individuals could be misclassified as seronegative at baseline; whether all post-secondary infections are asymptomatic; whether the titer required for 50% protection against symptomatic dengue (\({n}_{50}\)) in seronegative and multitypic individuals should be serotype specific; and whether the risk ratio of vaccine-associated hospitalization enhancement compared with symptomatic disease (\(\tau\)) should be serotype specific.

Simulation study to assess the model identifiability

To confirm the identifiability of the model parameters given the data, we conducted a simulation study. We sampled 20 sets of the parameters to generate simulated datasets from distributions that were wider or centered away from the prior distributions used for parameter inference. The only exception was the parameters describing the titer decay (\({hs}\), \({hl}\) and \({ts}\)), which were estimated directly from the observed data and were intentionally kept informative. We then simulated case data using each set of parameters. The case data were comparable to the observed data in magnitude and distribution across time, age group, serostatus and trial arm. Finally, we calibrated the survival model to the simulated case data and compared the simulated posterior distribution of each parameter to the true parameter distributions.

Inferential framework

The model was fitted to the baseline seropositivity data (\({N}_{{{\rm{SP}}}}\)) and number of VCD symptomatic (\({N}_{{{\rm{symp}}}}\)) and hospitalized (\({N}_{{{\rm{hosp}}}}\)) cases observed across the trial, stratified by baseline serostatus, age and serotype within a Bayesian framework. Parameter inference was carried out using the Hamiltonian Monte Carlo algorithm and the No-U-Turn sampler via Stan, a probabilistic programming language implemented via the cmdstanr package⁵⁰ (version 0.8.1) in R⁵¹ (version 4.2.2) and RStudio (version 2023.06.1+524). Four chains were run for 10,000 iterations each, and we discarded the first 5,000 iterations as burn-in. Convergence was assessed visually and using the \(\hat{R}\) statistic⁵². The estimated parameters, priors and prior sources are listed in Supplementary Table 3.

Transmission model

To evaluate the population-level impact of Qdenga, we applied the VE estimates obtained above to a 4-serotype stochastic compartmental model of dengue transmission, a deterministic version of which was published by Ferguson et al.⁴ (Supplementary Fig. 1b). The model includes mosquito population dynamics with a seasonally varying carrying capacity. It models primary to quaternary human infections, assuming that homotypic immunity is lifelong, and heterotypic immunity is temporary (12 months). It tracks the serotype-specific infection history but not the order of past infections. The model is stratified by vaccine status, serostatus and age, with single age classes for the first 20 years and 10 years thereafter. The model incorporates realistic non-stationary human demography, calibrated to the demographic estimates for Brazil and the Philippines published in the UN World Population Prospects 2022 (ref. ⁵³), which show that the Brazilian population is older than the Philippine’s population (Supplementary Fig. 29). For this study, we embedded the VE model detailed above into this transmission model. Full details of the transmission model are given below.

Mosquito dynamics

The transmission model accounts for mosquito population dynamics following the Ross–MacDonald model, with varying seasonal carrying capacity. Specifically, adult female mosquitoes lay eggs at rate \(\Gamma\) and larvae (\(L\)) mature into adult mosquitos (\(M\)) at rate \({\rm{{\rm E}}}\). Larvae are regulated by a power density-dependent mortality rate \(\mu\), with a seasonal carrying capacity \(K\), ensuring that mortality increases as the larval population increases. Adult mosquitoes are born susceptible (\({M}_{S}\)), experience a serotype-specific FOI \({\Lambda }_{k}\), and infected mosquitoes enter an exposed but not infectious compartment (\({M}_{E}\)) that lasts on average for a period \(\frac{1}{{\boldsymbol{\xi }}}\) days, the extrinsic incubation period. Finally, the adult mosquitoes that survive the extrinsic incubation period become infectious to humans (\({M}_{I}\)) until their death. Adult mosquitoes die at a rate \(\Delta\). The ordinary differential equations (ODEs) governing the mosquito dynamics are given as follows:

The mortality rate \(\mu\) of the larvae is defined as:

where \(\sigma\) is the low-density limit of larval death and \(N\) is the total human population.

In turn, the seasonal carrying capacity \(K\) is defined as:

Where \(\bar{K}\) is the mean carrying capacity, \({K}_{s}\) is the magnitude of the seasonal variation in carrying capacity, and \(t\) is time in days. We ran model simulations across nine transmission intensities, characterized by SP9 (see the model simulation section below for more details). For each value of SP9 between 10% and 90% in 10% steps, we determined the mean carrying capacity that led to that value of SP9 across the period 2020–24. The resulting carrying capacity values gave adult female mosquito densities per person (\(A\)) ranging from 1 to 5, compatible with values reported by Focks et al., who estimated a maximum of 1.5 Aedes aegypti pupae per person for ambient air temperatures of 28 °C (ref. ⁵⁴):

Values and data sources for mosquito parameters are provided in Supplementary Table 4.

The rate at which adult female mosquitoes lay eggs \(\Gamma\) is assigned to give the required mosquito population reproduction number (that is, representing the reproductive capacity of the mosquito population, rather than anything to do with disease transmission):

The probability of transmission from human to mosquito \({{\rm{{\rm B}}}}_{{hm}}\) is defined to give the required dengue \({R}_{0}\) as follows:

Where \(\eta\) is the human infectious period and \(\kappa\) is the biting rate per mosquito.

Finally, the FOI experienced by mosquitoes due to dengue serotype \(k\), \({\Lambda }_{k}\), is:

Where \(c{I}_{k}\) is the human incidence of infection (defined below in equation (50)) with serotype \(k\), and \(N\) is the total size of the human population.

Human transmission dynamics

As with the cohort survival model, let subscript \(\Omega\) denote the set of serotypes to which an individual has previously been exposed (with cardinality \(c\)). In the transmission model, we assume that individuals are born susceptible to all four serotypes \({S}_{\varnothing }\), after which they are subject to time-varying serotype-specific FOI \({\lambda }_{k}(t)\). Infection with any of the four serotypes confers long-lasting homotypic immunity and heterotypic immunity against the other serotypes for a period \(1/{\rm{\vartheta }}\) (12 months) (\({R}_{\Omega }\)). After this, individuals are assumed to be susceptible to infection with heterotypic serotypes. In total, an individual may be infected up to four times, and we track an individual’s infection history but not the serotype-specific order of infection. The model is solved daily (\(t\)).

The model is additionally stratified by the age group i (1–20 years of single age classes, then 10 year age classes up to 80 years), and vaccination status \(v\) (1, unvaccinated; 2, vaccinated), with individuals aged \(f\) vaccinated with coverage \({\rm{\upsilon }}\). Aging is modeled as an annual discrete event (for example, all 1-year-olds move to the 2-year-old age group) to track the vaccinated cohort accurately up to the age of 20. We assume that vaccination also occurs once per year at the same time as the ageing process and we do not separately model the receiving of the first and second dose. The birth rate \(\Psi ({\rm{t}})\) and the age-specific human mortality rate \({{\rm{\theta }}}_{{\rm{i}}}(t)\) are both calibrated to match the demography of either the Philippines or Brazil (see the model simulation section below for additional details). The ODEs governing human transmission dynamics are:

Here \({{\rm{\delta }}}_{{\rm{x}},{\rm{y}}}\) is the Kronecker delta function (equal to 1 when \({\rm{x}}={\rm{y}}\) and 0 otherwise), the first term of \(\frac{d{S}_{\Omega {vi}}}{{dt}}\) represents births, the second term represents vaccination, and the final term represents infections and deaths. The parameter \({\rm{RR}}_{{\inf}_{cvki}}\left(t\right)\) is the vaccine-associated RR of infection (see equation (51) below for details). The first term of \(\frac{d{R}_{\Omega {vi}}}{{dt}}\) represents vaccination, the second term represents infection, and the third term represents waning heterotypic immunity and mortality. All transitions between compartments (ageing, infection, waning immunity, death and vaccination) are drawn from binomial distributions to introduce model stochasticity.

The FOI on humans due to serotype \(k\) is:

where \({{\rm{{\rm B}}}}_{{mh}}\) is the per bite transmission probability from mosquitoes to humans.

The incidence of individuals exposed (but not yet infectious) to serotype \(k\) is given as:

The incidence of individuals who are infectious \(I\), symptomatic \(D\), or hospitalized \(H\) due to serotype \(k\) at time \(t\) is given as:

Here \(1/\zeta\) is the human incubation period, \({p}_{{sym}{p}_{c}}\) is the probability that an infected individual with serostatus \(c\) is symptomatic (= \(\rho \gamma\) for primary infections, \(\gamma\) for secondary infections and \(\rho \gamma \varphi\) for post-secondary infections, see equation (12) above for parameter descriptions), \(Q\) is the probability that a symptomatic individual requires hospitalization (fixed at 9% to match the average probability of hospitalization observed in Brazil and the Philippines during the phase III clinical trial), and \({{{\rm{RR}}}_{{symp}}}_{{cvki}}\left(t\right)\) and \({{{\rm{RR}}}_{{{\rm{hosp}}}}}_{{cvki}}\left(t\right)\) are the vaccine-associated RR of disease and hospitalization (defined in equations (8) and (9) above).

The cumulative incidence of individuals infectious to serotype \(k\), irrespective of serostatus, vaccine status or age group, as seen in equation (43) above is therefore:

Model parameter values

The VE parameters (\({hs}\), \({hl}\), \({ts}\), \(L\), \({n}_{50}\), \(\alpha\), \(\beta\), \(w\) and \(\tau\)) and probabilities of symptoms (\(\rho\), \(\gamma\), \(\varphi\)) were sampled from the posterior distribution of the survival model used to reconstruct \({\rm{RR}}_{\rm{symp}}\) and \({{{\rm{RR}}}}_{{{\rm{hosp}}}}\) as described in equations (8) and (9). All other parameter values used in the stochastic compartmental model of transmission are listed in Supplementary Table 4.

Model simulations

We equilibrated the transmission dynamics by running the transmission model in the absence of vaccination for 175 years, starting from 1850. For the first 100 years we assumed a static 1950s demography. From 1950 onwards we assumed time-varying demographies matching the UN World Population Prospects 2022 estimates for Brazil or the Philippines⁵³, to explore the sensitivity of the impact estimates to different population structures (Supplementary Figs. 29 and 30). We assumed that vaccination started in 2024.

We ran multiple scenarios to investigate the vaccines impact across: two demographies (Brazil and the Philippines); nine transmission intensity settings, defined using SP9, from 10% to 90% in steps of 10%; four vaccine coverages (20%, 40%, 60%, 80%); seven ages at vaccination (from 6 to 12 years of age); four hypotheses about the vaccine’s mechanism of action, obtained by assuming vaccine waning for either 5 years (approximately the period of the phase III trial) or 15 years post-vaccination (extrapolating beyond the trial), and that the vaccine protects against clinical disease only (VS), or also against infection (VI, see equation (51) below); and finally, with and without pre-vaccination screening of individuals’ serostatus assuming a diagnostic test with 94.7% specificity and 89.6% sensitivity²⁹.

We expressed transmission intensity in terms of SP9 because the latter can be directly measured from seroprevalence surveys. Conversely, the reproduction number (R₀) depends on additional assumptions (for example, the extent to which post-primary infections contribute to transmission, as described in the Imai et al. study⁵⁵) and is therefore difficult to estimate empirically for dengue. However, we note that specifying one of R₀, FOI and SP9 determines the values of the other two variables, as shown in Supplementary Fig. 30. Supplementary Fig. 31 shows model assumptions about the seasonality of dengue transmission intensity. We assumed that all serotypes have near-identical transmissibility but relaxed this assumption in a sensitivity analysis.

In the absence of serological tests able to identify the infecting serotype and infer asymptomatic infections, data on Qdenga’s efficacy against infection are limited, but there is evidence that it may offer limited protection against infection for some months²⁸. In the VI modeling scenarios we therefore assumed that the probability of infection, as for clinical outcome, can be explained by the antibody titer induced by previous infection⁸ or vaccination¹³. We assumed that the titers required for protection from infection are higher than those required to protect against symptomatic disease, and that (unlike enhancement of disease) there is no enhancement of infection risk^8,13. We therefore modeled a scenario of low–moderate VE against infection by scaling the posterior estimates for \({n}_{{50}_{{ck}3}}\)(the neutralizing titer conferring 50% protection against disease in the absence of infection) in equation (8) by \({{\alpha }_{\inf }}_{c}\), equal to 12 if \(c=0\) (primary infection) and 3 if \(1\le c\le 3\):

The resulting \({{\rm{VE}}}_{\inf }\) curves are plotted in Supplementary Fig. 14.

For each of the modeling scenarios, we sampled the VE parameters from the posterior distributions 200 times (Supplementary Table 1). Then, for each posterior parameter set we ran the stochastic transmission model 50 times (where the probability of infection and vaccination are drawn from binomial distributions), giving 10,000 simulations in total for each scenario (Supplementary Fig. 1b).

Each model simulation was run twice: once with VE set to our estimated values (‘vaccination’ scenario) and once with VE set to zero (‘no vaccination’ scenario). This enabled us to estimate the population- and individual-level impact of vaccination. Individual impact is measured as the proportion of cases averted in the first vaccinated cohort (over a period of 10 years) and population-level impact is the proportion of cases averted in the entire population, including non-vaccinated individuals (over a period of 10 years) (Supplementary Fig. 1b):

Here \({\rm{cases}}_{\rm{NV}}\) is the number of cases in the zero VE (‘no vaccination’) counterfactual scenario, and \({\rm{cases}}_{\rm{V}}\) is the number of cases in the VE > 0 (‘vaccination’) scenario.

We summarize uncertainty in two ways. The first is the overall uncertainty, calculated over the 10,000 realizations, which represents uncertainty in both the posterior parameter estimates and the stochastic simulations. The second is the parameter uncertainty, calculated as the uncertainty of the mean of the 50 stochastic simulations for each of the 200 posterior samples.

We note that the vaccination and counterfactual (that is, no vaccination) runs were matched exactly (that is, identical transmission dynamics were simulated for each such pair of simulations) up to the point of vaccine introduction. In the VS scenario (no VE against infection), transmission dynamics remained identical for the vaccination and no-vaccination scenarios following the introduction, but this was not the case if non-zero VE against infection was assumed (given that vaccination then modified transmission).

Transmission model sensitivity analyses

We explored the sensitivity of the population impact estimates to, first, the period of heterotypic cross-immunity, which we assumed to be 6 months (rather than 12 months as in our main analysis); second, the VE estimates obtained when assuming serotype-specific titers for 50% protection against symptomatic dengue in seronegative and multitypic individuals (rather than single seronegative and multitypic estimates used in our main analysis); third, the VE estimates obtained when we estimate a single risk ratio of vaccine-associated hospitalization enhancement compared with symptomatic disease (rather than serotype specific as in our main analysis); fourth, the lack of seasonality in mosquito carrying capacity in the transmission model; and last, the R₀ differences between serotypes, with relative values of 1.05, 1.15, 0.95 and 0.85 for DENV1–4, respectively.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.