Mathematical Modelling of Gonorrhoea Spread in Northern Ireland between 2012 and 2022

Abstract

The number of confirmed positive tests of various sexually transmitted infections has grown recently in the United Kingdom. The objective of this study is to propose a deterministic compartmental model to investigate gonorrhoea spread in Northern Ireland between 2012 and 2022. The differential equation based model includes both symptomatic and asymptomatic spread, spontaneous recovery and treatment compartments. After fitting our model to the monthly number of new positive tests, we found that the basic reproduction number is approximately $1.0030$. In addition, we derive the endemic equilibrium of the model, which exists if and only if ${\mathcal{R}}_0>1$. The sensitivity analyses of the basic reproduction number and the endemic values of the compartments of treated individuals indicate that infection spreading time can have a significant impact on gonorrhoea spread.

1. Introduction

1.1. Gonorrhoea

The number of confirmed gonorrhoea positive tests is growing in the United Kingdom (UK). Most recently, it was reported that gonorrhoea and syphilis sex infections reached record levels in England [1,2]. Furthermore, the number of people diagnosed with gonorrhoea in Scotland has doubled over the last 5 years [3], or in Northern Ireland it was the highest on record already in 2019 [4]. However, gonorrhoea is not a public health problem specific only to the UK. Among adults globally, The World Health Organization (WHO) estimates that there were 82.4 million new infections in 2020. Using this estimate as a baseline, it aims at a 90% reduction in the incidence of Neisseria gonorrhoeae infections by 2030—as a part of its global health sector strategies on HIV, viral hepatitis and STIs 2022–2030 [5].

Although it can spread perinatally from mother to baby during childbirth, gonorrhoea is a predominantly sexually transmitted infection (STI) caused by a species of bacteria Neisseria gonorrhoeae (NG) which affects only humans. In a sexually active population, NG can be transferred from an infectious individual to a susceptible one via sexual contact. The potential sites of infection, depending on sexual practices, include the genitals, mouth, or rectum. After exposure to NG, the incubation period (time from infection to symptoms) of gonorrhoea is usually two weeks [3], which, similarly to other STIs, is shorter than the infectious period [6] (Ch10.3, p. 329). A common symptom of a urogenital infection is pain or a burning sensation when passing urine. Other symptoms can include unusual vaginal or penile discharge and inflammation of the foreskin. Less frequent symptoms in women are pain or tenderness in the lower abdominal area and bleeding between periods, heavier periods and bleeding after sex, and in men, there can be pain or tenderness in the testicles. However, symptoms sometimes do not appear until many months later, if at all. Infection is possible during the incubation period or even without the presence of any of the symptoms. About 1 in 10 infected men and 5 in 10 infected women will not experience any obvious symptoms [7,8,9,10], which also means that the condition can go untreated for some time. The most common asymptomatic cases of gonorrhoea are rectal and pharyngeal infections and mostly identified in men who have sex with men (MSM) [11].  Spontaneous clearance of genital and extragenital NG is also possible, and a recent, large multicentre UK cohort based report indicates that clearance of infection occurred in $20.5\%$ of 405 confirmed cases over a median of 10 days [12]. However, the estimated probability of spontaneous clearance depends on the site of infection and varies between 0.05 and 0.57 [13,14,15,16]. The duration of infectivity is complex and dependent on available health services [17]. Although relatively rare, exposure to NG can lead to disseminated gonococcal infection with symptoms including bacteraemia, fever, dermatitis, tenosynovitis, septic arthritis, endocarditis and meningitis [18,19].

Gonorrhoea testing can be performed by a culture using a sample collected from the site of potential infection; this also can identify any antibiotic resistance [20]. However, since it has higher sensitivity than cultures, Nucleic Acid Amplification Testing is used much more often nowadays. A urine based test is available for men. Patients with gonorrhoea are advised to avoid having sex until they, and their partner, have been treated and tested for being cured to prevent re-infection or treatment failure.

Exposure to NG does not offer protection against reinfection, that is, sexual contact with a person infected with NG can lead to gonorrhoea in previously infected or treated individuals [21]. Although the probability of spontaneous clearance may increase over time, since the urogenital infection may move up to the upper genital tract, an untreated NG infection may lead to many severe complications, such as endometritis, penile oedema, and epididymitis, resulting in infertility. In addition, many untreated bacterial infections spreading from the vagina or the cervix to the reproductive organs can result in pelvic inflammatory disease (PID) or involuntary loss of life through ectopic pregnancy [20,22] (p. 588, [22]). Although the exact numbers are not known, in the UK, it is estimated that around 220,000 women are affected each year by PID, mainly from the age group between fifteen and twenty four [23]. PID is a serious common condition, meaning that in sexually active women, a low threshold for diagnosing should be applied and antimicrobial treatment should not be withheld while waiting for STI testing results [24]. Besides pain and discomfort, based on the chlamydia and gonorrhoea focused study in [25], a delay of 3 or more days while seeking treatment after the onset of lower abdominal pain can result in a 3-fold increased risk of infertility or ectopic pregnancy. Furthermore, NG, together with Chlamydia trachomatis, is classified as the second most likely microbial cause of subclinical PID [26].

Antimicrobial resistance in NG is a significant public health concern [27]. Gonorrhoea has consistently developed resistance to multiple different antimicrobials. During the 1940s, sulfonamides were used to treat gonorrhoea, but already by 1944, many gonococcal strains showed resistance. The mass production of penicillin after 1944 provided an effective treatment although requiring gradually increased doses over the years because of the novel NG strains, eventually leading to treatment failure. NG has the ability to transfer partial or whole genes during its entire life cycle and can also change its genome through all types of mutations resulting in decreasing efficacy of treatments with quinolones, macrolides and cephalosporins, even with a high rate of initial success. Most recently, the only available first line monotherapy was ceftriaxone, but because of the decreasing susceptibility of NG, dual-antimicrobial therapy was introduced in the US, the UK and Europe by adding low-dosage azithromycin [22,28]. The failure of the dual therapy was reported first in the UK in 2016 [29]. Based on surveillance data spanning over 15 years, the use of antibiotics, and epidemiological trends of antibiotic resistance in gonorrhoea in the UK are reviewed in [30] where possible actions to enhance surveillance of gonorrhoea antibiotic resistance are also discussed. These trends are regularly reported and reviewed in the UK. Similar trends and challenges in the US are reported in [31,32]. The failure of dual therapy with high-dosage azithromycin was reported in England [33] and in Sweden [34]. As a results of these trends, monotherapy with ceftriaxone is used currently in the US and in the UK, [35,36]. Because of the diminishing options for gonorrhoea treatment [37,38,39,40], alternative treatment methods are continually investigated [41]. The WHO is implementing a Global Action Plan to Control the Spread and Impact of Antimicrobial Resistance in NG with actions including effective prevention and control of gonococcal infections. This will use prevention messages and interventions and appropriate treatment regimens, strengthened surveillance systems for antimicrobial resistance and research into alternative treatments for gonococcal infections [42].

1.2. Mathematical Models of Spread of Gonorrhoea and Sexually Transmitted Infections

Various modelling approaches, such as Markov Models, Agent-Based Models, Network Models, System Dynamics Models, Hybrid Models, Microsimulation Models, and Partial Differential Equations, have been used in modelling different aspects of gonorrhoea spread. One of the earliest compartmental models of gonorrhoea in [43] uses a multigroup susceptible–infectious partition of a nonhomogeneous population where the subgroups are homogeneous regarding contact and recovery rates, and it is shown that the disease is endemic. In [44], based on modelling, it is argued that the existing epidemiological data and behavioural data best fit their oropharyngeal model when studying spread of gonorrhoea within men who have sex with men (MSM). Using stochastic simulations and time-dependent incidence matrices, the effects of sexual activity, mixing patterns and network structure were investigated in [17]. Distinguishing aspects of sexually transmitted infection spread have also been modelled. For instance, in [45], by considering explicitly a paring rate and a separation rate (that is the rates at which a sexual partnership forms and dissolves into two singles, respectively), it is shown that endemic equilibria can only exist if the separation rate is sufficiently large in order to ensure the necessary number of sexual partners. The effects of groups of individuals who are sexually very active and are efficient transmitters is studied in [46]. Furthermore, Ref. [47] uses an individual-based Susceptible–Infected–Susceptible–Protected model to determine robust distributions for the rate of new partnerships that involve condomless sex and can therefore facilitate the spread of sexually transmitted infections. In [48], it is shown that a two-sex Susceptible–Infected–Susceptible model cannot support the coexistence of two competing strains in realistic circumstances. The model in [49] is used to capture the dynamics of cefixime resistance observed in England between 2008 and 2015. In addition, it derives estimates of the basic reproduction number of gonorrhoea under various modelling considerations. More recently, Ref. [50] investigates the impact of possible vaccination against gonorrhoea.

