Assessment of the Potential of Vaccination to Combat Antibiotic Resistance in Gonorrhea: A Modeling Analysis to Determine Preferred Product Characteristics

Abstract Background Gonorrhea incidence is increasing rapidly in many countries, while antibiotic resistance is making treatment more difficult. Combined with evidence that two meningococcal vaccines are likely partially protective against gonorrhea, this has renewed interest in a gonococcal vaccine, and several candidates are in development. Key questions are how protective and long-lasting a vaccine needs to be, and how to target it. We assessed vaccination’s potential impact and the feasibility of achieving the World Health Organization’s (WHO) target of reducing gonorrhea incidence by 90% during 2018–2030, by comparing realistic vaccination strategies under a range of scenarios of vaccine efficacy and duration of protection, and emergence of extensively-resistant gonorrhea. Methods We developed a stochastic transmission-dynamic model, incorporating asymptomatic and symptomatic infection and heterogeneous sexual behavior in men who have sex with men (MSM). We used data from England, which has a comprehensive, consistent nationwide surveillance system. Using particle Markov chain Monte Carlo methods, we fitted to gonorrhea incidence in 2008–2017, then used Bayesian forecasting to examine an extensive range of scenarios. Results Even in the worst-case scenario of untreatable infection emerging, the WHO target is achievable if all MSM attending sexual health clinics receive a vaccine offering ≥ 52% protection for ≥ 6 years. A vaccine conferring 31% protection (as estimated for MeNZB) for 2–4 years could reduce incidence in 2030 by 45% in the worst-case scenario, and by 75% if > 70% of resistant gonorrhea remains treatable. Conclusions Even a partially-protective vaccine, delivered through a realistic targeting strategy, could substantially reduce gonorrhea incidence, despite antibiotic resistance.

We adapted a previously published stochastic compartmental model of gonorrhea [3] to allow for sexual heterogeneity by separating the population into two groups based on sexual activity levels g ∈ {H, L} and adding vaccine protection . The full model is illustrated in Figure S1 with notations summarised in Supplementary Tables S1 and S2.
The initial population size was N = 600,000, with rates of population entry and exit calibrated based on data from the UK Office for National Statistics (Supplementary Figure S2) [3,4]. The number of MSM in each activity group is denoted N g (t) (N L (t) + N H (t) = N (t)). Individuals remain in the sexual activity group to which they are initially assigned throughout the model simulations. The proportion of the population in each sexual activity group and average number of partners of those in each activity group were calibrated according to the Natsal-3 data. The 85% of MSM having up to five partners a year were assigned to the group L, with the remaining 15% who had more than five partners per year assigned to group H.
The rate of infection for uninfected individuals depends on the probability of infection per partnership (β), the average partner change rate of the activity group (p g ) and the sexual mixing pattern ( ), i.e the number of partnerships formed within and between activity groups. The parameter represents the sexual mixing coefficient and ranges from 0 to 1, representing a range of scenarios from proportionate mixing, where partners are chosen at random regardless of their sexual activity group, to assortative mixing, where partners are preferentially chosen from an individual's own sexual activity group [5].  Figure 1 in the main paper, but with arrows labelled with algebraic expressions, as defined in Table S1.   The rate of transmission from an individual in group h to an uninfected individual in group g is denoted θ gh and defined as follows: Individuals in group g are initially uninfected (U g ). They become infected with strain s ∈ {ABS, ABR}, denoting antibiotic susceptible (ABS) and ABR strains respectively. The model assumes that strains do not vary in transmissibility, propensity to cause symptoms or rate of natural clearance, and that the rate of infection from a contagious individual in group h to an uninfected individual in group g is θ gh (t). Infected individuals initially pass through an incubation period (I s g ), which they leave at rate σ. A proportion ψ of those infected then go on to develop symptoms (S s g ), whereas the rest remain asymptomatic (A s g ). Symptomatic individuals (S s g ) seek care at rate µ. Individuals without symptoms attend sexual health screening at rate η. Natural clearance of asymptomatic infection happens at rate ν. Treated individuals recover from the infection and become uninfected again at rate ρ, with the exception of a proportion φ of the individuals infected with an ABR strain (T ABR g ) for whom treatment fails, leading to persistent infection for which the same dynamics are assumed as for asymptomatic cases. Gonococcal infection can occur in the rectum, pharynx and/or urethra, resulting in different rates of onward transmission and probabilities of developing symptoms [6]. Surveillance data are not stratified by infection site (rectum, pharynx, urethra), hence the parameters are averages across sites. The model assumes that no natural immunity is conferred by infection [7,8]. The ABR strain is imported into the highly sexually active group at rate ω/N H (t) per person per year from time t ABR =2020 onwards.

