Model structure
We developed a new model of gonorrhea transmission and vaccination, building on previous work10,12 which represents England’s sexually active MSM population, stratified by sexual behavior, current infection status, history of gonorrhea diagnosis, and vaccination status/vaccine sentiment. The model accounts for rates of SHC attendance due to partner notification.
Within the model, individuals who are uninfected (U) can become infected through sexual contact with a contagious individual. Newly-infected individuals pass through an incubating state (I), before developing symptoms (S), or remaining asymptomatic (A). Infected individuals who attend an SHC (either due to care-seeking for symptoms or asymptomatic screening) are treated and enter a treatment state (T). Treated individuals are cured, returning to the U state; infection does not confer natural immunity, and previously-infected individuals are as susceptible as those never infected.
Heterogeneity in sexual behavior is represented by low (L) and high (H) sexual activity groups; high-activity MSM have more partners and attend SHCs for screening more frequently.
The population is stratified by vaccination status/sentiment. Individuals who are fully-vaccinated (Vi, Ve) or revaccinated (Ri, Re) have lower rates of acquisition than partially-vaccinated individuals (P), who have lower rates of acquisition than those who are unvaccinated (X) and those whose vaccine protection has waned (W). We also consider, in supplementary analysis, a scenario in which a proportion of individuals are unwilling to be vaccinated (H).
In order to consider vaccination strategies in which vaccination is offered to those with a clinical history of gonorrhea, the model distinguishes between individuals diagnosed in the last D years (represented by d = 1) from those with no previous diagnosis or who were last diagnosed more than D years ago (represented by d = 0). All treated individuals move into the U state in the “Diagnosed in the last D years” (d = 1) stratum. Over time, individuals in the “Diagnosed in the last D years” stratum move into the corresponding infection state in the “Not diagnosed in the last D years” (d = 0) stratum.
Eligibility for revaccination is dependent upon the time since the last vaccination. We assume that to allow for variation in the duration of protection, repeat vaccination is offered to individuals after a period shorter than the average duration of protection. The model distinguishes between recently vaccinated individuals (Vi), who are ineligible for revaccination, and vaccinated individuals who were vaccinated more than EV years ago (Ve), who are eligible for revaccination despite still being vaccine-protected. As a two-dose primary vaccination provides protection for DV years, individuals remain in Ve for a period of (DV–EV) years before their protection wanes. In the same vein, the model distinguishes between recently revaccinated individuals (Ri, ineligible) and revaccinated individuals who received their booster vaccine more than ER years ago (Re, eligible). As booster vaccination provides protection for DR years, individuals remain in Re for a period of (DR–ER) years before their protection wanes. We assume that individuals who only received one primary dose are offered another course of two-dose primary vaccination (according to the implemented vaccination strategy), and those who accept can receive either one or two further doses. Those in P who receive one dose remain in P, while those who receive two doses move to Vi.
Figure 4 shows a diagram of the model’s structure, and the meaning and value of model parameters are detailed in Table 4. We denote the number of individuals in infection state Z, activity group a, vaccination status v, and diagnosis history d at time t as \({Z}_{a}^{v,d}(t)\) (or simply \({Z}_{a}^{v,d}\) when context is clear). The model is described by Equations (1–5).
The population is divided into compartments representing different states of infection, and stratified according to diagnosis history and vaccination status. Note that there are separate sets of compartments for low and high sexual activity groups, which have identical arrangements of compartments and vaccination/diagnosis history strata. a Infection states and stratification by diagnosis history. Individuals enter the population uninfected (U); those leaving through aging leave from any state. Uninfected individuals can acquire infection, passing through an incubating state (I) before either developing symptoms (S) or remaining asymptomatic (A). Symptomatic individuals seek treatment due to symptoms, while a proportion of asymptomatic infections are identified via screening. Diagnosed individuals enter the treatment state (T). Treated individuals and untreated asymptomatic individuals who recover naturally return to U. Those diagnosed in the last D years (where D = 1 or 2) are distinguished from those never diagnosed or who were last diagnosed over D years ago. All treated individuals move into U in the “Diagnosed in the last D years” stratum. Over time, individuals in the “Diagnosed in the last D years” stratum move into the corresponding infection state in the “Not diagnosed in the last D years” stratum. Under all considered vaccine-targeting strategies, vaccination is offered to unprotected individuals diagnosed with (and treated for) gonorrhea. Under some strategies, vaccination is also offered to some uninfected individuals at screening (criteria in Table 1). Those who accept enter U in a vaccine-protected stratum. b Shows stratification according to vaccination status: unvaccinated (X), partially-vaccinated (P), fully-vaccinated (split into those who are ineligible for revaccination due to recent vaccination, Vi, and those who are eligible for revaccination, Ve), waned (W), and revaccinated (split into those who are ineligible due to recent revaccination, Ri, and those who are eligible for repeat vaccination, Re). Unvaccinated individuals in X who accept vaccination and receive one dose enter P, those who accept vaccination and receive two doses enter Vi, and those who decline vaccination remain in X. Individuals in P remain eligible for two-dose primary vaccination, those who choose to receive only one dose remain in P and those who receive two doses enter Vi. Individuals in Vi become eligible for single-dose booster vaccination after a period of EV years, moving to Ve. When vaccine protection wanes, partially vaccinated individuals return to X, while fully vaccinated individuals in Ve enter W. Individuals in Ve and W who accept single-dose booster vaccination enter Ri. Individuals in Ri become eligible for booster vaccination after a period of ER years, moving to Re. When protection from booster vaccination wanes, individuals move from Re back to W.
The terms αv (the rate of new entrants into the population entering vaccination-status stratum v), rv (the probability of accepting vaccination and moving between vaccination-status strata), βv (the per-partnership probability of gonorrhea acquisition), and \({\theta }_{a,v,d}^{Z,wane}\) (the rate of moving between vaccination-status strata due to vaccination waning) take different formulations for different vaccination-status strata.
Dependent on the vaccination-targeting strategy implemented, vaccination is offered to uninfected individuals attending for screening (who were not notified partners), to uninfected individuals attending for screening due to partner notification, or to infected individuals when they are diagnosed. The rates of movement between vaccination-status strata by these are denoted by \({\theta }_{a,v,d}^{VoS}\), \({\theta }_{a,v,d}^{VaN}\), and \({\theta }_{a,v,d}^{VoD}\) respectively, which take different formulations for different vaccination strata and depend upon the vaccination strategy implemented. Terms that differ between vaccination strata are detailed in Equations (6–12).
$$\frac{d{U}_{a}^{v,d}}{dt}=\left\{\begin{array}{l}{q}_{a}{\alpha }_{v}+{\theta }_{a,v,d}^{U,wane}-({\lambda }_{a}^{v}+\gamma ){U}_{a}^{v,d}+\nu {A}_{a}^{v,d}+{r}_{v}({\theta }_{a,v,d}^{VoS}+{\theta }_{a,v,d}^{VaN})+\frac{1}{D}{U}_{a}^{v,1}\,{\text{, if}}\,d=0\quad \\ {\theta }_{a,v,d}^{U,wane}-({\lambda }_{a}^{v}+\gamma ){U}_{a}^{v,d}+\nu {A}_{a}^{v,d}+{r}_{v}({\theta }_{a,v,d}^{VoS}+{\theta }_{a,v,d}^{VaN}+\mathop{\sum}\limits _{d^{\prime} }{\theta }_{a,v,d{\prime} }^{VoD})-\frac{1}{D}{U}_{a}^{v,d}+\rho \mathop{\sum}\limits _{d^{\prime} }{T}_{a}^{v,d{\prime} }\,{\text{, if}}\,d=1\quad \end{array}\right.$$
(1)
$$\frac{d{I}_{a}^{v,d}}{dt}=\left\{\begin{array}{l}{\theta }_{a,v,d}^{I,wane}+{\lambda }_{a}^{v}{U}_{a}^{v,d}-(\sigma +\gamma ){I}_{a}^{v,d}+\frac{1}{D}{I}_{a}^{v,1}\,{\text{, if}}\,d=0\quad \\ {\theta }_{a,v,d}^{I,wane}+{\lambda }_{a}^{v}{U}_{a}^{v,d}-(\sigma +\gamma ){I}_{a}^{v,d}-\frac{1}{D}{I}_{a}^{v,d}\,{\text{, if}}\,d=1\quad \end{array}\right.$$
(2)
$$\frac{d{A}_{a}^{v,d}}{dt}=\left\{\begin{array}{l}{\theta }_{a,v,d}^{A,wane}+(1-\psi )\sigma {I}_{a}^{v,d}-(\nu +{\eta }_{a}+\gamma ){A}_{a}^{v,d}+\frac{1}{D}{A}_{a}^{v,1}\,{\text{, if}}\,d=0\quad \\ {\theta }_{a,v,d}^{A,wane}+(1-\psi )\sigma {I}_{a}^{v,d}-(\nu +{\eta }_{a}+\gamma ){A}_{a}^{v,d}-\frac{1}{D}{A}_{a}^{v,d}\,{\text{, if}}\,d=1\quad \end{array}\right.$$
(3)
$$\frac{d{S}_{a}^{v,d}}{dt}=\left\{\begin{array}{l}{\theta }_{a,v,d}^{S,wane}+\psi \sigma {I}_{a}^{v,d}-(\mu +\gamma ){S}_{a}^{v,d}+\frac{1}{D}{S}_{a}^{v,1}\,{\text{, if}}\,d=0\quad \\ {\theta }_{a,v,d}^{S,wane}+\psi \sigma {I}_{a}^{v,d}-(\mu +\gamma ){S}_{a}^{v,d}-\frac{1}{D}{S}_{a}^{v,d}\,{\text{, if}}\,d=1\quad \end{array}\right.$$
(4)
$$\frac{d{T}_{a}^{v,d}}{dt}=\left\{\begin{array}{l}{\theta }_{a,v,d}^{T,wane}+{\eta }_{a}{A}_{a}^{v,d}+\mu {S}_{a}^{v,d}-(\rho +\gamma ){T}_{a}^{v,d}+\frac{1}{D}{T}_{a}^{v,1}\,{\text{, if}}\,d=0\quad \\ {\theta }_{a,v,d}^{T,wane}+{\eta }_{a}{A}_{a}^{v,d}+\mu {S}_{a}^{v,d}-(\rho +\gamma ){T}_{a}^{v,d}-\frac{1}{D}{T}_{a}^{v,d}\,{\text{, if}}\,d=1\quad \end{array}\right.$$
(5)
$${\alpha }_{v}=\left\{\begin{array}{ll}(1-h)\alpha \quad &\,{\text{if}}\,v=X\\ 0\quad &\,{\text{if}}\,v\in \{P,{V}_{i},{V}_{e},W,{R}_{i},{R}_{e}\}\\ h\alpha \quad &\,{\text{if}}\,v=H\end{array}\right.$$
(6)
$${\beta }_{v}=\left\{\begin{array}{ll}\beta \quad &\,{\text{if}}\,v\in \{X,W,H\}\\ (1-{e}_{p})\beta \quad &\,{\text{if}}\,v=P\\ (1-{e}_{v})\beta \quad &\,{\text{if}}\,v\in \{{V}_{i},{V}_{e},{R}_{i},{R}_{e}\}\end{array}\right.$$
(7)
$${\theta }_{a,v,d}^{Z,wane}=\left\{\begin{array}{ll}\frac{1}{{D}_{P}}{Z}_{a}^{P,d}\quad &\,{\text{if}}\,v=X\\ -\frac{1}{{D}_{P}}{Z}_{a}^{P,d}\quad &\,{\text{if}}\,v=P\\ -\frac{1}{{E}_{V}}{Z}_{a}^{{V}_{i},d}\quad &\,{\text{if}}\,v={V}_{i}\\ \frac{1}{{E}_{V}}{Z}_{a}^{{V}_{i},d}-\frac{1}{{D}_{V}-{E}_{V}}{Z}_{a}^{{V}_{e},d}\quad &\,{\text{if}}\,v={V}_{e}\\ \frac{1}{{D}_{V}-{E}_{V}}{Z}_{a}^{{V}_{e},d}+\frac{1}{{D}_{R}-{E}_{R}}{Z}_{a}^{{R}_{e},d}\quad &\,{\text{if}}\,v=W\\ -\frac{1}{{E}_{R}}{Z}_{a}^{{R}_{i},d}\quad &\,{\text{if}}\,v={R}_{i}\\ \frac{1}{{E}_{R}}{Z}_{a}^{{R}_{i},d}-\frac{1}{{D}_{R}-{E}_{R}}{Z}_{a}^{{R}_{e},d}\quad &\,{\text{if}}\,v={R}_{e}\\ 0\quad &\,{\text{if}}\,v=H\end{array}\right.$$
(8)
$${\theta }_{a,v,d}^{VoS}={\eta }_{a}{\delta }_{a,d}^{VoS}\times \left\{\begin{array}{ll}-{r}_{1}{U}_{a}^{X,d}\quad &\,{\text{if}}\,v=X\\ {r}_{1}(1-{r}_{2}){U}_{a}^{X,d}-{r}_{1}{r}_{2}{U}_{a}^{P,d}\quad &\,{\text{if}}\,v=P\\ {r}_{1}{r}_{2}{U}_{a}^{X,d}+{r}_{1}{r}_{2}{U}_{a}^{P,d}\quad &\,{\text{if}}\,v={V}_{i}\\ -{r}_{b}{U}_{a}^{v,d}\quad &\,{\text{if}}\,v\in \{{V}_{e},W,{R}_{e}\}\\ {r}_{b}{U}_{a}^{W,d}+{r}_{b}{U}_{a}^{{V}_{e},d}+{r}_{b}{U}_{a}^{{R}_{e},d}\quad &\,{\text{if}}\,v={R}_{i}\\ 0\quad &\,{\text{if}}\,v=H\end{array}\right.$$
(9)
$${\theta }_{a,v,d}^{VaN}={\delta }^{VaN}\times \left\{\begin{array}{ll}-{r}_{1}{\phi }_{a}^{X,d}\quad &\,{\text{if}}\,v=X\\ {r}_{1}(1-{r}_{2}){\phi }_{a}^{X,d}-{r}_{1}{r}_{2}{\phi }_{a}^{P,d}\quad &\,{\text{if}}\,v=P\\ {r}_{1}{r}_{2}{\phi }_{a}^{X,d}+{r}_{1}{r}_{2}{\phi }_{a}^{P,d}\quad &\,{\text{if}}\,v={V}_{i}\\ -{r}_{b}{\phi }_{a}^{v,d}\quad &\,{\text{if}}\,v\in \{{V}_{e},W,{R}_{e}\}\\ {r}_{b}{\phi }_{a}^{W,d}+{r}_{b}{\phi }_{a}^{{V}_{e},d}+{r}_{b}{\phi }_{a}^{{R}_{e},d}\quad &\,{\text{if}}\,v={R}_{i}\\ 0\quad &\,{\text{if}}\,v=H\end{array}\right.$$
(10)
$${\theta }_{a,v,d}^{VoD}=\rho {\delta }^{VoD}\times \left\{\begin{array}{ll}-{r}_{1}{T}_{a}^{X,d}\quad &\,{\text{if}}\,v=X\\ {r}_{1}(1-{r}_{2}){T}_{a}^{X,d}-{r}_{1}{r}_{2}{T}_{a}^{P,d}\quad &\,{\text{if}}\,v=P\\ {r}_{1}{r}_{2}{T}_{a}^{X,d}+{r}_{1}{r}_{2}{T}_{a}^{P,d}\quad &\,{\text{if}}\,v={V}_{i}\\ -{r}_{b}{T}_{a}^{v,d}\quad &\,{\text{if}}\,v\in \{{V}_{e},W,{R}_{e}\}\\ {r}_{b}{T}_{a}^{W,d}+{r}_{b}{T}_{a}^{{V}_{e},d}+{r}_{b}{T}_{a}^{{R}_{e},d}\quad &\,{\text{if}}\,v={R}_{i}\\ 0\quad &\,{\text{if}}\,v=H\end{array}\right.$$
(11)
Letting \({C}_{a}^{v,d}={I}_{a}^{v,d}+{A}_{a}^{v,d}+{S}_{a}^{v,d}\), \({C}_{a}={\sum }_{v}{\sum }_{d}{C}_{a}^{v,d}\), \({N}_{a}^{v,d}={U}_{a}^{v,d}+{I}_{a}^{v,d}+{A}_{a}^{v,d}+{S}_{a}^{v,d}+{T}_{a}^{v,d}\), \({N}_{a}={\sum }_{v}{\sum }_{d}{N}_{a}^{v,d}\), with ca denoting the partner change rate and with ϵ denoting the level of assortative mixing between activity groups, the force of infection to an individual in activity group a and vaccination stratum v is given by:
$$\begin{array}{lll}{\lambda }_{a}^{v}={c}_{a}{\beta }_{v}\left(\epsilon \frac{{C}_{a}}{{N}_{a}}+(1-\epsilon )\frac{\sum _{a{\prime} }{c}_{a{\prime} }{C}_{a{\prime} }}{\sum _{a{\prime} }{c}_{a{\prime} }{N}_{a{\prime} }}\right)\end{array}$$
(12)
Partner notification
Implementing a strategy which offers vaccination to uninfected notified partners within the model requires quantifying the rate that uninfected individuals in activity group a, vaccination stratum v, and diagnosis history stratum d attend SHCs through partner notification. This is derived as follows:
-
1.
The overall rate of diagnosis from activity group b is given by \(\sum _{v}\sum _{d}\mu {S}_{b}^{v,d}+{\eta }_{b}{A}_{b}^{v,d}\)
-
2.
For each individual who is diagnosed, κ individuals are notified.
-
3.
Of the notified individuals, a proportion p are infected.
-
4.
Of uninfected notified partners, ωU(b, a) individuals are from activity group a, where ωU is given by:
$${\omega }_{U}=\left[\begin{array}{cc}\epsilon +(1-\epsilon )\frac{{c}_{L}{U}_{L}}{\sum _{a{\prime} }{c}_{a{\prime} }{U}_{a{\prime} }}&(1-\epsilon )\frac{{c}_{H}{U}_{H}}{\sum _{a{\prime} }{c}_{a{\prime} }{U}_{a{\prime} }}\\ (1-\epsilon )\frac{{c}_{L}{U}_{L}}{\sum _{a{\prime} }{c}_{a{\prime} }{U}_{a{\prime} }}&\epsilon +(1-\epsilon )\frac{{c}_{H}{U}_{H}}{\sum _{a{\prime} }{c}_{a{\prime} }{U}_{a{\prime} }}\\ \end{array}\right]$$
(13)
-
5.
The proportion of uninfected notified partners in activity group a who are from vaccination stratum v and with diagnosis history d is given by \(\frac{{U}_{a}^{v,d}}{{\sum }_{v{\prime} \in {\mathcal{V}}}{\sum }_{d{\prime} \in {\mathcal{D}}}{U}_{a}^{v{\prime} ,d{\prime} }}\).
-
6.
So, the overall rate of uninfected individuals in subgroup (a, v, d) being notified from a contact in activity group b is given by
$$\kappa \times (1-p)\times \frac{{U}_{a}^{v,d}}{\sum _{v{\prime} }\sum _{d{\prime} }{U}_{a}^{v{\prime} ,d{\prime} }}\times {\omega }_{U}(b,a)\times \sum _{v{\prime\prime} }\sum _{d{\prime\prime} }\left({{\mu {S}_{b}^{v}}^{\prime\prime} ,d^{\prime\prime} }+{\eta }_{b}{A}_{b}^{v{\prime\prime} ,d{\prime\prime} }\right)$$
(14)
-
7.
So, the overall rate of uninfected individuals in subgroup (a, v, d) being notified by someone in any group is given by:
$${\phi }_{a}^{v,d}=\kappa \times (1-p)\times \frac{{U}_{a}^{v,d}}{\sum _{v{\prime} }\sum _{d{\prime} }{U}_{a}^{v{\prime} ,d{\prime} }}\times \sum _{a{\prime\prime} }\sum _{v{\prime\prime} }\sum _{d{\prime\prime} }{\omega }_{U}(a{\prime\prime} ,a)\left(\mu {S}_{a{\prime\prime} }^{v{\prime\prime} ,d{\prime\prime} }+{\eta }_{a{\prime\prime} }{A}_{a{\prime\prime} }^{v{\prime\prime} ,d{\prime\prime} }\right)$$
(15)
This approach allows us to capture partner notification within a compartmental model provided that estimates of κ and p can be obtained.
Note that infected individuals who attend SHCs through partner notification are already implicitly captured within our framework, as vaccination is offered on diagnosis, and the model is calibrated to diagnosis numbers in England (which includes diagnoses through partner notification).
Calibration and accounting for uncertainty
Epidemiological parameters are estimated via Bayesian calibration to data from the UK’s GUMCAD and GRASP surveillance systems7,40. We account for uncertainty in estimated parameters by using 1000 samples from the joint posterior distribution. Full details of the model calibration, including plots of model trajectories in comparison with the observations used for calibration, are provided in previous work10.
We introduce two parameters to incorporate partner notification: κ, the number of individuals attending SHC via partner notification per index case (0.39 [95% CrI:0.31–0.47])7, and p, the prevalence of gonorrhea infection among notified partners (38% [95% CrI:31–45%])41.
We account for uncertainty in the interpretation of vaccine uptake data42 by considering two vaccine-sentiment scenarios. In our main analysis, all MSM have the same probability of accepting vaccination, as in most modeling analyses. Vaccine uptake (the probability of acceptance when offered) is set at the level observed for MSM HPV vaccination in England42, i.e., 40.8% (95% CrI:40.6–41.0). In supplementary analysis, we consider an alternative scenario in which a proportion of individuals are vaccine-willing and the remainder are unwilling and never accept vaccination12. In this scenario, the proportion of individuals who are willing to be vaccinated is set to equal the probability of acceptance in our main scenario, to obtain equal initial rates of vaccination in the two scenarios12. In both scenarios, 61.7% (95% CrI:61.2–62.1) of those who receive a first dose return to obtain a second dose42. A model schematic of vaccination-status strata, including a stratum for those unwilling to be vaccinated, is shown in Supplementary Fig. 9.
Vaccine-targeting strategies
We compare five vaccine-targeting strategies that could be implemented at SHCs, where the vast majority of gonorrhea diagnoses in England are managed43: VoD, VaH, VaN, VaH, and VaN combined (VaH + VaN), and VaR (eligibility detailed in Table 1). We start the model at equilibrium before each vaccination scenario is implemented.
Primary vaccination requires two doses, but some vaccinees only receive one dose due to not returning for the second primary dose. Those who receive only one dose get two-thirds the protection of two doses20. Revaccination is offered after a period of time since the last vaccination according to the eligibility criteria of the implemented strategy: those who received a two-dose primary vaccination are offered a single booster dose, while those who received only one primary dose are offered another course of primary vaccination.
We evaluate the impact of each strategy by calculating the yearly number of gonorrhea diagnoses averted, compared to a baseline scenario without vaccination. We evaluate the efficiency of each strategy by calculating the cumulative number of gonorrhea diagnoses averted per dose administered.
In scenario analysis, as before10 we vary the level of protection from two primary doses (20 and 40%) and the duration of protection (1.5 years after primary then 3 years after booster vaccination, or 4 years after primary and booster vaccination, or 7.5 years after primary and booster vaccination).
Health-economic analysis
Health-economic analysis takes the perspective of the UK National Health Service (NHS), considering (i) costs of gonorrhea testing, treatment, and vaccination, and (ii) loss of health due to symptoms, measured in quality-adjusted life years (QALYs). Costs and QALYs are calculated as in Whittles et al. 202210; a brief description of the approach is included below, and health-economic parameters are detailed in Table 5. For the convenience of the reader, equations for model outputs required to evaluate impact, efficiency, and cost-effectiveness are detailed in Table 6.
Costs (£2021-22 GBP) and QALYs are discounted at 3.5% per annum44. In our main analysis, we consider a vaccination cost of £18/dose administered, which includes a £10 administration fee10,45. We also consider a vaccination cost of £85/dose administered—the vaccine’s list price46 plus £10 administration fee.
We calculate the net monetary benefit (NMB) of each vaccination strategy, i.e., the monetary value of averting QALY losses minus the net costs of vaccination, compared with no vaccination. QALY loss is the product of the Quality-of-Life disutility of symptomatic infection47 and the average duration of symptoms. For vaccination to be considered cost-effective in the UK, there needs to be at least 50% probability that NMB is positive with a QALY valued at £20,000 and at least 90% probability with a QALY valued at £30,00044. Therefore, we calculate NMB with a QALY valued at £20,000 (NMB £20k) and at £30,000 (NMB £30k).
We calculate both impact and cost-effectiveness over a 10-year time horizon; in supplementary analysis, we also consider a 20-year time horizon.