1.3. Northern Ireland

Northern Ireland (NI) is on the isle of Ireland, and it shares an open border with the Republic of Ireland, Figure 1. The population in NI grew from 1,685,300 in 2011 to 1,903,200 in 2021 [51]. Although NI is a part of the UK, it has had a devolved government, the Northern Ireland Executive, since 1998. In particular, the Department of Health is responsible for the areas of health and social care, public health and public safety. It also funds the healthcare system, the Health and Social Care (HSC), which is considered as a part of the overall National Health Service (NHS) in the UK. The Public Health Agency (PHA) is the executive agency responsible for the provisions of health and social wellbeing improvement, health protection, public health support to commissioning and policy development, HSC research and development.

The paper is organised as follows. In Section 2, we describe the data set, and we introduce our model, the basic reproduction number and the methods used in the paper. In Section 3, we present our findings illustrated by simulations. In Section 4, we discuss further some of the findings presented in Section 3 together with possible directions of future work. In Section 5, we briefly summarise our findings.

2. Materials and Methods

2.1. Data

Although the population increased between January 2012 and September 2022 in NI, we assume that the population was closed, that is, in our work, we kept it constant. Furthermore, because of the nature of the disease, we consider the subpopulation aged 15 and older as of 2020. Therefore, using the data from [52], we set $N=\mathrm{1,522,883}$. Furthermore, the birth rate was obtained as $\frac{25269}{1522883\times 365}$ where 25,269 is the number of newborns in 2012 [53] the largest in a year over the studied period. (Notice that it means that we assume that every newborn survives up to the age of 15).

We used the monthly number of positive tests reported in Genitourinary Medicine clinics in Northern Ireland between January 2012 and September 2022. The data set was obtained under the Freedom of Information Act 2000, and it is available from the PHA, see [54].

2.2. The Model

Although in modelling sexually transmitted infections, the number and pattern of sexual contacts play an important role, in our proposed model to incorporate various aspects of gonorrhoea spread in Northern Ireland, we assume random mixing to keep the model analysis relatively simple. Nevertheless, our proposed route of transmission can serve as a core to develop more detailed models by incorporating the above mentioned approaches to gain better insights into the spread of sexually transmitted infections. Furthermore, to the best of our knowledge, this is the first study of its kind.

Based on the transmission and the features of Gonorrhoea described in Section 1, we derived the transmission diagram shown in Figure 2. The interpretation of the non-negative parameters is given

Table 1. Parameters. The interpretations and the ranges of the parameters in Figure 2.

Parameter	Meaning	Value	Range
$\beta$	Transmission rate (days −1)	$\frac{2}{10}$	0–1
$\alpha$	Infectiousness reduction in asymptomatic	$\frac{7}{10}$	0–1
${\sigma }_A^{-1}$	Average time from infectiousness to positive test of asymptomatic (days) [10]	40	2–150
${\sigma }_I^{-1}$	Average time from infectiousness to positive test of symptomatic (days)	14	2–14
${\gamma }^{-1}$	Average latent period (days) [10]	2	2–5
${\mu }_{T_c}^{-1}$	Average contactless period of complicated (days)	49	30–130
${\mu }_T^{-1}$	Average contactless period of symptomatic (days)	14	7–30
p	The probability of being asymptomatic [9,10,22] (p. 588)	$\frac{7}{20}$	0.2–0.8
$p_A$	The probability of recovery without treatment, asymptomatic [12,13,14,15,16]	$\frac{1}{4}$	0.05–0.57
$p_I$	The probability of recovery without treatment, symptomatic [12,13,14,15,16]	$\frac{1}{4}$	0.05–0.57
${\delta }_I^{-1}$	Time to spontaneous clearance, symptomatic (days) [12,13,14,15,16]	30	5–130
${\delta }_A^{-1}$	Time to spontaneous clearance, asymptomatic (days) [12,13,14,15,16]	30	5–130
$\mathrm{\Lambda }$	Birth rate (days −1)	$\frac{25269}{1522883\times 365}$
d	Death rate (days −1)	$\frac{25269}{1522883\times 365}$

. More specifically, people can go through the following disjoint groups. Compartment S contains the individuals who are susceptible to the disease. In compartment E, we have the infected individuals who are not infectious yet but become gonorrhoea spreaders, eventually. The infection is spread by symptomatic and asymptomatic individuals in compartments I and A, respectively. Furthermore, T contains the treated individuals who tested positive primarily because of gonorrhoea symptoms. Finally, individuals who are diagnosed with gonorrhoea while being tested because of signs and symptoms of any other disease warranting testing for STIs, such as PID, are in compartment $T_c$. We refer to $T_c$ as the compartment of individuals treated with complications. We introduce these two compartments assuming that individuals after a positive test avoid having sex until they and their partner have been treated and given the all-clear. (By assuming that the majority of $T_c$ is a result of asymptomatic cases, we do not consider transition from I to $T_c$.) The assumed progressions between compartments are as follows. We assume that the individuals from S move to E with rates of either $\frac{\beta I}{N}$ or $\frac{\alpha \beta I}{N}$ after having some sexual contacts with symptomatic or asymptomatic individuals, respectively. Individuals in E become spreaders, either symptomatic or asymptomatic, after ${\gamma }^{-1}$ units of time. That is, individuals from E either remain asymptomatic with a probability of p and move into A with a rate of $p\gamma$ or develop symptoms with a probability of $1-p$ and move with a rate of $(1-p)\gamma$ into I. The probabilities of spontaneous recovery in I and A are modelled by $p_I$ and $p_A$, respectively; and individuals from these compartments move back to S with rates $p_I{\delta }_I$ and $p_A{\delta }_A$, respectively. Individuals from I and A, with rates $(1-p_I){\sigma }_I$ and $(1-p_A){\sigma }_A$, move to T and $T_c$, respectively. Finally, after receiving treatments, individuals from T and $T_c$ move back to S with rates ${\mu }_T$ and ${\mu }_{T_c}$, respectively.

Using the transmission diagram, we derive the following set of differential equations to model the dynamics of gonorrhoea transmission, considered for non-negative values of the compartments.

$$\begin{array}{cc}\hfill \dot{S}& =\mathrm{\Lambda }N-\frac{\beta S}{N}\left(I+\alpha A\right)+p_A{\delta }_{A}A+p_I{\delta }_{I}I+{\mu }_{T}T+{\mu }_{T_c}{T}_c-dS\hfill \\ \hfill \dot{E}& =\frac{\beta S}{N}\left(I+\alpha A\right)-\gamma E-dE\hfill \\ \hfill \dot{A}& =p\gamma E-(1-p_A){\sigma }_{A}A-p_A{\delta }_{A}A-dA\hfill \\ \hfill \dot{I}& =(1-p)\gamma E-(1-p_I){\sigma }_{I}I-p_I{\delta }_{I}I-dI\hfill \\ \hfill {\dot{T}}_c& =(1-p_A){\sigma }_{A}A-{\mu }_{T_c}{T}_c-d{T}_c\hfill \\ \hfill \dot{T}& =(1-p_I){\sigma }_{I}I-{\mu }_{T}T-dT\hfill \\ \hfill \dot{C}& =(1-p_A){\sigma }_{A}A+(1-p_I){\sigma }_{I}I\hfill \end{array}$$

(1)

The additional variable C does not affect the disease dynamics but it is used in the parameter fitting process, and its role is explained in Section 2.4. With any non-negative initial condition, the compartments in (1) remain non-negative. For instance, if there is a $t_0>0$ such that $S\left(t_0\right)=0$, then $\dot{S}\left(t_0\right)=\mathrm{\Lambda }N+p_A{\delta }_{A}A+p_I{\delta }_{I}I+{\mu }_{T}T+{\mu }_{T_c}{T}_c\ge 0$. Similar arguments apply for the other compartments.

Figure 2. Transmission diagram. The compartments are for susceptible, S, exposed, E, infectious with symptoms, I, infectious without symptoms, A, treated without complications, T, and treated with complications, $T_c$, individuals. The interpretation of the non-negative parameters is given Table 1.