Interventions
We considered three potential vaccination strategies in the model: • Vaccination immediately before entry into the population -where a proportion p vax E > 0 of adolescents are vaccinated before they become sexually active; • Vaccination on diagnosis -where the vaccine is given to a proportion p vax D > 0 of MSM diagnosed with gonorrhea; • Vaccination on attendance -where the vaccine is given to a proportion p vax S = p vax D > 0 of MSM that present to sexual health services, either through screening or due to symptoms.
We varied vaccine efficacy in the analysis, with factor γ denoting the relative reduction in susceptibility to gonorrhea conferred by the vaccine, which is known as degree-type protection [9]. As such, vaccinated MSM (Û g ) can still become infected, albeit less frequently, on sexual contact with contagious individuals (C). These vaccinated-yet-infected individuals progress through the stages of infection (Î g ,Ŝ g ,Â g ) similarly to those who are unvaccinated. The contagious population for each strain s is denoted C s g = I s g + S s g + A s g +Î s g +Ŝ s g +Â s g . Under the vaccination-on-diagnosis and vaccination-on-attendance strategies, treatment fails for a proportion φ of those infected with the ABR strain, who remain contagious, despite receiving the vaccine. Vaccine protection wanes at rate ζ per year, after which time individuals become fully susceptible to infection once more.

Model equations and stochastic simulations
In order to illustrate the model dynamics, we describe a deterministic version of the model by the differential equations S5 -S20, which are easier to read than the stochastic specification of the model, which follows.
dT s g (t) dt = µS s g (t) + ηA s g (t) − ρ + χ(t) T s g (t) + ζT s g (t) (S12) − γ θLLCL(t) + θLH CH (t) + ζ + χ(t) Û L(t) We used a stochastic version of the model described in Equations S5 -S20, time-discretised using the Gillespie tau-leap method with steps of 1 day. The process was initialised on 31 December 2007, with A ABS g (0) asymptomatic infections in each sexual activity group g ∈ {L, H} and the remainder of the population uninfected. Simulation proceeded through iterations of the steps below for each day. Since the model is stochastic, the outbreak can become extinct even when its basic reproduction number is >1. Transition variables b g 1 , b g 2 , d g 1 , . . . , d g 32 are drawn from the following distributions: θHLC ABR  For ease of notation, we have expressed the equations above using the rates we have previously defined, rather than their respective probabilities. Therefore, for example, d ∼ Multinom(S, r) should be understood as d ∼ Multinom(S, 1 − e −r ).
The compartments of the model were updated as follows:

Bayesian Inference
We considered the recorded cases (Y D (t), t = 2008, ..., 2017) as the observed realisations of an underlying unobserved Markov process: the total number of diagnosed gonococcal infections (Z D (t)) (Supplementary Figure S3). We allowed for a 10% margin of under-reporting, such that 90% of all gonorrhea diagnoses are recorded in the surveillance data -consistent with evidence of the proportion of gonorrhea diagnoses made in a General Practice setting and therefore not reported [2].
Additionally, to ensure that the rate at which MSM attend sexual health screening reflects the most recent data, we considered the number of sexual health clinic attendances in 2017 (Y A (t), t = 2017) as the observed realisation of an underlying unobserved variable -the total number of MSM attending sexual health services (Z A (2017)), again allowing for a 10% margin of under-reporting.
... ... An expression for the likelihood of the observed data given the model is not analytically tractable, so we used a particle filter to compute an unbiased estimate. This estimated likelihood was then incorporated into a Monte Carlo Markov Chain (MCMC) algorithm, which was used to sample from the joint posterior distribution of the parameters given the observed data [10].
As the number of parameters increases, so does the difficulty of exploring the joint posterior space. To be parsimonious we reduced the number of model parameters by making the initial number of symptomatic cases, S L (2007) and S H (2007), to be zero. After a few simulated days these variables quickly reach the stochastic values implied by the model and parameters. Therefore, the first few days of each simulation can be thought of as the equivalent to a Markov Chain convergence or burn-in step.
We implemented the model fitting using a pMCMC algorithm from the R package pomp [11], which we modified to enable parallel computation. The estimation of the likelihood was based on 1,000 particles, which we found to be sufficiently robust to produce consistent estimates of the likelihood. Four pMCMC with dispersed starting points were run for 110,000 iterations, with the first 10% discarded as burn-in. The chains were compared using the R package coda [12]; the multivariate Gelman-Rubin diagnostic was ≤1.1 for all inferred parameters, suggesting that they appear to have converged on the same posterior distribution [13]. For robustness, we combined the posterior samples from all four chains and ensured all parameters had an effective sample size greater than 200.