Figure 2. Transmission diagram. The compartments are for susceptible, S, exposed, E, infectious with symptoms, I, infectious without symptoms, A, treated without complications, T, and treated with complications, $T_c$, individuals. The interpretation of the non-negative parameters is given Table 1.

2.3. The Basic Reproduction Number

The basic reproduction number, ${\mathcal{R}}_0$, providing the average number of new infections generated by a typical infected individual in a completely susceptible population, is a dimensionless key epidemiological parameter. That is, if ${\mathcal{R}}_0>1$, the presence of an infectious individual leads to an outbreak. However, ${\mathcal{R}}_0<1$ does not necessarily mean the eradication of a well-established communicable disease [55,56]. To find the ${\mathcal{R}}_0$ of gonorrhoea in NI in terms of the parameters in our model, we used a general method developed in [57].

Existing mathematical models, using different approaches, modelling assumptions and population data, provide a range of estimates of ${\mathcal{R}}_0$. For instance, Ref. [58] provides three different estimates as $1.34$, $2.84$ and $3.81$. Using core groups, in [59], the estimates are $0.77$ and $1.01$ for individuals testing once and more than once, respectively. Furthermore, Ref. [49] finds that while modelling the dynamics of cefixime-susceptible and cefixime-resistant strains, that ${\mathcal{R}}_0$ is approximately $1.19$ for the cefixime-susceptible strain, while it is between $0.73$ and $1.61$, depending on the frequency of cefixime prescription, for the cefixime-resistant strains.

2.4. Parameter Estimation

Since the accessible data sets contains monthly positive tests and the transmission model does not have an explicit equation for the number of positive tests, we use the cumulative number of positive tests at time t, $C_D\left(t\right)$, and proceeded as follows. Since we assumed that the new positive tests are registered only in either in compartment T or $T_c$, the number of new positive tests at any given time t in our model is $(1-p_A){\sigma }_{A}A\left(t\right)+(1-p_I){\sigma }_{I}I\left(t\right)$. Therefore, the cumulative case number at any time time t is

$$C\left(t\right)=C\left(0\right)+{\int }_0^t\left((1-p_A){\sigma }_{A}A\left(s\right)+(1-p_I){\sigma }_{I}I\left(s\right)\right)ds.$$

(2)

By taking the derivative with respect to time on both sides, we obtain the last equation in (1).

We fitted our model to the data by leveraging the Matlab function lsqcurvefit where the initial point in the argument of lsqcurvefit had two parts. The first part was composed of the values given in Table 1, and the upper and lowers bounds of the parameters were set to values listed in the column titled Range. Furthermore, since we had information only about the monthly positive tests, the initial values of the compartments in (1) were considered as parameters with the values listed in

Table 2. Initial values. The values of the compartments (1) at $t=0$ to initiate the parameter fitting process in Matlab by using the function lsqcurvefit where $C_D\left(0\right)=30$.

Compartment	Meaning	Value	Range
S	Susceptible individuals	$N-40-C$
E	Exposed individuals	10	5–20
A	Asymptomatically infectious individuals	10	10–100
I	Infectious individuals with symptoms	10	1–100
$T_c$	Treated individuals complications	10	$0-C_D\left(0\right)$
T	Treated individuals without complications	$C-T_c$	$0-C_D\left(0\right)$
C	Cumulative case numbers	$C_D\left(0\right)$

. That is, the parameter vector in lsqcurvefit is of dimension 17. During the fitting process, we generated the numerical solution of (1) by using the Matlab function ode23 and returned the numerical of values of C to lsqcurvefit.

2.5. Confidence Intervals

We also formed $95\%$ confidence intervals of the estimated values. These intervals were created by using a parametric bootstrap method described in [60]. Namely, from $C\left(t\right)$, we derived the fitted monthly incidence of infection, ${\left\{{\mathrm{\Xi }}_j\right\}}_{j=1}^{129}$, and generated a new random sequence of daily incidence of infection, ${\left\{{\delta }_j\right\}}_{j=1}^{129}$, from the Poisson distribution specified by the rate parameter ${\mathrm{\Xi }}_j,j=1,\cdots ,129$. Then, we fitted C of (1) to $\Delta$, the cumulative values of ${\left\{{\delta }_j\right\}}_{j=1}^{129}$, by applying the steps detailed in Section 2.4. By repeating these steps 10 times, we obtained estimates of the parameters and the initial values of (1) which we used in constructing the corresponding confidence intervals. More specifically, because of the sample size and the unknown standard deviation of the sample, we used the inverse cumulative distribution function of the Student’s t distribution.

2.6. Parameter Sensitivity Analysis

Since most of the expressions presented in Section 3 are rather complicated, we computed the corresponding partial rank correlation coefficients of the parameters contributing to the formula of interest to obtain a better understanding and some interpretations of them. More specifically, by using Latin Hypercube Sampling, [61], we generated parameter sets of 30,000 samples from the considered parameter ranges. Using these parameter values, we obtained 30,000 different values of the expression in questions. Finally, employing pcc from the software R [62], we obtained Partial Rank Correlation Coefficients (PRCC) to assess the effects of the parameters on the given formula.

3. Results

3.1. Fitted Parameters and Initial Values

In Figure 3, we plotted the values of C of the numerical solution of (1), with the fitted parameter values and initial values discussed below, together with the data, D.

3.1.1. Parameters

In Table 3, we provide the fitted parameter values with their associated confidence intervals. As can be seen, the time from exposure to being infectious is approximately $2.6$ days with approximately half of the infectious population being asymptomatic. Furthermore, the spreading period (the time from being infectious to getting a positive test), on average is approximately 33 days for an asymptomatic individual, and it is 20 days for a symptomatic one. Furthermore, the probability of spontaneous recovery is $0.42$ after 14 and 9 days in the asymptomatic and symptomatic groups, respectively. Finally, following a positive test, symptomatic and asymptomatic individuals re-enter the susceptible population after about 30 and 128 days, respectively.

3.1.2. Initial Values

Because the available epidemiological data contain only the number of positive tests, we also estimated the initial values of compartments $E,A,I,T_c$ and T during the fitting process, and their approximations are summarised in

Table 4. Approximated initial values. The values are obtained in Matlab by using the function lsqcurvefit. The confidence intervals are constructed as described in Section 2.5.

Comparment	Meaning	Value	CI
E	Exposed individuals	$11.8904$	11.3232–12.1985
A	Asymptomatically infectious individuals	$35.6136$	28.3610–36.4412
I	Infectious individuals with symptoms	$23.4546$	21.9749–29.3474
$T_c$	Treated individuals complications	$3.0066$	0.1901–8.1531
T	Treated individuals without complications	$26.9934$	Not computed

. Namely, we estimate that the numbers of exposed, infectious without symptoms, infectious with symptoms, treated with complications, and treated primarily because of gonorrhoea symptoms individuals were 12, 37, 23, 3 and 27, respectively, in January 2012.

Table 3. Fitted parameter values. The parameter values are obtained in Matlab by using the function lsqcurvefit where $C_D\left(0\right)=30$. The confidence intervals are constructed as described in Section 2.5.

Parameter	Meaning	Value	CI
$\beta$	Transmission rate	$0.1007$	0.0968–0.1020
$\alpha$	Infectiousness reduction in asymptomatic	$0.3656$	0.3562–0.3927
${\sigma }_A^{-1}$	Time from infectiousness to positive test of asymptomatic (days)	$32.7496$	30.4953–37.4233
${\sigma }_I^{-1}$	Time from infectiousness to positive test of symptomatic (days)	$19.3535$	16.4988–19.6343
${\gamma }^{-1}$	Average latent period (days)	$2.5620$	2.4638–2.6256
${\mu }_{T_c}^{-1}$	Average contactless period of complicated (days)	$127.8659$	48.5581–190.5469
${\mu }_T^{-1}$	Average contactless period of symptomatic (days)	$29.9861$	11.8992–39.4454
p	The probability of being asymptomatic	0.5651	$0.5512-0.5755$
$p_A$	The probability of recovery without treatment, asymptomatic	$0.4132$	0.3959–0.4165
$p_I$	The probability of recovery without treatment, symptomatic	$0.4191$	0.3989–0.4213
${\delta }_A^{-1}$	Time spontaneous clearance, asymptomatic	$14.4923$	14.3744–15.1723
${\delta }_I^{-1}$	Time spontaneous clearance, symptomatic	$8.5805$	8.4850–8.8217

Table 3. Fitted parameter values. The parameter values are obtained in Matlab by using the function lsqcurvefit where $C_D\left(0\right)=30$. The confidence intervals are constructed as described in Section 2.5.

Parameter	Meaning	Value	CI
$\beta$	Transmission rate	$0.1007$	0.0968–0.1020
$\alpha$	Infectiousness reduction in asymptomatic	$0.3656$	0.3562–0.3927
${\sigma }_A^{-1}$	Time from infectiousness to positive test of asymptomatic (days)	$32.7496$	30.4953–37.4233
${\sigma }_I^{-1}$	Time from infectiousness to positive test of symptomatic (days)	$19.3535$	16.4988–19.6343
${\gamma }^{-1}$	Average latent period (days)	$2.5620$	2.4638–2.6256
${\mu }_{T_c}^{-1}$	Average contactless period of complicated (days)	$127.8659$	48.5581–190.5469
${\mu }_T^{-1}$	Average contactless period of symptomatic (days)	$29.9861$	11.8992–39.4454
p	The probability of being asymptomatic	0.5651	$0.5512-0.5755$
$p_A$	The probability of recovery without treatment, asymptomatic	$0.4132$	0.3959–0.4165
$p_I$	The probability of recovery without treatment, symptomatic	$0.4191$	0.3989–0.4213
${\delta }_A^{-1}$	Time spontaneous clearance, asymptomatic	$14.4923$	14.3744–15.1723
${\delta }_I^{-1}$	Time spontaneous clearance, symptomatic	$8.5805$	8.4850–8.8217

3.2. Simulations

In Figure 4, we visualised the numerical solution of the disease-spread modelling compartments of (1) using the fitted parameters and the initial conditions in Table 3 and Table 4. In Figure 5, we visualised the estimated monthly number of new cases obtained from the values of C from the numerical solution of (1). Together with that, we also plotted the data set, D, and its trend component obtained in Matlab by using its function trenddecomp. As can be seen from the figure, the simulated numerical values of the monthly new cases follow closely the trend of the data set. That is, our model does not capture the oscillation, probably due to seasonal patterns, of the data. We will mention some possible, potentially oscillation-inducing, extensions of the model in Section 4. Furthermore, the data show a big drop between from the beginning of 2020 until early 2021 followed by a sharp rise, which potentially might be due to COVID-19-related changes in human behaviour and the health care system [63].

The bands around the numerical solutions in Figure 4 and Figure 5 represent the uncertainty in modelling and testing. To generate the bands around the solutions—by using Latin Hypercube Sampling [61]—we prepared 5000 samples of the parameters from the confidence intervals given in Table 3 and generated the numerical solutions of (1) with the initial values in Table 4. Out of those solutions, we kept the 239 solutions which satisfied

$$D\left(t\right)\left(1-\frac{P}{100}\right)\le C\left(t\right)\le D\left(t\right)\left(1+\frac{P}{100}\right),t=129$$

(3)

where $P=20$ in this study, that is, we assume $20\%$ overall uncertainty.

Figure 4. Numerical solutions. The 30-day evaluations of the numerical solutions of compartments $S,E,A,I,T_c\mathrm{and}T$ of (1).

Figure 4. Numerical solutions. The 30-day evaluations of the numerical solutions of compartments $S,E,A,I,T_c\mathrm{and}T$ of (1).

Figure 5. The monthly number of new positive tests. The blue dots are the monthly new cases from the data, D. The green dots are from Matlab’s trenddecomp of the data in blue. The red dots are the monthly new cases from the numerical solution C of (1) using the parameters and the initial conditions in Table 3 and Table 4.

Figure 5. The monthly number of new positive tests. The blue dots are the monthly new cases from the data, D. The green dots are from Matlab’s trenddecomp of the data in blue. The red dots are the monthly new cases from the numerical solution C of (1) using the parameters and the initial conditions in Table 3 and Table 4.

3.3. The Basic Reproduction Number

In the absence of the infectious agent, the compartments in (1) settle on the so-called disease-free equilibrium,

$$DFE=(N,0,0,0,0,0)$$

in our model, which plays an important role in finding ${\mathcal{R}}_0$. (Notice that we do not use C as it has no impact on the disease dynamics.) Namely, using the method from [57] and assuming that $\mathrm{\Lambda }=d$, we have

$${\mathcal{R}}_0=\underset{{\mathcal{R}}_{0,A}}{\underbrace{\frac{\alpha \beta \gamma }{(\gamma +d)}\frac{p}{d+p_A{\delta }_A+{\sigma }_A(1-p_A)}}}+\underset{{\mathcal{R}}_{0,I}}{\underbrace{\frac{\beta \gamma }{(\gamma +d)}\frac{1-p}{d+p_I{\delta }_I+{\sigma }_I(1-p_I)}}}.$$

(4)

The derivation of the formula can be found in Section 6.1. As the formula indicates, the average of the total number of secondary cases is generated via the routes of asymptomatic and symptomatic infections. The corresponding contributions are denoted by ${\mathcal{R}}_{0,A}$ and ${\mathcal{R}}_{0,I}$, respectively.

Using the values from Table 3 in (4), we have the estimates

$${\mathcal{R}}_{0,A}\approx 0.4477\mathrm{and}{\mathcal{R}}_{0,I}\approx 0.5552,$$

that is,

$${\mathcal{R}}_0\approx 1.0030.$$

Furthermore, we included the basic reproduction number in the process of constructing the confidence intervals by evaluating (4), and we obtained the corresponding confidence interval for ${\mathcal{R}}_0$ as $[1.0029,1.0032]$.

While generating the bands in Figure 4, we also evaluated (4) with the kept 239 parameters to obtain

$$0.9996\le {\mathcal{R}}_0\le 1.0056.$$

(5)

Notice that although the estimates in (5) are obtained by using parameter values from the confidence intervals of the parameters, they provide a wider range than its confidence interval for the possible values of ${\mathcal{R}}_0$. However, it is natural since the parameters were selected randomly from the confidence intervals, and parameter combinations can certainly provide evaluations of ${\mathcal{R}}_0$ outside of the confidence interval provided above.

In Figure 6, we visualised the results of the parameter sensitivity analysis of ${\mathcal{R}}_0$. As it can be seen, after the transmission rate $\beta$, ${\sigma }_A$ and ${\sigma }_I$ are the most influential contributors to ${\mathcal{R}}_0$. Furthermore, since these partial rank correlation coefficients are negative, their increasing value decrease the basic reproduction number. That is, by using the interpretation of these parameters, it means that decreasing the time between becoming infectious and testing positive in the groups of the symptomatic or asymptomatic individuals implies a reduction in ${\mathcal{R}}_0$. Notice that the reduction in the other spreading times, ${\delta }_A^{-1}$ and ${\delta }_I^{-1}$, also have ${\mathcal{R}}_0$ lowering impacts. However, since the contributions of these are smaller compared to ${\sigma }_A$ and ${\sigma }_I$, we do not investigate their roles further.

To investigate further the way in which $\beta$, ${\sigma }_A^{-1}$ and ${\sigma }_I^{-1}$ affect ${\mathcal{R}}_0$, we visualised the level curves of ${\mathcal{R}}_0$ in Figure 7. In these plots, only the variables on the axes were allowed to change while the other parameters were set to the fitted values provided in Table 3. These plots indicate that for any fixed value of $\beta$, a reduction in time between becoming infectious and being tested positive has a basic reproduction number reducing effect. We illustrate these observations further in Section 3.4 and discuss their implications in Section 4.

3.4. The Endemic Equilibrium

If there is an outbreak—such that ${\mathcal{R}}_0>1$—with no intervention, the population will approach the so-called endemic equilibrium, $EE$, as time goes by.

$$\begin{array}{cc}\hfill S_{\infty }& =N\frac{1}{{\mathcal{R}}_0}\hfill \\ \hfill E_{\infty }& =N\frac{{\mathcal{R}}_0-1}{{\mathcal{R}}_0}\frac{1}{1+{\mathcal{R}}_{0,A}\frac{d+\gamma }{\alpha \beta }\frac{d+{\mu }_C+{\sigma }_A(1-p_A)}{d+{\mu }_C}+{\mathcal{R}}_{0,I}\frac{d+\gamma }{\beta }\frac{d+{\mu }_T+{\sigma }_I(1-p_I)}{d+{\mu }_T}}\hfill \\ \hfill A_{\infty }& =E_{\infty }\frac{\gamma +d}{\alpha \beta }{\mathcal{R}}_{0,A}\hfill \\ \hfill I_{\infty }& =E_{\infty }\frac{\gamma +d}{\beta }{\mathcal{R}}_{0,I}\hfill \\ \hfill {T_c}_{\infty }& =E_{\infty }\frac{(\gamma +d){\sigma }_A(1-p_A)}{\alpha \beta (d+{\mu }_C)}{\mathcal{R}}_{0,A}\hfill \\ \hfill T_{\infty }& =E_{\infty }\frac{(\gamma +d){\sigma }_I(1-p_I)}{\beta (d+{\mu }_T)}{\mathcal{R}}_{0,I}\hfill \end{array}$$