Prior distributions of parameters
Performing Bayesian inference requires the specification of plausible parameter priors. For parameters whose values are uncertain ( , , β, ψ) we adopted highly uninformative uniform priors (Table S2). The remaining four parameters ν, σ, µ and ρ were assigned informative priors based on literature review, as summarised in Table S2.
The length of time for which asymptomatic gonorrhea infection persists has not been well measured. Duration of carriage likely depends on the site of infection, with estimates varying widely from 12 weeks at the pharynx, up to a year for rectal infection [14,15]. In a study of 18 asymptomatic men infected with (urethral) gonorrhea, no infections were cleared in the 165 days from diagnosis to treatment [16]. Modelling studies commonly use an average duration of untreated infection of 6 months [5,17,18]. This is supported by evidence from genomic data, where bacterial genomes isolated from known contact pairs had a maximum time to most recent common ancestor of 8 months [19]. We therefore assigned the duration of untreated carriage, 1 ν , a prior corresponding to a mean duration of carriage between three months and one year with 95% prior weight.
The UK national guidelines on safer sex advice recommend MSM test for gonorrhea annually, and that highly sexually active MSM get tested every three months [20,21]. However, in a social media survey of 2,668 MSM in Scotland, Wales and Ireland, two-thirds reported STI testing less than once a year [22]. We assigned a uniform distribution to the rate of screening, η with 95% prior weight between 0.4 and 4. This corresponds to a range of scenarios: from testing at quarterly intervals, in adherence with the British Association for Sexual Health and HIV guidelines, through to the situation described by [22]. The lower bound for η was calculated by setting the probability of testing less than once a year equal to 2 3 under an exponential distribution with rate η, so that P[T > 1] = e −η = 2 3 , and solving to find The length of time spent in the incubation, symptomatic, and treatment stages of infection are observed to be short [23][24][25]. We set priors accordingly: the mean duration of incubation, 1 σ , was between three and eight days; the mean duration of symptomatic infection (before seeking care), 1 µ , was between a day and two weeks; and the mean duration of treatment, 1 ρ , was between five and nine days, in accordance with the BASHH recommendation to abstain from sex for seven days after treatment [26].

Estimation of model parameters
The posterior estimate of the probability of transmission per partnership per year β was 59% with a wide range of credible values from 28% to 97%. The parameter was positively correlated with the level of assortativity in sexual mixing , which was estimated to be between 0.06 and 0.75 with 95% credibility. This corresponds to the trade-off required to maintain a force of infection compatible with the observed data, as described by equations S1-4.
The model predicts that 25% (95% CrI: 13%, 35%) of infections become symptomatic, and that this proportion (ψ) is strongly negatively correlated with the rate of asymptomatic screening (η). This corresponds to the tradeoff between duration of carriage, 1 ν+η , and the proportion of infections entering the carriage state (1 − ψ) The parameters corresponding to the time spent in the incubation, symptomatic, and treatment stages of infection (σ, µ, ρ) had posterior distributions that closely matched their prior distributions, suggesting that the priors were appropriate and that there is little additional information on these parameters to be gleaned from this data set. These parameters did not show a strong correlation with any others; as expected given the short duration of the incubation, symptomatic and treatment stages variation in their parameters makes little difference to the dynamics of infection.
The posterior distribution of ν has a slightly lower mean than the prior, implying a longer mean duration of carriage: 174 (95% CrI: 95, 443) days compared to 157 (95% CrI: 87, 364). The prior and posterior credible intervals intersect to a large extent so there is not significant evidence of a departure from the prior based on the data.     Figure S4. Posterior distributions of parameters. The diagonal plots show histograms of the posterior distributions for all sampled parameters. The blue lines show prior distributions, and the red lines indicate posterior mean and 95% credible intervals. The plots below the diagonal show scatter plots based on 1,000 samples from the posterior, illustrating the relationships between pairs of estimated parameters. An orange background indicates a correlation higher than 0.8, a yellow background indicates a correlation between 0.5 and 0.8, a green background indicates a correlation between 0.3 and 0.5, and a white background indicates a correlation less than 0.3. The plots above the diagonal show the corresponding correlation coefficients with the 95% credible intervals in parentheses.

Overview of simulation analysis
To examine the effect of vaccine efficacy and duration of protection, and the effect of different deployment strategies, we considered 1,000 vaccine profiles (i.e. combinations of level and duration of protection), selected using Latin Hypercube Sampling, which ranged in efficacy from 1-100% and offered protection lasting between 1 and 20 years. For each vaccine profile we performed 100 Monte Carlo simulations of the epidemic between 2008 and 2030. we repeated the analysis for each of the ABR strains of varying controllability, under each of the three vaccination strategies.
We assessed the potential impact of a MeNZB-like vaccine by considering 100 vaccine profiles conferring protection from infection ranging from 21%-39% with an assumed duration of two to four years. We repeated the analysis for varying levels of vaccine coverage to assess how the level of uptake achieved impacts the success of each strategy.