(6)

The derivation of these formulas can be found in Section 6.2.

Figure 6. Parameter sensitivity analysis ${\mathcal{R}}_0$ in (4). The partial rank correlation coefficients of ${\mathcal{R}}_0$ in (4) is computed by using the method described in Section 2.6 with 30,000 random parameter samples from the ranges in Table 1. Positive means increasing parameter value is increasing ${\mathcal{R}}_0$, and negative means increasing parameter value is decreasing ${\mathcal{R}}_0$.

Figure 6. Parameter sensitivity analysis ${\mathcal{R}}_0$ in (4). The partial rank correlation coefficients of ${\mathcal{R}}_0$ in (4) is computed by using the method described in Section 2.6 with 30,000 random parameter samples from the ranges in Table 1. Positive means increasing parameter value is increasing ${\mathcal{R}}_0$, and negative means increasing parameter value is decreasing ${\mathcal{R}}_0$.

Figure 7. The contour plots of ${\mathcal{R}}_0$ in (4). Panel (a), $\beta$ and ${\sigma }_I^{-1}$ changes in the ranges given in Table 1, while the other parameters are kept at the fitted values as given in Table 3. Panel (b), $\beta$ and ${\sigma }_A^{-1}$ changes in the ranges given in Table 1 while the other parameters are kept at the fitted values as given in Table 3.

Figure 7. The contour plots of ${\mathcal{R}}_0$ in (4). Panel (a), $\beta$ and ${\sigma }_I^{-1}$ changes in the ranges given in Table 1, while the other parameters are kept at the fitted values as given in Table 3. Panel (b), $\beta$ and ${\sigma }_A^{-1}$ changes in the ranges given in Table 1 while the other parameters are kept at the fitted values as given in Table 3.

Notice that, $E_{\infty }>0$ if and only if ${\mathcal{R}}_0>1$. Furthermore, from the formulas in (6), it is clear that the same holds for the endemic values of the other compartments except the two treatment compartments. Namely, if ${\sigma }_A=0$ or ${\sigma }_I=0$, that is, asymptomatic or symptomatic individuals are not tested, the above formulas are still valid, and ${T_c}_{\infty }=0$ or $T_{\infty }=0$. Furthermore, these assumptions have an ${\mathcal{R}}_{0,A}$ or ${\mathcal{R}}_{0,I}$ increasing effect, respectively, and they also increase ${\mathcal{R}}_0$, as can be seen from (4).

Substituting the fitted parameter values into the formulas above, we obtain

$$EE\approx \left(151835421910414722372424\right).$$

(7)

To visualise how our model reaches this state, we plotted the numerical solutions of the infected compartments on a time interval necessarily long enough in Figure 8 by using the estimates from Table 3 and Table 4. The visualised simulation is based on the, very unlikely, hypothetical scenario of unchanging circumstances over a very long period of time. Nevertheless, besides providing the same endemic values as (7), it indicates a very slow growth of the infected compartments.

From a disease management viewpoint, in our model, it is only compartments $T_{\infty }$ and ${T_c}_{\infty }$ which are visible to and imply direct costs for the health care system, hence we investigated the impacts of the parameters on $T_{\infty }$, ${T_c}_{\infty }$ and $T_{\infty }+{T_c}_{\infty }$. (Although the cost of an STI test is not negligible, it is not accessible in our model; however, by incorporating tested compartments—together with the up to two-week waiting time [11]—simultaneously in the asymptomatic and symptomatic routes, one can investigate the effects of testing for the disease.) To investigate the impacts of changing parameters on $T_{\infty }$, ${T_c}_{\infty }$ and $T_{\infty }+{T_c}_{\infty }$, we computed the partial rank correlation coefficients of the parameters in (6) for these groups and visualised them in Figure 9. Similarly to the sensitivity analysis of ${\mathcal{R}}_0$ in Section 3.3, parameters affecting the transmission rates, $\beta$ and $\alpha$, have a large positive impact on the endemic values in question. Furthermore, when $T_{\infty }$ and ${T_c}_{\infty }$ are considered separately, the probability of being asymptomatic, p, has the largest impact on their values, and, as expected from the model derivation in Section 2.2, an increase in p decreases $T_{\infty }$ and increases ${T_c}_{\infty }$. However, the impact of p reduces significantly when the total number of treated individuals is considered in the endemic state. In addition, the probability of being asymptomatic may correlate with human behaviour through multiple exposure; however, this aspect of gonorrhoea is not part of our model. Furthermore, the transfer rates ${\mu }_T$ and ${\mu }_{T_c}$ show significant contributions to the changes in the endemic values. However, their roles, for instance the $T_{\infty }$ reducing contactless period from the group of treated individuals with symptoms, is somewhat expected on the one hand. On the other hand, it would be unrealistic to expect from a health care system to force patients out of treated status, in order to reduce the level of reported endemic states. Furthermore, the spontaneous recovery times, ${\delta }_A^{-1}$ and ${\delta }_I^{-1}$, show noticeable contributions to the endemic values of the treatment compartments. Similarly to ${\mu }_T^{-1}$ and ${\mu }_{T_c}^{-1}$, these times can be shortened via various human behaviour influencing interventions as mentioned in Section 4. However, a good understanding, even if only based on simulations, of the interplays between these parameters would require a systematic investigation, and it would add significantly to the length of this study. Therefore, in what follows, we do not consider further how p, ${\mu }_T$, ${\mu }_{T_c}$, ${\delta }_A^{-1}$ and ${\delta }_I^{-1}$ affect the endemic values. Instead, we observe that, similarly to the sensitivity analysis of ${\mathcal{R}}_0$ in Section 3.3, parameters ${\sigma }_A$ and ${\sigma }_I$, in addition to the transmission rates, have a somewhat noticeable impact on $T_{\infty }$, ${T_c}_{\infty }$ and $T_{\infty }+{T_c}_{\infty }$. Furthermore, they can be relevant parameters from an outbreak (endemic in this case) control viewpoint. In addition, as discussed below, their roles in the endemic values are not as straightforward as in the values of the basic reproduction number, even if the value of $\beta$ is frozen. For instance, Figure 9 suggests that a decreasing spreading time, whether symptomatic or not, reduces $T_{\infty }$; although, the impact of ${{\sigma }_I}^{-1}$ is relatively small. On the other hand, while decreasing ${{\sigma }_I}^{-1}$ reduces, decreasing ${{\sigma }_A}^{-1}$ increases ${T_c}_{\infty }$, with a relatively small impact of the latter. As for $T_{\infty }+{T_c}_{\infty }$, it is only ${{\sigma }_I}^{-1}$ with a considerably significant impact on it, with the response of a decrease from the total number of treated individuals to decreasing values of ${{\sigma }_I}^{-1}$.

To investigate effects of the parameters $\beta ,{\sigma }_A^{-1}$ and ${\sigma }_I^{-1}$ on $T_{\infty }$, ${T_c}_{\infty }$ and $T_{\infty }+{T_c}_{\infty }$, the surfaces defined by them are plotted in Figure 10 when the other parameters are kept at the fitted values provided in Table 3. These plots indicate some non-trivial responses of the sizes of the treatment compartments to parameter changes. More specifically, they show both monotonic and non-monotonic responses of $T_{\infty }$, ${T_c}_{\infty }$ and $T_{\infty }+{T_c}_{\infty }$ to changes in ${\sigma }_A^{-1}$ and ${\sigma }_I^{-1}$ as visible from their equally spaced evaluations, in black, for some fixed values of $\beta$.

To gain a better understanding of the roles that $\beta ,{\sigma }_A^{-1}$ and ${\sigma }_I^{-1}$ play in the endemic values of the treated compartments, the level curves of $T_{\infty }$, ${T_c}_{\infty }$ and $T_{\infty }+{T_c}_{\infty }$ are shown in Figure 11. In the observations below, $\beta$ is considered as fixed at some value, and the endemic value in question as a function of one variable. As these plots reveal, the impact of $\beta ,{\sigma }_A^{-1}$ and ${\sigma }_I^{-1}$, on the considered endemic values is indeed more complex than their monotonic effects on ${\mathcal{R}}_0$. A monotonic, endemic size-decreasing impact of the reduced time from being infectious to a positive test is only observable with respect to ${\sigma }_A^{-1}$ on $T_{\infty }$, and with respect to ${\sigma }_I^{-1}$ on ${T_c}_{\infty }$ and $T_{\infty }+{T_c}_{\infty }$, for any fixed value of $\beta$. These observations are in line with PRCC-based findings. As for the other plots in Figure 11, it seems that each might be related to the corresponding PRCC plots by taking the ratios between the areas on which the endemic value in question is increasing and decreasing. However, we do not discuss this issue further here. Instead, we note that, after inspecting the contour plots in Figure 11, it seems that the benefits of reducing either ${\sigma }_I^{-1}$ or ${\sigma }_A^{-1}$ are present in each of the configurations in terms of lowering the number of treated patients, at least for small values of $\beta$. However, these plots also indicate that in each case, there is a critical value of $\beta$ when these benefits are lost if the transmission rate is fixed at a value greater than the critical value. For example, the $\beta -{\sigma }_I^{-1}$ plot of $T_{\infty }$ suggests the existence of a ${\beta }_0\in (0.3,0.4)$ such that $T_{\infty }$ is a monotonically decreasing function of ${\sigma }_I^{-1}$ for any fixed $\beta >{\beta }_0$. Furthermore, in the plot of ${T_c}_{\infty }$ against $\beta$ and ${\sigma }_A^{-1}$, it seems that for any $0.15<\beta <0.3$, the vertical line at that value intersects most of the level curves at least twice. Similar observations can be made about $T_{\infty }\left({\sigma }_I^{-1}\right)$, ${T_c}_{\infty }\left({\sigma }_A^{-1}\right)$ and the total number of treated individuals with respect to ${\sigma }_A^{-1}$. That is, for some values of the transmission rate $\beta$, the observations indicate that although a reduction in time between being infectious and testing positive initially may elevate the number of treated individuals, a further decrease in the time interval in question may lower those numbers eventually. For instance, by fixing $\overline{\beta }\approx 0.21$ in the $\beta -{\sigma }_I^{-1}$ level curves of $T_{\infty }$ and taking the estimated value $\widehat{{\sigma }_I^{-1}}\approx 20$, we can see that $\mathrm{60,000}\le T_{\infty }\le \mathrm{80,000}$ which changes to $\mathrm{80,000}\le T_{\infty }\le \mathrm{100,000}$ when ${\sigma }_I^{-1}\approx 14$, and back to $\mathrm{60,000}\le T_{\infty }\le \mathrm{80,000}$ when ${\sigma }_I^{-1}\approx 4$. However, it is not clear from this plot whether there is a critical value of $\widetilde{{\sigma }_I^{-1}}$ such that $T_{\infty }\left({\sigma }_I^{-1}\right)<T_{\infty }\left(\widehat{{\sigma }_I^{-1}}\right)$ at least for some $0<{\sigma }_I^{-1}<\widetilde{{\sigma }_I^{-1}}$.

To investigate the effects $\beta ,{{\sigma }_A}^{-1}$ and ${{\sigma }_I}^{-1}$ on the endemic values of treated compartments further, in Figure 12, we re-plotted the level curves of parameter configurations, together with functions of one variable for the indicated fixed values of the transmission rate, when the phenomenon mentioned above is potentially present in Figure 11. As these plots reveal, in each of the investigated cases, there is a ${{\sigma }_A}_0^{-1}$ or ${{\sigma }_I}_0^{-1}$ such that the endemic value in question is below its evaluation at the fitted values of the spreading times in Table 3, for any ${{\sigma }_A}^{-1}<{{\sigma }_A}_0^{-1}$ or ${{\sigma }_I}^{-1}<{{\sigma }_I}_0^{-1}$, for any fixed $\beta$. However, in most of the cases, a significantly shortened period between the start of infectivity and a positive test is required for such an impact to manifest. For instance, as panel (d) suggests, ${{\sigma }_A}^{-1}$ should be reduced from 33 to 12 days, when $\beta \approx 0.185$. On the other hand, for the same transmission rate, a reduction in ${{\sigma }_I}^{-1}$ from 20 to 5 days would lead to a slightly lower value of $T_{\infty }$ compared to its level at the fitted parameter values, as panel (b) suggests. Furthermore, on panels (b), (d) and (f), we can observe that the endemic values of the treatment compartments have maximum values which increase by the increasing values of $\beta$. Furthermore, the larger the value of $\beta$, the smaller the value of the independent variable, the spreading time, at which the maximum value is taken. For instance, in panel (f), while the maximum value of $T_{\infty }+{T_c}_{\infty }$ is $4.3453e+05$ at ${{\sigma }_A}^{-1}\approx 25.2462$ when $\beta \approx 0.1850$, it is $5.6714e+05$ at ${{\sigma }_A}^{-1}\approx 16.9880$ for $\beta \approx 0.2350$. During the preparation of Figure 12, we also computed the estimates of the components of the basic reproduction number which are provided in the caption of Figure 12. As can be seen, the values of ${\mathcal{R}}_0$ are above 1 except in the case of $T_{\infty }$ when its smallest value is approximately $0.95$. Furthermore, the small white region in the bottom left corner of panel (a) marks the set $M=\{(\beta ,{\sigma }_I):{\mathcal{R}}_0<1\}$. This implies that decreasing values of ${\sigma }_I$ are not only reducing $T_{\infty }$ but also lower ${\mathcal{R}}_0$ below 1 eventually, for any fixed $\beta \in M$.

Figure 11. The level sets of the treated compartments. The contour plots of $T_{\infty }$ (panels (a,b)),${T_c}_{\infty }$ (panels (c,d))$\mathrm{and}{T}_{\infty }+{T_c}_{\infty }$ (panels (e,f)) using (6).

Figure 11. The level sets of the treated compartments. The contour plots of $T_{\infty }$ (panels (a,b)),${T_c}_{\infty }$ (panels (c,d))$\mathrm{and}{T}_{\infty }+{T_c}_{\infty }$ (panels (e,f)) using (6).

The observations above based on numerical simulations, and their rigorous verification is outside of the scope of this paper. Nevertheless, we discuss some aspects of these observations in relation to epidemic control in Section 4.

4. Discussion

We start our discussion by restricting our attention to the treatment compartments, since, in our model, they are the only ones that may be of concern to any health care system. Although the system would evolve very slowly to its endemic state, as Figure 8 suggests, the constant presence of the relatively large number treatment-requiring individuals might warrant some interventions to lower the endemic values of the infected compartments. Hence, we followed steps similar to the investigations of the same intention regarding ${\mathcal{R}}_0$. The sensitivity analysis of the endemic values showed some, to a degree expected, dependence on parameters. For instance, when considered separately, the number of treated individuals without complications would decrease as a response to an increased transfer rate from that group. A less intuitive observation is that the size of the same compartment would be increased by people spending less time in the other, disjointed, treatment compartment. An analogous observation, by using the reversed roles of parameters, can be formulated about the number of individuals treated with complications. However, these effects diminish when the total number of treated individuals is considered. These observations require further analysis, which is not in the scope of this paper. Furthermore, it is unlikely that the treatment period could be shortened below a necessary value in order to lower the sizes of the groups of treated individuals. Therefore, to assess the possibility and the effects of intervention strategies, we focused on the ways in which the transfer rate, $\beta$, and the timing of positive tests influence the size of those groups. The level curves of values of the treatment compartments, projected on the corresponding parameter planes, revealed some nontrivial dependence on parameters of the groups of treated individuals. Namely, for certain fixed values of the transmission rate, we found that the period from being infectious to a positive test in certain situations changes the number of the individuals in the groups of treated individuals non-monotonically. More specifically, the reduced period initially increases the sizes of some of the compartments before lowering those values. Furthermore, in most of the situations, a significantly earlier test is needed to achieve the desired results. However, this phenomenon is present only for large values of the transmission rate, and because the fitted value of $\beta$ is seemingly not in those regions, we expect a monotonic decrease in the epidemic sizes of the treated compartments when the fitted parameter values are used in simulations. To illustrate this effect, we ran simulations using the fitted parameter values and initial conditions, listed in Table 3 and Table 4, respectively. As the numerical solutions, plotted in panels (a) and (b) of Figure 13, indicate, by shortening the period from being infectious to a positive test even only by one day, either in the asymptomatic or in the symptomatic group, the system is moved away from evolving into its endemic state. Furthermore, if the shortening of periods is achieved simultaneously in both groups, panel (c) of Figure 13, not only are the peak values of the sizes of the treated groups lowered, but the time needed to reach levels similar to the ones shown in panels (a) and (b) is also shortened. These findings are in line with the conclusion of [64] claiming that it is unlikely for gonorrhoea to evolve into and maintain an endemic state in a heterosexual population even with the addition of anal and oral sex under the assumption of prompt treatment of symptomatic infections. In panels (d)–(f) of Figure 13, we visualised the effects of the reduction in the period from being infectious to a positive test by an additional day. Furthermore, in the caption of Figure 13, we listed the components of the basic reproduction number. These findings suggest that strategies similar to the efforts targeting the basic reproduction number can be implemented not only to lower the endemic values of the treated compartments but also to eradicate gonorrhoea. However, since NG has acquired or developed almost all known physiological mechanisms of antimicrobial resistance to all antimicrobials recommended and/or used for treatment, the statement on its eradication might be far too optimistic [22,28,65]. Notice that our model assumes perfect treatment for gonorrhoea which is probably the most important of its limitations. We used this assumption, because the data we used provided information about positive tests only but not about the results of retests, and we did not model this aspect of the disease. However, treatment failure can be added to our model by considering removed compartments, as detailed below, and routes from them back to the treatment compartments without registering them as new positive tests.

There are several further limitations of our model. The first and foremost is that it considers the homogeneous spread of an STI, and as such, it would require the consideration of structuring the population at different levels and from different aspects such as age, sexual orientation and socioeconomic background, etc. Furthermore, Ref. [44] argues that the site of infection, in particular when the transmission in men who have sex with men is modelled, should be considered in modelling gonorrhoea spread. These limitations might be addressed by incorporating the considerations we have already mentioned in Section 1.2. Nevertheless, our model can also provide insights—as a first, at population, level of assessment—into the trends of the spread of venereal diseases with similar profiles such as Chlamydia or trichomoniasis, both linked to serious complications if untreated. The findings above suggest that the control of gonorrhoea spread might require different strategies in different circumstances, and those findings are in line with traditionally recommended infectious disease impact mitigating protocols. In [66], it is argued that the most effective ways of controlling STIs are (a) reducing rates of new sex partner acquisition, (b) reducing susceptibility of exposed individuals, and (c) reducing duration of infectivity of spreaders. However, in the case of gonorrhoea, because of its features mentioned in Section 1, (a) and (c) are the only seemingly available strategies, and our findings also single these out as potential control activities. To achieve and maintain a desired and significant impact on those disease spread affecting parameters, school-based sexual health education programmes are probably the most effective ones, in particular if they follow the recommendations in [67]. The efficacy of these programmes can probably be reinforced later on via information reiteration at the time of a positive test. However, these strategies might be insufficient when the number of positive tests is raising sharply, or in the endemic state when the implementation of the two main forms of combatting any infectious disease, screening and contact tracing, should be considered [58]. In particular, screening programs act at the individual level by testing and treating a random sample of the population, or a subpopulation of high-risk individuals. In addition, because of the large number of asymptomatic cases, Ref. [9] argues for routine screening at high-volume sites for STIs as the most effective way of controlling them. For instance, in a screening program for gonorrhoea in the USA in the 1970s, the detection of an additional $10\%$ of NG infections in women resulted in a $20\%$ decrease in infections in the years immediately following, an observation which also led to the first description of core groups (a small subpopulation) [68]. Furthermore, the likelihood of false positive tests may lead to harm in the form of stigmatisation, needless use of antibiotics and avoidable financial costs [69]. However, the timing of a positive test is affected not only by the infected individuals but also by the capacity of the health care system. For instance, the detrimental impact on sex infections of staff shortages at testing sites was reported in Britain in 2005 [70,71]. In addition to the already mentioned measures, awareness raising media campaigns might be utilised and be added to our model [72]. In the deterministic setting, these modelling approaches might induce oscillatory behaviour in the solutions of the formulated differential Equations [73,74]. Furthermore, our model does not separate the treatment-time and the after-treatment contactless period. By adding chains of removed compartments, following the treatment compartments, this aspect can be incorporated into our model which in turn might provide a better insight into the post-treatment attitudes of infected individuals. Such an extension might lead to oscillations if the added chains are long enough [75]. On the other hand, by structuring the susceptible group according to their risk of being infected, the model might exhibit a sub-threshold endemic state making the design of control strategies even more challenging [76]. Furthermore, the positive-testing-individual-based partner notification system is considered as an important aspect of controlling STIs, but its modelling requires further consideration.

Without the intention of investigating the possible effects that the COVID-19 pandemic might have on gonorrhoea spread in detail, we also fitted our model to the restricted period between January 2012 and January 2020 of the data. This process resulted in $\beta \approx 0.1001$, ${\gamma }^{-1}\approx 2.0247$, ${\mu }_T^{-1}\approx 7.2451$, ${\mu }_{T_c}^{-1}\approx 30.0452$, $\alpha \approx 0.3756$, $p\approx 0.5749$, ${\sigma }_A^{-1}\approx 38.4660$, ${\sigma }_I^{-1}\approx 15.4248$, $p_I\approx 0.4022$, $p_A\approx 0.4088$, ${\delta }_I\approx 8.5672$ and ${\delta }_A\approx 15.0948$. By comparing these estimates with the ones in Table 3, we can see that, apart from ${\mu }_T^{-1}$ and ${\mu }_{T_c}^{-1}$, the estimated parameter values are relatively close to each other. Furthermore, by using these estimates, we have that ${\mathcal{R}}_{0,A}\approx 0.5087$, ${\mathcal{R}}_{0,I}\approx 0.4962$ and ${\mathcal{R}}_0\approx 1.0049$. The estimated initial values are $S\approx 1522798.088$, $E\approx 11.2773$, $A\approx 20.9249$, $I\approx 27.6042$ and $T_c\approx 15.1047$.

5. Conclusions

In this paper, we derived and investigated, both analytically and numerically, a deterministic, continuous-time compartmental description of gonorrhoea to model the dynamics of its spread by sexual means in Northern Ireland between January 2012 and September 2022. In particular, we considered both the symptomatic and asymptomatic routes of disease transmission. Furthermore, we incorporated spontaneous recovery and treated groups. By fitting our model to the reported number of monthly positive tests, we obtained estimates, with confidence intervals of $95\%$, of the gonorrhoea spread-dynamics affecting parameters. We found that it takes approximately 3 days to become infectious after exposure, and slightly more than half of the infectious population remains asymptomatic. Furthermore, half of the infections both in the groups of asymptomatic and symptomatic individuals remains undetected due spontaneous recovery of symptomatic and asymptomatic individuals after 14 and 8 days, respectively. We also found that for infectious individuals, from the the onset of infectiousness, it takes either around 33 days if they are asymptomatic or around 20 days if they are symptomatic to test positive. In addition, individuals do not re-enter the routes of disease transmission for approximately 30 days if the test was taken because of the presence symptoms gonorrhoea, and 128 days otherwise. We also derived a formula for the basic reproduction number, ${\mathcal{R}}_0$, and using the formula, together with the estimated parameter values, we estimate that ${\mathcal{R}}_0\approx 1.003$, $\mathrm{CI}:$ 1.0029–1.0032.

Since the estimated value of ${\mathcal{R}}_0>1$ indicates a gonorrhoea outbreak, we investigated the possibility of potential interventions to lower the basic reproduction number. First, we carried out a sensitivity analysis of ${\mathcal{R}}_0$ to find parameters of large contributions to its values. We found that parameters with the largest increasing contributors are the transmission rates which might be reduced by sexual education and awareness campaigns with an emphasis on safe sex practises. In addition to the transmission rates, we found the timing of testing as a significant contributor to the value of ${\mathcal{R}}_0$. More specifically, lowering the period between becoming infectious and obtaining a positive test decreases the value of the basic reproduction number. These periods, similarly to the transmission rates, might be influenced through sexual education and awareness campaigns urging people to avail of sexual health testing if they are sexually active or intending to become sexually active soon. By using the counterplots of ${\mathcal{R}}_0$ projected on the corresponding parameter spaces, we observed that the basic reproduction number lowering effects are monotonic in each of the cases when only one them is allowed to change.

We also derived the formulas of the endemic values of the compartments of our model. The formulas show that the endemic equilibrium exists if and only if ${\mathcal{R}}_0>1$. Upon substituting the fitted parameter values into the formulas, we found that in the endemic state, the number of individuals in the susceptible, exposed, infections without symptoms, infectious with symptoms, treated with complications and treated without complications is $1518354,219,1041,472,2372\mathrm{and}424$, respectively. We also found that the these values agree with the corresponding values of the numerical simulation.

In addition to the already discussed possible extensions of our model, we plan to investigate the effects of environmental and socioeconomic factors on future local and global spatial spread of gonorrhoea. To achieve these goals, whenever necessary, stratified (by age, susceptibility, or contact patterns, etc.) populations will be considered. In addition, health economic factors will be integrated into the disease transmission frameworks.

6. The Derivation of (4) and (6)

6.1. Derivation of the Basic Reproduction Number

By following the method and using the notations from [57], we consider the vector $x={(E,I,A,T,T_c,S)}^T$. Furthermore, we define functions

$$\mathcal{F}:{\mathbb{R}}^6\ni {\left(E,I,A,T,T_c,S\right)}^T\mapsto {\left(\frac{\beta S}{N}\left(I+\alpha A\right),0,0,0,0,0\right)}^T\in {\mathbb{R}}^6$$

and

$$\mathcal{V}:{\mathbb{R}}^6\ni {\left(E,I,A,T,T_c,S\right)}^T\mapsto \left(\begin{array}{c}\gamma E+dE\\ -(1-p)\gamma E+(1-p_I){\sigma }_{I}I+p_I{\delta }_{I}I+dI\\ -p\gamma E+(1-p_A){\sigma }_{A}A+p_A{\delta }_{A}A+dA\\ -(1-p_I){\sigma }_{I}I+{\mu }_{T}T+dT\\ -(1-p_A){\sigma }_{A}A+{\mu }_{T_c}{T}_c+d{T}_c\\ -\mathrm{\Lambda }N+\frac{\beta S}{N}\left(I+\alpha A\right)-p_A{\delta }_{A}A-p_I{\delta }_{I}I-{\mu }_{T}T-{\mu }_{T_c}{T}_c+dS\end{array}\right)\in {\mathbb{R}}^6$$

At the disease free equilibrium (DFE), $\overline{x}={\left(0,0,0,0,0,N\right)}^T$, the Jacobian matrix of $\mathcal{F}$ and $\mathcal{V}$ is

$$D\mathcal{F}\left(\overline{x}\right)=\left(\begin{array}{cccccc}0& \beta & \alpha \beta & 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)$$

and

$$D\mathcal{V}\left(\overline{x}\right)=\left(\begin{array}{cccccc}\gamma +d& 0& 0& 0& 0& 0\\ -\gamma (1-p)& {\sigma }_I(1-p_I)+p_I{\delta }_I+d& 0& 0& 0& 0\\ -\gamma p& 0& {\sigma }_A(1-p_A)+p_A{\delta }_A+d& 0& 0& 0\\ 0& -{\sigma }_I(1-p_I)& 0& d+{\mu }_T& 0& 0\\ 0& 0& -{\sigma }_A(1-p_A)& 0& d+{\mu }_{T_c}& 0\\ 0& \beta -p_I{\delta }_I& \alpha \beta -p_A{\delta }_A& {\mu }_T& {\mu }_{T_c}& d\end{array}\right),$$

respectively.

Since in our model there are 5 infected compartments, we have

$$F=\left(\begin{array}{ccccc}0& \beta & \alpha \beta & 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right)$$

and

$$V=\left(\begin{array}{ccccc}\gamma +d& 0& 0& 0& 0\\ -\gamma (1-p)& {\sigma }_I(1-p_I)+p_I{\delta }_I+d& 0& 0& 0\\ -\gamma p& 0& {\sigma }_A(1-p_A)+p_A{\delta }_A+d& 0& 0\\ 0& -{\sigma }_I(1-p_I)& 0& d+{\mu }_T& 0\\ 0& 0& -{\sigma }_A(1-p_A)& 0& d+{\mu }_{T_c}\end{array}\right)$$

and $J_4=d>0$. Therefore, the conditions in [57] are satisfied, and

$$V^{-1}=\left(\begin{array}{ccccc}\frac{1}{\gamma +d}& 0& 0& 0& 0\\ \frac{\gamma (1-p)}{(\gamma +d)({\sigma }_I(1-p_I)+p_I{\delta }_I+d)}& \frac{1}{{\sigma }_I(1-p_I)+p_I{\delta }_I+d}& 0& 0& 0\\ \frac{\gamma p}{(\gamma +d)({\sigma }_A(1-p_A)+p_A{\delta }_A+d)}& 0& \frac{1}{{\sigma }_A(1-p_A)+p_A{\delta }_A+d}& 0& 0\\ \frac{\gamma (p-1)(p_I-1){\sigma }_I}{(\gamma +d)(d+{\mu }_T)({\sigma }_I(1-p_I)+p_I{\delta }_I+d)}& \frac{{\sigma }_I(1-p_I)}{(d+{\mu }_T)({\sigma }_I(1-p_I)+p_I{\delta }_I+d)}& 0& \frac{1}{d+{\mu }_T}& 0\\ \frac{\gamma p(1-p_A){\sigma }_A}{(\gamma +d)(d+{\mu }_{T_c})({\sigma }_A(1-p_A)+p_A{\delta }_A+d)}& 0& \frac{{\sigma }_A(1-p_A)}{(d+{\mu }_{T_c})({\sigma }_A(1-p_A)+p_A{\delta }_A+d)}& 0& \frac{1}{d+{\mu }_{T_c}}\end{array}\right)$$

and

$$F{V}^{-1}=\left(\begin{array}{ccccc}\frac{\beta \gamma (1-p)}{(\gamma +d)({\sigma }_I(1-p_I)+p_I{\delta }_I+d)}+\frac{\alpha \beta \gamma p}{(\gamma +d)({\sigma }_A(1-p_A)+p_A{\delta }_A+d)}& \frac{\beta }{{\sigma }_I(1-p_I)+p_I{\delta }_I+d}& \frac{\alpha \beta }{{\sigma }_A(1-p_A)+p_A{\delta }_A+d}& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right).$$

Because of our assumptions on the parameters and since $F{V}^{-1}$ upper triangular,

$${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)=\frac{\beta \gamma (1-p)}{(\gamma +d)({\sigma }_I(1-p_I)+p_I{\delta }_I+d)}+\frac{\alpha \beta \gamma p}{(\gamma +d)({\sigma }_A(1-p_A)+p_A{\delta }_A+d)}$$

where $\rho \left(A\right)$ is the spectral radius of matrix A.

6.2. Derivation of the Endemic Equilibrium

To find the endemic equilibrium, one can solve first

$$\begin{array}{cc}\hfill 0& =p\gamma E-(1-p_A){\sigma }_{A}A-p_A{\delta }_{A}A-dA\hfill \\ \hfill 0& =(1-p)\gamma E-(1-p_I){\sigma }_{I}I-p_I{\delta }_{I}I-dI\hfill \\ \hfill 0& =(1-p_A){\sigma }_{A}A-{\mu }_{T_c}{T}_c-d{T}_c\hfill \\ \hfill 0& =(1-p_I){\sigma }_{I}I-{\mu }_{T}T-dT\hfill \end{array}$$

(8)

to find $A_{\infty },I_{\infty },{T_c}_{\infty }\mathrm{and}{T}_{\infty }$ in terms of $E_{\infty }$. Using this solution in

$$\begin{array}{cc}\hfill 0& =\mathrm{\Lambda }N-\frac{\beta S}{N}\left(I+\alpha A\right)+p_A{\delta }_{A}A+p_I{\delta }_{I}I+{\mu }_{T}T+{\mu }_{T_c}{T}_c-dS\hfill \\ \hfill 0& =\frac{\beta S}{N}\left(I+\alpha A\right)-\gamma E-dE,\hfill \end{array}$$

(9)

one can obtain $S_{\infty }$ in terms of ${\mathcal{R}}_0$. Finally, one can derive $E_{\infty }$ in terms of ${\mathcal{R}}_0$ from

$$N=S_{\infty }+E_{\infty }+A_{\infty }+I_{\infty }+{T_c}_{\infty }+T_{\infty }.$$