Computational Approaches of Modelling Human Papillomavirus Transmission and Prevention Strategies: A Systematic Review

System dynamics models implement simulations as a set of time dependent differential equations that describe interactions between different subpopulations and the transition between compartments [52]. Within fixed time steps, simulations recursively solve the differential equations and capture the dynamics of epidemic transmission over time [53]. For instance, the system dynamics version of SIR model is defined as:

The differential equations become:

The initial condition is generally defined as:

Therefore, equations (5)-(7) have a solution for $\Delta t=1$, and as $\Delta t\to 0$, the solution approaches the solution of equations (1)-(3) [53]. System dynamics approach strongly emphasises cohort-based modelling; for more precise simulation results, the population is often stratified into small cohorts according to age, gender, ethnicity, sexual activity, etc. System dynamics models tend to be more flexible to include variables without available empirical data support, and the process of modelling is both iterative and participatory [52].

However, a point to be noted that system dynamics models mentioned in empirical publications [52, 48] only include models implemented with differential equations like the Kermack and McKendrick model introduced above, which are known as deterministic models. In contrast to deterministic approach, stochastic approach introduces randomness to modelling, using transitional probabilities instead of rates to model epidemic process. To a degree, deterministic models can be regarded as an approximation of stochastic models [54, 29]. Common stochastic models include discrete time Markov chain (DTMC) models, continuous time Markov chain (CTMC) models, and stochastic differential equation (SDE) models; these approaches varies in terms of the underlying assumptions of time and state variables. In a DTMC model, both time and state variables are discrete. In a CTMC model, state variables remain discrete, while time is continuous. SDE models are based on a transmission process, and both time and state variables are continuous. The principal difference between deterministic models and stochastic models is that even though the deterministic model converges to an endemic equilibrium, the corresponding stochastic model can converge to the disease-free state. The stochastic approaches also have unique properties including the probability of an outbreak, the quasi stationary probability distribution, the final size distribution of an epidemic and the expected duration of an epidemic, which can not be found in deterministic models [54].

Here, we assign stochastic models to system dynamics models as a different transition setting approach.

In epidemics, infectivity is not homogeneous; only a small number of infective individuals known as the super spreaders can transmit the infection to many other susceptible individuals [55]. To represent the heterogeneity, Keeling et al.[56, 57, 58] proposed that many infections can be considered to be transmitted through a limited network of contacts with individuals denoted by nodes and contacts between individuals denoted by links. A node in a network is not necessarily an individual; it can also be a cohort or a subpopulation. In unweighted networks, the transmissibility between two nodes is homogeneous; in weighted networks, the transmissibility between two nodes is denoted by the weight of the link. In undirected networks, the infection can be transmitted both ways along the links; in directed networks, the infection is transmitted following the direction of the links. Generally, it is assumed that a node’s transmissibility depends on its degree, the number of links it has [29].

Topology of the underlying network varies according to transmission routes and infection types [58]. The networks commonly investigated in the literature of computational epidemiology include lattices, small-world networks, scale-free networks and random networks, etc [59, 29]. Of special importance to study epidemic transmission dynamics is scale-free networks. In a network, the degree distribution is denoted as ${p_{k}}$. In a random network, ${p_{k}}$ is distributed exponentially with $p_{k}$ approaching to zero rapidly. While ${p_{k}}$ of a scale-free network presents a ‘fatter tail’, and the quantity $p_{k}$ relatively slowly decreases to zero as $k\to\infty$. This corresponds to a situation that there is an active core group and also super spreaders with high degrees in the context of epidemics [55]. Scale-free networks exhibit power-law degree distribution whereby the probability mass function (pmf) of degree $P(k)$ is distributed as:

where $A$ is a coefficient, and $\gamma$ is known as the scale-free exponent (also referred as the power-law exponent), with its value typically ranging from 2.0 to 3.0 for real-world scale-free networks [60, 55, 61, 62].

The network-growth based approach assigns each node to a neighbourhood that they can only infect their neighbours and be infected by the neighbours directly [59]. The study of Pastor et al. [60] exemplified that on a wide range of scale-free networks, the critical threshold of epidemic as $R_{0}=1$ in deterministic models can be absent, and the infection may benefit from the network topology regardless of the transmission rate. Moreover, both a low average degree and a high value of clustering of a network reduce the basic reproduction number $R_{0}$, comparing to the initial growth rate, which can influence the estimation of epidemic parameters in early stages of the epidemic [56, 57]. Therefore, different types of networks can result in varied epidemic behaviours; the underlying topology of a contact network can profoundly affect the infection transmission [56, 58, 57, 55, 48, 15].

Network-based models have advantages of representing heterogeneously distributed population and interaction patterns among individuals; numerical simulations can help investigate temporal dynamics of epidemic spread [15]. When the contagion is transmitted much faster than the evolution of the contact network, the network can remain static. On the other hand, the network can be adaptive or dynamic to be comparable with the epidemic process. There is increasing attention in amplifying complex network modelling to simulate spatio-temporal epidemic transmission dynamics, and network-based interventions are an active research area in computational epidemiology [48, 15, 29, 55].

In this section, only static networks are considered. Models with underlying static network topology are assigned to network-based models, while the ones built with adaptive or dynamic networks are recognised as agent-based models, which will be introduced next.

In agent-based models, each individual within the population is modelled as a distinctive agent with attributes and behaviours; interactions between agents provide opportunities to transmit contagion among the population [63]. Interdependent agents make interactions to generate “histories”, and collective properties and behaviours of the overall system will be manifested in simulations [52, 64]. Agent-based models offer the ability of implementing dynamic interactions between agents, rather than creating static networks, and they can also simulate sophisticated collective phenomena such as neighbourhood effects, markets, social capital, etc [48].

Agent-based models represent agents with various attributes and behaviours, a variety of structural and relational constrains and the entwines between local behaviours and global outcomes [65]. Characteristics and behaviours of agents can be implemented in programming languages, such as R, Java, C++ and Python, etc. and also in environments like NetLogo[66], Repast Suite [67], AnyLogic [68] etc., which provide visualisation of real-time dynamics. The contact patterns and mobility of agents are determining factors of contagion dissemination. Agents change their location through the simulation and exhibit varied behaviours spatio-temporally, which can shed light on why social distancing is required during a pandemic [63]. To better mimic the reality, geographic information systems (GIS) have been utilised to retrieve demographic and environment data to facilitate visualisation of epidemic outbreaks in geographical landscapes [47, 15].

Agent-based models have been shown to be a useful analysis tool in public health research, and are suited in particular to study epidemics dynamics and infectious disease transmission at multiple scales varying from individual communities to the global population[52]. Their flexibility and generality are their main strengths. Agent-based models hold considerable promise for multi-disciplinary fields, which is not possible for other modelling approaches to handle with similar ease, and the underlying formalism prevents researchers from making assumptions with limited knowledge of the situation [64].

Integrating large-scale census datasets can considerably improve the accuracy of epidemiological models; precise implementation of the underlying contact network also plays a crucial role in the simulations [15, 16, 69]. However, the implementation of agent-based models requires data availability, and the simulation need high-performance computational tools. Moreover, with the wide data collection of human behaviour and mobility, the underlying algorithms may be too complex to be formalised [15].

A key point to be noted here is how we make a deliberate distinction in this review between papers which undertake agent-based modelling, and papers which undertake network growth-based modelling, to capture HPV disease dynamics. Indeed, a contact network could be inferred from an agent-based model based on the interactions between individual agents, so it may seem that such an approach uses both agent-based and network-based modelling. However, the primary modelling approach in such cases would be agent-based modelling, and the contact network would be an observable feature of such a model. Such models are classified as ‘agent-based models’ in this review. By contrast, network growth-based models are assumed to be those where the relevant community is primarily captured by a network-growth model, with emphasis on interaction patterns and topological features, rather than individual attributes. In other words, agent-based models are assumed to be those which use high-fidelity data, are computationally intensive, and focus on capturing attributes and transactions of people at individual level (and may be able to indeed express some interaction patterns as networks), while network-growth models are assumed to be those models which do not primarily focus on individual attributes but capture interaction patterns in a topological sense. Therefore, network growth models would typically be less computationally intensive, and they achieve this by capturing interactions only topologically, not spatially. To stress that the modelled system is created by network growth, and not simply represented as a network, in such models, we call them ‘network growth-based models’ rather than simply ‘network-based models’.

HPV infection can be transmitted through sexual contacts; however, traditional approaches of reducing sexually transmitted infections are not typically applicable to HPV infection [1, 8]. Limiting number of sexual partners, practising abstinence and delaying coitarche can reduce the risk but not eliminate it [1]. HPV infection has high transmissibility and is prevalent in all sorts of sexual contacts [1, 8]. Similarly, circumcision alone cannot prevent HPV infection. Medical data shows high prevalence of HPV infection among circumcised men and their sexual partners, and no evidence has proved that circumcision can reduce the incidence. Condom use is similarly ineffective [8]: even when a condom is used, the exposed part of genital skin may get infected [1].

WHO issued a call to eliminate cervical cancer in 2018 and launched the global strategy to accelerate the elimination in 2020. The global strategy includes three complementary pillars: 1) 90% HPV vaccination coverage of eligible girls; 2) 70% screening coverage using a high-performance test; and 3) At least 90% of women with a positive screening test or a cervical lesion to be treated properly [13].

The recognised best way of preventing cervical cancer is to immunise girls before they are sexually active. Comparing to natural infection, HPV vaccination is more immunogenic and triggers higher polyclonal antibody response, and the avidity of the response will not increase after boosting [2]. HPV vaccines are administrated to prevent cervical premalignant lesions and cancers induced by high-risk HPV infections. The first HPV vaccine was licensed in 2006, and so far, six prophylactic vaccines have been approved. The details of currently available HPV vaccines are listed in Table 2.2. All HPV vaccines contain virus-like particles (VLPs) similar to HPV types 16 and 18, which are the two most prevalent genotypes [8, 2]. Empirical experiments suggest that all licensed HPV vaccines perform with high efficacy in populations that are not exposed to HPV infections before [2]. Following the onset of sexual activity, adolescents and young women are at the highest risk of acquiring HPV infection, and the risk decreases as age increases [8, 4, 3]. Therefore, HPV vaccination primarily targets girls aged 9-14 years before their sexual debut. After assuring the high coverage on the primary target population, the secondary target population can include females aged 15+ years, boys in general, and specifically men who have sex with men (MSMs). Till 2022, 125 countries have added HPV vaccination to their national immunisation program for girls, and 47 of them include boys [2].

The current HPV vaccination cannot provide full protection against all high-risk HPV types, and its impact on women who are vaccinated after the recommended age or not vaccinated is limited. It is argued that, as a primary prevention against cervical lesions and cancer, HPV vaccination cannot replace the role of screening later in life [13, 2].

To prevent cervical cancer, multiple screening methods are available to identify women who are at risk of, or have cervical precancerous lesions and early invasive cancer [13, 2]. The conventional approach is the Papanicolaou (Pap) smear test, also known as the cytology test. The cytology test was introduced in the 1940s to detect atypical cells in the cervical epithelium through preparation and interpretation of slides using microscopy. This initiated the transition of switching the emphasis to preventing invasive cancer before it happens [8, 13, 2]. In the last 15 years, newer approaches have been introduced. Visual inspection with acetic acid (VIA) identifies lesions that need treatment or further evaluation on the cervix without magnification. Nucleic acid amplification tests (NAAT) (DNA or mRNA) are performed on molecular level to detect early HPV infection [13, 2]. Newer techniques and tests have been developed for more precise screening results recently [13]. To accelerate the elimination of cervical cancer, WHO strongly recommends that women in general should receive HPV DNA tests as the primary screening test other than cytology or VIA in screening every 5-10 years starting at aged 30 years; for immunocompromised women, HPV DNA tests should be performed every 3-5 years starting at aged 25 years. If the primary screening test result is positive, responsive treatment can be administered based on the test result, which is called the “screen-and-treat approach”. Otherwise, females with positive results can receive a second test as a “triage” test with or without histologically confirmed diagnosis, which is called the “screen, triage and treat approach” [2, 13].

In this paper, we classify studies based on prevention approaches as follows:

• •

  Screening: screening is the only implemented prevention strategy in the model used by the study.

• •

  Vaccination: vaccination is the only implemented prevention strategy in the model used by the study.

• •

  Hybrid: both screening and vaccination are implemented in the model used by the study.

System dynamics models inherit the frameworks of mathematical models and incorporate them with computer simulations [47, 48]. System dynamics models are developed following the compartmentalisation concept that divides the population into exclusive compartments of different health statuses, as mathematical models do. The transitions that individuals move from one compartment to the next are captured as the epidemic dynamics [49, 29]. Generally, all individuals are susceptible to the infection in the beginning. When the infection seed is introduced to the population, the epidemic process starts. Susceptible individuals encounter infectious individuals, which enables the dissemination of the infection; susceptible individuals will be moved to the infectious compartment if they are infected during the encounter. The infectious individuals will be immune to the infection after they recover; these individuals will be moved from the infectious compartment to the recovered compartment. When there is no more individuals that can be infected, the epidemic will gradually die out. This is the basic susceptible(S)-infectious(I)-recovered(R) compartmental framework [49, 29]. If no long-term immunity can be acquired from the infection, infected individuals will be moved to the susceptible compartment after they recover, which is the SIS approach. In addition, if the encounter is supposed to happen at a rate, an E(exposed) compartment will be added between S and I compartments. The epidemic dynamics are similar to what is shown in Fig 2. Vaccinated individuals will be marked as an individual V(vaccinated) compartment in the process. More compartments can be added if needed [49].

Due to the relatively high volume of system dynamics models included in the literature, the models are further classified according to the following criteria:

• •

  Model Basis

• •

  Transition Setting

According to model basis, system dynamics models can be divided into cohort based models and individual based models. If the epidemic process is defined by a set of differential equations, for instance, as shown in Eq. [1]-[3], where the transitions are decided by constant rates, the model is deterministic. In other words, all individuals in the same compartment share the same transitional rates. For stochastic models, the differential equations are replaced by transitional probabilities with the introduction of randomness. Based on the criteria, the system dynamics models are divided into cohort deterministic models[72, 74], cohort stochastic models[77, 78, 79, 81, 84], and individual stochastic models[21, 97, 98, 99]. This is summarised in Table 3.1.1.

In cohort-based deterministic models [74, 72], the population is stratified into cohorts based on gender, age, etc. The progression and regression rates stay constant, and the epidemic process is modelled by differential equations. Berhe and Al-arydah [74] developed a heterosexual deterministic SIVS HPV model. The female population was divided into susceptible unvaccinated females, infected females and susceptible vaccinated females, and the male population was compartmentalised similarly. Susceptible unvaccinated females and males received impulsive vaccination periodically, and a certain proportion of susceptible vaccinated individuals might lose their immunity from the vaccination and become susceptible unvaccinated each year. During the period between two consecutive pulse shots, the model worked as a standard SIVS model. In a similar work, Jarynowski and Serafimovic [72] extended the classical SIR model. The male population was divided into Susceptible (Sm), Infectious (HPVm) and Recovered (Rm). The female population included “Long-term colonised (StageIIw)”, “Vaccinated (Vw)” and “Having cancer (Cancer)”, besides Sw, HPVw and Rw. All relationships shared the same infectivity, and partner change rates were stochastic according to sexual active levels and age. Progression and regression rates between compartments were implemented to be age-specific, and the epidemic process was again modelled by differential equations.

Cohort-based models [77, 78, 79, 81, 84] also utilised transitional probabilities from stochastic models. Burchell et al. [77] employed a heterosexual stochastic Monte Carlo simulation approach to implement sexual contact dynamics. Accordingly, each cohort shared the same intercourse frequency, and condom use frequency was assigned to each couple. HPV transmissibility was set as male-to-female and designated per coital act in 12 different levels. 27 HPV types were considered in the model. Van de Velde et al. [78] modelled four HPV genotypes: HPV-16, HPV-18, other high-oncogenic risk types and low-risk oncogenic types and incorporated dual between the four genotypes. The model implemented the natural history of HPV infection and SCC in a woman cohort during their lifetimes with six mutually exclusive compartments: susceptible, immune, infected, CIN1, CIN2/3, and SCC. In the model, individuals could develop lifelong immunity after infection or CIN clearance, or become susceptible again. The transitional probabilities were designed to vary depending on age and HPV genotype. It should be noted that here the model included a separate compartment, which derived from CIN2/3, representing individuals having hysterectomy. This cohort would not become susceptible or resistant to the infection. Coupé et al. [79] also followed health trajectories of women health over their lifetimes. They implemented 14 HPV infections in parallel and allowed concurrent infections. Progression probabilities were type-specific. In the model, females could develop CIN1 without HPV infection, and viral clearance and lesion regression were possible before CIN3. Progression probabilities were stochastic according to HPV genotype, age, and compartments. No immunity was implemented since this is an SIS model. Bogaards et al. [81] developed a dynamic model to mimic heterosexual partnership formation and also implemented 14 high-risk HPV infection. Only women were included in the SIRS model, whose structure is similar to the model developed by Van de Velde et al. [78] without the hysterectomy compartment. Though the model incorporated relationship formation, transmissibility was independent of sexual activity and infection stage. Landy et al. [84] ran HPV16/18 and non-16/18 high-risk HPV transmission independently among women aged 12-80 in an SIS model. HPV infections were classified into new HPV infection and persistent HPV infection with the same parameter for progression, and high-grade CIN would progress into asymptomatic cancer first then symptomatic cancer. The transitional probabilities were specific to age and type, and before asymptomatic cancer, infection and CIN could possibly clear spontaneously.

Similar to cohort-based stochastic models, individual-based stochastic models replace certain transitional rates with heterogeneous probabilities. Kim et al. [21] and Coupé et al. [97] model women from prior to their sexual debut till the very end of their life and include compartments of HPV infection and natural history of cervical cancer. Kim et al.[21] simulated four HPV strains: low-risk HPV, HPV-16, HPV-18 and other high-risk HPV; Coupé et al. [97] employed a Markov chain that runs independently for 14 high-risk HPV genotypes. Both models implemented transmission probability and progression rates as type dependent and regression rates as homogeneous. Similarly, Kim et al.[21] included different stages of cervical cancer: local cancer, regional cancer and distant cancer, which can not spontaneously clear or be cured in the model. Ekwunife and Lhachimi [98] remodelled a patient level simulation of natural history for cervical cancer disease from an empirical model-based study [100]. The simulation commenced when women were 9 years old and terminated when they are 99 years old. The model unified 12 oncogenic HPV genotypes as a single HPV infection and signified nine mutually exclusive heath states. Here, before cancer state, all states could lead to mortality in general, and cervical cancer can be cured or lead to death due to cancer or in general. Transitional probabilities included both certain rates and stochastic probabilities depending on age. Datta et al.[99] developed a heterosexual SIRS-V model with compartments including susceptible, infected, recovered and vaccinated. Natural immunity acquired from HPV infection is defined here as having a short-duration; only vaccination can provide longer-lasting protection. The population is divided into cohorts according to age, sex, sexual orientation and previous sexual experience. Partnership formation is implemented, and nine independent HPV genotype transmissions are modelled independently. Transmission probability and recovery rate differ between females and males. The asymmetric transmissibility between partners depend on individual characteristics and HPV types, and the recovery rate only depends on sex.

Network-based approach extends system dynamics to graphs to emphasise the relationships among sets of actors [48, 52]. Network analysis has been widely utilised in public health especially disease transmission in this century [52]. Representing interpersonal relationships with abstract network can help analyse the underlying topology and understand individual behaviour and adaption, which give insights into contagion transmission and prevention strategies [19]. Back in the early years of the HIV epidemic, epidemiologists incorporated the fundamental SIR model with the local social network information to visualise the disease transmission and investigate feasible interventions [52].

In network-based models, individuals or subpopulations are denoted by nodes, and the interplay between them are represented by links. Generally, the transmissibility of a node depends on the degree of the node, which is the number of links it has [29]. Two networks with the same size can display various structures, and different topologies can represent varied transmission routes and infection types [58]. The networks commonly investigated in the literature of computational epidemiology include lattices, small-world networks, scale-free networks and random networks [59, 29]. Empirical studies demonstrated that the underlying topology of a contact network can impact the infection transmission over the network [56, 57, 58, 55, 15]. For instance, the epidemic threshold $R_{0}$ approaches zero in a scale-free network of infinite size so that large outbreaks can emerge with fairly low transmissibility [48].

In the existing literature that we review, there is only one model [104] implemented which has an underlying static network; to better introduce the model, variations of this [101, 102, 103, 105] are included, as shown in Table 3.1.2. Díez-Domingo et al. [104] constructed a random network model to mimic sexual contacts. The original network was designed as a heterosexual network, then nodes representing men who have sex with men (MSM), who might also have female sexual partners, were later added into the model [103, 102, 105], In the sexual contact network, each node denoted an individual, and basic characteristics including age and gender were randomly assigned according to population histogram. The total number of lifetime sexual partners (LSP) representing sexual behaviours was distributed to each node proportionally considering assortativity, which is the attribute quantifying the tendency of nodes connecting with similar nodes in networks [112, 113, 114, 62, 115]. Most individuals had four or less LSPs. Formation of relationships depended on age and LSP for both partners, and the weight function pairing nodes $i$ and $j$ was defined as follows:

where $nLSP[i]$ denoted the expected LSP of node $i$. $Age[i]$ was the age of node $i$, and $1.8$ was the average age difference of couples in Spain back then. To form an assortative network, the lower the weight was, the more likely the relationship would form. Nodes with similar age and similar behaviours represented by LSP had more chances to pair with each other; couples where one owned few LSP and the other had way more LSP would be uncommon [104, 105]. During the network formation, women were prioritised to be paired in descending order of LSP in order to optimise the pairing process, and each MSM was paired with a woman to avoid isolated network components [101]. In the network, sexual contact frequency was age-dependent [104].

The primary model [104] simulated two HPV strains: high-risk HPV and low-risk HPV, and later the infection types were stratified to specific HPV genotypes [105]. Transmission probability depended on the age groups where the two individuals belong to, infection type and transmission route. Men were estimated to have higher probabilities to transmit contagion to women. Here, randomness was introduced to the model. Certain probabilities were set as epidemic parameters over the network. In order to constrain that a contagion was transmitted by sexual contacts, random numbers were generated to see if they were under the given thresholds decided by the epidemic parameters. If the random numbers met the requirements, the transmission was successful. When a susceptible individual was infected during sexual contacts, the infection duration (clearance time) would be randomly assigned from the expected intervals [102]. No natural immunity or reinfection was mentioned in the model.

Agent-based models use bottom-up implementation to simulate epidemic transmission dynamics by examining how individual’s attributes, their interactions and the environment result in collective emergent behaviour [52, 15]. In agent-based models, an agent represents a distinctive individual with traits and behaviours, and the agents interact with each other, which provides opportunities to transmit contagion among the population [63]. Similar to network-based models, agent-based models also utilise underlying network infrastructure to represent connections among the agents; however, to demonstrate the interaction dynamics, agent-based models construct adaptive or dynamic networks as opposed to static networks in network-based models [48]. In the dynamic system, the local contact patterns and mobility change consecutively as the agents’ behaviours vary throughout the simulation, which further impact the global contagion transmission [65, 63]. In agent-based models, there is no formalised definition of the global system behaviour like system dynamics models, which result in a form of decentralised modelling [52]. To understand the collective properties and behaviours of the global system, local interactions are captured during simulations [52, 64]. Agent-based models also offer the benefit of modifying agents’ attributes, behaviours, interactions and epidemic parameters to assess their impact on producing the population-level outcome and generate more complex phenomena such as social capital, neighbourhood effects, etc. [48, 63]. With the advantage of solving problems with a multi-disciplinary approach and providing realistic approximations, agent-based models are becoming increasingly popular in computational epidemiology [64, 63].

In the existing literature that we review, three agent-based models [106, 20, 23] are proposed, as shown in Table 3.1.2.

Olsen and Jepsen [106] constructed a closed “mini-society” of 25,000 individuals that constituted a heterosexual contact network in the simulation. The agents were aged 10-78 years, and the distribution mimicked the general population in Denmark back then. The gender ratio was approximately 50%. Here, the agents’ interactions are implemented by if/then logic depending on agents’ attributes instead of equations, and the sexual act frequency is assigned by a random-gamma distribution with a mean of 9.48 per month. It should be noted that this model allowed concurrent partners, and each relationship was assigned a duration from a normal distribution based on age in months, as follows:

where the older the agents were, the longer the relationship maintained. Four HPV genotypes, HPV-6/11/16/18 corresponding to quadrivalent HPV vaccination, were allocated to each age group according to the prevalence in the Danish population. The risk and duration of infection were implemented type-specific and stochastic. As low-risk HPV strains, HPV types 6 and 11 were implemented with a different progression history from the high-risk strains; agents with HPV-6/11 infection could only progress to genital warts instead of cervical lesions. Female agents with high-risk HPV infection could progress to CINs and further cervical cancer. No natural immunity is implemented in the model; therefore, after viral clearance, the risk of being reinfected by the same infection would not change. During the health trajectories, the progression and regression probabilities were defined by estimated parameters in years.

Van de Velde et al. [20] developed a stochastic individual-based HPV transmission model with sequential partnership formation and dissolution, which was later named as HPV-ADVISE [108, 109, 110, 111]. Individuals joined the simulated population prior to sex debut, and they were assigned characteristics including age, gender, sexual activity level, onset of sexual activity, current health status and current partnership status upon entering the simulation [20, 107, 108, 109, 110, 111]. Without the HPV infection, the model was referred as the demographic model. In the demographic model, the population stayed open and stable, in which individuals enter and die following the mortality rate specific to age and gender to mimic the national demography [20]. During the simulation, potential sexual behaviour, health-seeking behaviour (vaccination and screening), etc. are recorded for each agent [108]. The underlying sexual contact pattern was represented by a stochastic pair formation and separation process. Sequential monogamous stable and instantaneous relationships were implemented in the model; moreover, this model also considered about casual sexual relationships between men and female sexual workers (FSW). Partnership formation and termination were initiated by female agents; the rates of new partner acquisition and partnership termination was related to females’ age, sexual activity level and current partnership status [20]. Single females formed a new partnership at a rate, which could become stable with a probability. To form the partnership, age and sexual activity level on both sides would be considered, which implied assortativity in the sexual contact pattern; the type of partnership is based on the proportion of female agents’ having sex before marriage. In a stable relationship, the rate of separation, total number of sexual contacts and partnership duration solely depended on the female’s age and sexual activity level. For instantaneous partnerships, a small value is assigned as the total number of sexual acts with a fairly short duration [20, 108]. The model implemented transmission and natural progression history of 18 HPV types including both low-risk and high-risk strains, and each type was modelled individually without co-infections; therefore, an individual can be infected by multiple genotypes simultaneously [107]. In the natural history of HPV infection, infected females could clear the infection and become immune or susceptible, or remain infected or progress to more severe stages; females with CIN might regress to a less severe stage or even fully recover to immune or susceptible status [20]. The infection transmission was implemented per sex act, and the model allowed male-to-female and female-to-male transmission probabilities to be different. Since the transmission probability per sex act was gender and type-specific, and the frequency of sex acts was decided by age and level of sexual activity of the female; therefore, the transmission probability per relationship was age, gender and level of sexual activity-specific. Infection clearance rates were related to infection type, patient’s age and gender. Only females could develop life-long immunity from infection by stochastic processes [20]. To be noted, though 18 HPV types were represented in the model, quadrivalent vaccine was implemented against HPV-16, 18, 6 and 11. Three characteristics of the vaccine were modelled: 1) vaccine uptake probability; 2) vaccine efficacy per sex act; and 3) duration of protection sampled from a normal distribution. The baseline scenario was the vaccination coverage reached 100% with 95% efficacy, and the duration was either lifelong or averaged 20 years [20, 108].

Matthijsse et al. [23] employed an empirical STD model, STDSIM [116] to investigate HPV transmission dynamics of HPV-16 and 18. In STDSIM, each individual is assigned with both constant characteristics (e.g. sex, birth, etc.) and variable ones (e.g. heath status, number of sexual partners, etc.). The probability distributions of sex debut for men and women were defined separately, and for each individual in the simulated population, the age of first sexual contact was determined by Monte Carlo sampling. In the simulation, all events were arranged in the appropriate time order, which was called an event-driven approach. Characteristics of individuals or relationships were updated when an event happened. The occurrence of an event might generate new events, which would occur immediately or later in the simulation, or delay or cancel scheduled events. The event-driven approach in STDSIM was more efficient than a time-based approach in checking if events occurred at fixed time steps and cutting off unnecessary checks. STDSIM is the only model taking immigration and emigration into account. The demography module of the model included birth, death, immigration and emigration to represent the general population. The sexual behaviour module depicted how relationship were formed and terminated. Pairing took place on the basis of age. Individuals becoming available for a relationship and seeking for a partner were determined by probabilities, and their preferences of the relationship duration were defined by a function of age, personal interest and whether they have other sexual partners already, in which concurrent relationships were allowed. The frequency of sexual contacts within a relationship depended on the age of the man. STDSIM also modelled commercial sex, specifying the frequency of each man visiting commercial sex workers (CSW) [116, 23]. Matthijsse et al. [23] simulated HPV-16 and HPV-18 transmission in the network on the basis of STDSIM following the demography in the Netherlands. They assumed male-to-female and female-to-male transmissions sharing equal transmissibility per coital act, and the transmission parameters were calibrated to be age and type-specific. Comparing to empirical experiments on other STDs, they assumed condom uptake did not have protection against HPV, which agreed with prevention approaches mentioned in Section 2.2.1. Explicitly, they assumed women with hysterectomy were immune to HPV infections and modelled the fraction of hysterectomy following the general age distribution. Durations of infection and acquired immunity followed a Weibull distribution that was a continuous probability distribution. They employed a grid search to choose the optimal combinations of transmission probability and distributions of acquired immunity and infection durations. As a result, duration of acquired immunity with substantial variation played a crucial part in HPV epidemiology; without it, a model could not reproduce the age distribution of HPV prevalence.

WHO introduced the global strategy to accelerate elimination cervical cancer in 2020, in which 70% of eligible women should receive HPV screening at least twice in their lifetimes [13]. The conventional screening approach has been the cytology test, introduced in the 1940s, which has reduced the mortality from cervical cancer five-fold over the past 50 years; however, due to limited resources, cytology has not been promoted in low- and middle- income countries (LMICs) [8, 13]. Newer technologies have provided new options for screening. VIA allows examining by naked eyes or magnifying by colposcope or camera [13]. Since cervical cancer and all precursor lesions are attributed by high-risk HPV infection, NAAT can identify possible precursors on molecular level and prevent unnecessary treatments of plausible lesions, and it is more reproducible and can be applied under all settings [8, 13]. In empirical experiments, NAAT detected the same amount of cases as cytology but earlier [8]. As a result, WHO recommends that women in general should receive HPV DNA tests every 5-10 years starting at aged 30 years; for immunoscompromised women, HPV DNA tests should be performed every 3-5 years starting at aged 25 years [2].

In the literature, there are three publications [42, 95, 96] focusing on screening in LMICs, and all employed the same HPV transmission model [21]. In the three publications [42, 95, 96], types, frequencies, ages, visits and following treatments of screening were compared to evaluate different screening strategies; to analyse screening’s efficiency in reducing cervical cancer risks, immediate treatment were implemented after positive screening results came out, which are shown in Table 3.2.

Campos et al. [42] compared three types of screening tests: 1) the careHPV (provider-collected [cervical] and self-collected [vaginal sampling]); 2) VIA; and 3) conventional cytology. Here, the careHPV is an economically sustainable HPV DNA test introduced in places where resources are limited [117]. The model [42] assumed, women receiving screening during the first visit and returning for the results during the second visit, for careHPV; eligible women with positive screening results would mainly receive cryosurgery on the same day. In the VIA group, most eligible women who were screen-positive would receive immediate cryosurgery, but a part of them would either come back for a subsequent visit or drop out; those who were not eligible would be referred to a secondary facility for further diagnosis and treatment. Cytology required multiple visits. Women needed to be screened for the first visit and come back for results during the second one. Screen-positive women would be performed diagnostic colposcopy and biopsy during the third visit and visit a fourth time for treatment if necessary. Each screening test was performed between the ages of 25 and 50 years, at frequencies of once to three times in a lifetime, at 5 or 10 year intervals. 70% of women were randomly selected at each target age; women who missed one screening would not be selected next time. The results suggested that careHPV both cervical sampling and vaginal sampling were more effective and cost-effective in all settings comparing to the other two approaches. Within the age range from 30 to 45 years, screening three times could have the optimal and cost-effective performance, reducing cancer risk by 50%.

In the same year, Campos et al. [95] investigated the health and economic value of HPV screening algorithms in El Salvador, focusing on both screening and treatment. The simulated population in the model was women aged 30-65 years. First, for HPV-based screening, careHPV plan, women were supposed to receive screening every five years. The researchers compared two treatment plans following the positive screening results: 1) referral to colposcopy (Colposcopy Management [CM]) for women with CIN1 or higher stages at the designated hospital; 2) VIA triage followed by cryotherapy for eligible women in the clinic or colposcopy for ineligible women at the hospital (Screen and Treat [ST]). For Pap screening, women are supposed to uptake screening every two years followed by a referral to colposcopy for positive results (Pap). Here, they assumed 65% of the simulated women receiving administered screening and the left 35% never receiving screening. HPV-based screening demonstrated more efficiency in reducing cervical cancer risk by  60%, and ST was recognised as the most cost-effective strategy.

In 2017, Campos et al. [96] compared the 2- and 1- visit schedules of HPV DNA test to estimate the health and economic impact of reducing screening visits in LMICs. The model was calibrated to clinic data from India, Nicaragua and Uganda. In the simulation,each women were administrated to be screened three times in their lifetime starting at age 30 with five years interval, which reflected the results from their empirical study [42]. For 2-visit schedule, women would receive screening at the first visit and return for the results. the rate of loss to follow-up (LTFU) was varied from 10% to 70% with 10% increment, and each level was compared to 1-visit schedule. For 1-visit schedule, women would receive results on the same day of screening. Both strategies provided same-day cryotherapy at the result visit, and the performance of the two strategies were assumed comparable. The simulation results suggested the potential of developing 1-visit screening strategy, which was also called screen-and-treat.

HPV vaccination can induce avid antibody or polyclonal antibody response against designated HPV genotypes, which human bodies cannot develop from infection. The current prophylactic vaccines are indicated to be administrated before individuals are exposed to HPV, as vaccination can reach the optimal efficacy before the onset of sexual activity [2].

In the literature, the following publications utilised computational models to evaluate different HPV vaccination strategies. According to the position paper published by WHO [2], the primary target population of receiving HPV vaccination should be girls aged 9-14 years before their sexual debut, as prevention of cervical cancer is best achieved by immunising adolescent girls. When the vaccination coverage is high enough to achieve herd immunity, which is the threshold representing the proportion of the population that need to be immunised to the infection to sufficiently contain the further transmission [118]. Once the coverage of HPV vaccination in girls reaches 80%, the herd immunity can also benefit boys. Females over 15 years old, boys, men or MSM are listed as the secondary target populations [8, 2]. To deliver a more comprehensive investigation on HPV vaccination in computational models, we categorised the literature according to target vaccine uptake populations, as shown in Table 3.2.

Al-arydah and Smith [70], and Al-arydah and Malik [71] modelled vaccinating females only. They considered vaccination of teenage girls, as suggested by [2], but also modelled catch-up vaccination of females up to the age 26. In particular, Al-arydah and Smith [70] developed an aged-dependent two-sex deterministic SIR model to study the effect of the bi-valent HPV vaccination program targeting HPV-16 and 18. Their model employed a time and age dependent system of partial differential equations (PDEs) for simulating HPV transmission. The vaccine in their model is given to females in two stages: the childhood vaccine for ages under 13 and the adult vaccine for ages 13-26. It is assumed that the possibility of receiving the childhood vaccine does not depend on age as it is free of charge, while the chance of paying for the adult vaccine increases with age since adult women will have higher income and more awareness of the need for the vaccine as they grow older and gain more sexual experience. Their results suggested that the effects of age dependency in vaccine uptake are complex, but that vaccinating a single age cohort in one gender can result in eventual control of the disease across all age groups. Meanwhile Al-arydah and Malik [71] introduced an age-structured SIS model based on a nonlinear partial differential equation system to study the effect of a bi-valent HPV vaccination program targeting HPV-6 and -11. HPV-6 and -11 are common low-risk HPV genotypes and can cause genital warts, and recovery cannot provide permanent immunity. Their results suggest that early and catch-up female vaccine program eliminates the infection in both male and female populations over 30 years. They introduced the optimal control to an age-dependent model based on ordinary differential equations and solved it numerically to obtain the most cost-effective method for introducing the catch-up vaccine into the population. Though the framework is essentially mathematical, the numerical solutions used computational methods. Their results showed that, to ensure the eradication of HPV types 6 and 11, not only vaccines with high efficacy are needed, but it is also needful to increase the proportion of early females vaccinated and use optimal catch-up vaccines. The results imply that the catch-up vaccine can be administered optimally by using the maximum possible rate at the earliest possible age right away after females become sexually active.

Díez-Domingo et al. [104] devised a heterosexual contact network to simulate the HPV vaccination campaign against genital warts in Australia in 2007. In the project, two scenarios were simulated:

• •

  Scenario 1: 83% routine vaccination coverage of girls under 14 years old plus 73% catch-up vaccination coverage of women aged 14-26 years;

• •

  Scenario 2: 73% routine vaccination coverage of girls under 14 years old plus 52% catch-up vaccination coverage of women aged 14-26 years.

The vaccination efficacy was assumed 96.5%. In the simulation results, prevalence of genital warts in both women and men dropped after the introduction of vaccination, which revealed that non-vaccinated men could also benefit from the herd immunity by vaccinating females only.

The following three publications [83, 99, 105] considered adding boys to the immunisation program. Datta et al. [99] constructed an individual-based stochastic model simulating nine types of HPV transmission in the UK. They tested three current HPV vaccines, the bivalent, quadrivalent and nonavalent vaccines, under the “historical vaccination” background that only girls were vaccinated during 2008-2016. In the results, vaccinating girls only was the most cost-effective strategy; adding boys to an already successful girls-only immunisation plan was not cost-effective since men could be protected by the herd immunity as well. Villanueva et al. [105] adjusted an existing heterosexual network [104] with MSM implemented and simulated three vaccination scenarios:

• •

  1. 75% of 12- to 13-year-old girls were vaccinated;

• •

  2. 75% of 12- to 13-year-old boys were vaccinated;

• •

  3. 75% of 12- to 13-year-old both girls and boys (gender neutral) were vaccinated;

Similar to the results above, vaccinating boys was not cost-effective; however, MSM population could not be protected by the herd immunity when only girls were vaccinated. Vaccinating boys only could have partial control of infections in the MSM population but would put females at substantial risk. Therefore, the most effective campaign would be vaccinating both genders. Bogaards et al. [83] investigated the most effective sex-specific immunisation strategy and compared scenarios of vaccinating girls only, boys only and girls and boys. The objective was to vaccinate the sex with the highest prevalence of infection in order to achieve the largest reduction in population prevalence; however, the results were counterintuitive. Vaccinating the sex with the highest prevaccine prevalence yielded different results in SIS and SIR systems; vaccinating individuals at the highest infection risk was most effective in the SIR system, while in the SIS system, individuals with longer infectious periods needed prioritising. Moreover, vaccinating the sex with the higher prevaccine prevalence could not reach the expected performance if MSM were included in the model, and to reduce most prevalence in MSM, only boys could be vaccinated. However, in the conclusion, increasing vaccination coverage among girls was acknowledged as more effective since most existing HPV vaccination programs have achieved sufficient coverage with female-only vaccination.

In the following two models [74, 103], both females and males were vaccinated. Berhe and Al-arydah [74] proposed a novel SIVS HPV model with impulsive vaccination at regular intervals to examine the impulsive vaccine control strategy. Fractions of HPV naive girls and boys were administrated pulse vaccination before the onset of sexual activity. As a system dynamics model, the model’s global stability of disease-free periodic solution agreed with the condition based on basic reproduction number (for not having an epidemic) $R_{0}\leq 1$ as mentioned in Section 2.1. The results suggested that a larger pulse vaccination rate administrated over a longer period might help in eradicating the disease. In addition, the authors state that the model was limited without implementation of acquired immunity. Muñoz-Quiles et al. [103] analysed the prevalence decline in HPV infection under scenarios of 75% and 90% vaccination coverage in an existing sexual contact network [104]. The results suggested that once the vaccination coverage in the sexual contact network reached 70%-75%, improving the coverage showed limited variation in prevalence reduction. Meanwhile, the results again confirmed the importance of gender-neutral vaccinating program in providing health equity for MSM.

Three further works which model HPV vaccination for both males and females are Madhu and Al-arydah [73], Al-arydah (2021) [75], and Al-arydah (2023) [76]. Madhu and Al-arydah [73] introduced a two-sex age-structured cohort deterministic SIS model to describe the dynamics of HPV disease with childhood and catch-up vaccines, and discussed the impact of the age-difference of partners on the HPV prevalence and optimal vaccines. They proved that an optimal vaccine solution exists for this scenario and it is unique. They showed that in this scenario, 77% coverage of childhood vaccination is essential and controls HPV in all ages in a 20-year period. They demonstrated that childhood vaccination is more efficient than adult vaccination in controlling HPV diseases in youth and adults. Meanwhile Al-arydah (2021) [75] formulated a two-gender SIS model to search for optimal childhood and catch-up vaccination strategies against HPV-6 and -11 over a 20-year period. The populations of the two genders were modelled to grow logistically. They showed that to minimise the effective reproduction number under resource constraints in terms of HPV vaccinations, the priority for vaccines should be given to girls. They also determined minimum threshold levels of catch-up vaccine coverage separately for both genders, which is necessary to avoid an HPV epidemic (to keep the effective reproduction number less than unity). Similarly, Al-arydah (2023) [76] formulated an SIS cohort deterministic model that incorporates population behaviour and the impact of education on HPV dynamics. The goal was, specifically, assessing health education success in controlling HPV through increasing vaccine uptake and reducing risky behaviour in individuals. The model again assumed logistically growing populations. Al-arydah estimated the time series solution of the model for uneducated, lowly-educated and highly-educated populations. They showed that the higher the education level, the more the vaccine uptake, the less the risky behaviour, and the shorter the time needed to control HPV. They demonstrated that late education for adults is more effective than early education for children as late education is provided in tandem with vaccination programs, and that could increase vaccine uptake for parents and their children. Essentially they argued that increasing health education is a viable low-cost strategy to support vaccines that have low efficacy.

To achieve the maximum prevention of cervical cancer, HPV vaccination and screening should be practised together [2]. As the primary prevention intervention, HPV vaccination can induce persistent and avid antibody response that cannot develop from virus clearance. However, the current HPV vaccines cannot cover all oncogenic HPV genotypes and prevent some cancer-related infections, and vaccination can only provide partial protection to females. Therefore, the prevention of cervical cancer needs the secondary assistance of screening [8, 2]. With the diagnosis of screening, responsive treatment can be administrated to defer further progress of cervical lesions. Even if the global vaccination coverage were achievable today, the role of screening as the secondary prevention of cervical cancer cannot be replaced [8, 13, 2].

Considering the large quantity of HPV transmission models with hybrid prevention strategies implemented in the literature, we further categorised the related publications according to the following criteria:

• •

  Screening oriented

• •

  Vaccination oriented

• •

  Combination of vaccination and screening

The details are listed in Table 3.2.

Five publications [87, 80, 97, 84, 92] investigated cost-effectiveness of screening with different vaccination settings in the model.

Goldhaber-Fiebert et al. [87] and Burger et al. [92] investigated the cost-effectiveness of cervical cancer screening for both vaccinated and unvaccinated females. Goldhaber-Fiebert et al. [87] used an empirical individual-based stochastic model [21] to examine screening strategies for unvaccinated older women and girls eligible for HPV vaccination. In the implementation, prevention strategies included screening along, vaccination alone and vaccination accompanied by screening, where girls were vaccinated before age 12 years. The researchers examined cytology and HPV DNA testing and varied the combination of the initial age (18, 21, or 25 years) and interval (every 1, 2, 3, or 5 years). The testing strategies included:

• •

  Cytology with HPV test triage: HPV DNA testing after cytology for suspicious cytologic results;

• •

  HPV test with cytology triage: cytology after HPV DNA testing for positive HPV DNA diagnosis;

• •

  Combined HPV DNA testing and cytology.

The simulated individuals were allowed to switch testing strategies once at age 25, 30 or 35 years. The results suggested the most cost-effective screening strategy for both vaccinated and unvaccinated female is to use HPV DNA testing as triage for equivocal results in young women and primary screening approach in older women basing on age. Burger et al. [92] employed the same primary model [21] to find a cost-effective screening program in Norway. They compared the current cytology program with HPV DNA testing in for women who were (not) vaccinated before adolescence. Similar to Goldhaber-Fiebert et al.’s conclusion, Burger et al. [92] suggested to switch primary screening from cytology to HPV DNA testing in older women for cost-effectiveness.

The following three papers [80, 84] investigated cost-effectiveness of screening after vaccination. Coupé et al. [80] applied a cohort-based stochastic model [79] to evaluate the optimal approach to screen women after vaccination. They compared cytology and HPV DNA test as primary screening instrument. As results, screening five times with HPV DNA test or seven times with cytology between 30 and 60 years were comparable in cost-effectiveness. Later, Coupé et al. [97] compared screening scenarios after different types of vaccination. The first one is a bivalent HPV vaccine with partial cross-protection against HPV-31, 33, 45 and 58; the second one is a broad spectrum vaccine against HPV-5 to 13. In the simulation, screening once per lifetime was more cost-effective than receiving a broad spectrum vaccination. For females receiving the bivalent vaccine, they were suggested to take screening four times per lifetime, while for the ones taking broad spectrum vaccination, one screen was more cost-effective. Landy et al. [84] proposed a microsimulation model to determine the appropriate screening routine for females after vaccination. The simulation results suggested that females vaccinated against HPV-16 and 18 needed three screening tests during their lifetime, females vaccinated against HPV-16, 18, 31, 33, 45, 52 and 58 required twice, and the unvaccinated females should receive seven lifetime screens.

The majority of papers with hybrid prevention strategy implemented are vaccination-oriented, as shown in Table 3.2. To better deliver the implementation in the models, the related publications are divided into the following subcategories based on vaccinated population:

• •

  Girls

• •

  Females (girls and young women)

• •

  Females and males

In the models [85, 88, 106, 98, 82, 108, 109], adolescent girls were prioritised to be vaccinated as guided by WHO [2]; on the premise of ensuring the vaccination coverage in girls, young women [90, 79, 78] and boys or males [91, 107, 20, 110, 111] were included.

Goldie et al. [85] and Kim et al. [88] used an individual-based simulation model [21] to explore the health and economic impact of HPV vaccination in Brazil and Vietnam. In the model [21], girls were vaccinated, and females received cervical cancer screening two or three time per lifetime since they were 30. HPV DNA test and cytology were compared as the primary screening approach. Both Goldie et al. [85] and Kim et al. [88] suggested the combination of vaccination and screening could achieve the optimal reduction in cervical cancer risk; however, for LMICs, without assistance, they might not afford the vaccines, and screening only ended up the most cost-effective approach. Olsen and Jepsen [106] constructed an agent-based heterosexual contact model to analyse the incremental impact of introducing HPV vaccination compared to screening alone in Denmark. In the model, the screening program covered 70% of women aged 23-59 years in Denmark, who received screening every three years. Quadrivalent HPV vaccination was allocated to 70% of 12-year-old girls without catch-up program or booster dose in the analysis. The results suggested that though the introduction of quadrivalent vaccination would incur extra cost, the risk of cervical cancer and cost of treatment would be reduced. Ekwunife and Lhachimi [98] remodelled a microsimulation framework to estimate the cost-effectiveness of introducing HPV vaccination to Nigeria. They conducted scenarios in the simulation including:

• •

  Current opportunistic cervical cancer screening (CS): VIA once per lifetime with 8.7% coverage;

• •

  National cervical cancer screening (NS): VIA once per lifetime at age 30 for each woman;

• •

  National HPV vaccination with opportunistic screening (CS + NV): two doses of bivalent HPV vaccines for preadolescent girls at age 12 in addition to opportunistic screening;

• •

  National screening and national HPV vaccination (NS + NV): two doses of bivalent HPV vaccines for preadolescent girls at age 12 and VIA for women at age 30 once per lifetime

The CS + NV strategy was the only cost-effective and robust option when the cost of vaccines were relatively low enough, which was a problem for LMICs.

Bogaards et al. [82] studied the heterosexual transmission of HPV infection [81] to disentangle the direct and indirect effects of HPV transmission. They implemented the current cytologic screening program in the Netherlands in the model and simulated vaccination uptake among girls. The results demonstrated that the reduction in indirectly averted cancer cases peaked at 50-70 % vaccination coverage, which means vaccinating preadolescent girls could reduce transmission of vaccine-preventable HPV infection and lower infection rates in unvaccinated women. Van de Velde et al. [108] employed the HPV-ADVISE model [20] to compare population-level impact of the bivalent, quadrivalent and nonavalent HPV vaccines. In the model [108], both vaccination and screening were implemented as potential heath-seeking behaviour. There were five levels of screening behaviours representing the average interval between two routine cytology test when there was no equivocal result, ranging from every 1.25 years (Level 0) to never (Level 4); responsive treatment approaches were implemented to corresponding test results. Female’s screening rates depended on screening behaviour level, previous screening test results and age. Sensitive analysis was conducted on the three different types of vaccines. In the results, the three types of vaccines displayed different advantages. The bivalent vaccine was more effective at preventing high level lesions and cervical cancer in the long run, while the quadrivalent vaccine substantially reduced anogenital warts soon after introduction. They also suggested switching to the nonavalent vaccine could further decrease the risk of precancerous lesions and invasive cervical cancer.

Malagón et al. [109] calibrated the HPV-ADVISE model [20] to Canadian epidemiological data. The simulated population was stratified into four levels according to sexual activity (low = L0, high = L3), where the most sexual active groups represented a minority of the population but accounted for almost half of the prevalence. All girls were vaccinated at age 12 years with 100% efficacy against HPV-16 and 18, and the protection will last 20 years on average in the base case scenario. In the simulation, one-way sensitivity analysis was performed on vaccine uptake by sexual activity level. The results implied that the impact of uniform vaccination coverage on population effectiveness might be limited since the coverage was not high enough to produce herd immunity in high risk group; HPV vaccine uptake should be distributed heterogeneously according to risk of HPV-16 and 18 infection.

The following three models [90, 79, 78] considered the inclusion of young women in the immunisation program.

Kim and Goldie [90] and Coupé et al. [79] investigated the cost-effectiveness of HPV vaccination in females. Kim and Goldie [90] adopted a dynamic heterosexual model [21] to study the health and economic implications of HPV vaccination in the U.S. In the model, all prevention strategies contained screening, which was implemented following the local reported pattern. The researcher allocated catch-up vaccination from age 13 to 18, 21, or 26 in addition to routine vaccination at age 12 for girls and young women. To be noted, in the model, they assumed that unvaccinated females were also less likely to receive screening. The results suggested that the cost-effectiveness could be optimised by improving vaccination coverage in preadolescent girls and age threshold of catch-up vaccination to 18 or 21 years of age. The health impact of HPV vaccination also depended on the immunity duration, and screening policies needed revising as the introduction of vaccination. Coupé et al. [79] developed a simulation model accommodating both screening and vaccination. In the model, females aged 30-60 were invited to cytology every five years, and 85% 12-year-old girls were assumed to receive three doses of HPV vaccines. In addition, they separately investigated the impact of administrating catch-up vaccination at 30 years of age. The following scenarios were conducted:

• •

  No prevention

• •

  Screening only

• •

  Vaccination only

• •

  Combination of screening and vaccination

The simulation results confirmed the health impact of vaccination in reducing risk of cervical cancer and supported that vaccination should be introduced along with screening to achieve better effectiveness.

Van de Velde et al. [78] developed a stochastic model of natural history of HPV infection to predict the impact of quadrivalent HPV vaccination and measure parameter uncertainty. The probabilities of being screened was set as age-specific, and the sensitivity of test was lesion-specific. Sensitive analysis was performed on the quadrivalent HPV vaccination regarding to factors as follows:

• •

  Vaccine efficacy: 85%, 90%, 95%, 100%

• •

  Vaccine duration: life, 30 years, 30 years plus booster at age 30 years, 30 years plus booster every 10 years

• •

  Vaccinated age: 12, 15, 20, 30 years

As results, the effectiveness of vaccination might be influenced by vaccinated age, vaccine efficacy and immunity duration; vaccinating girls at 12 years of age with high efficacy and lifetime protection could reduce considerate risk of HPV infection.

Five other publications [91, 107, 20, 110, 111] accounted males in the target vaccinated population.

First, Kim and Goldie [91] and Brisson et al. [107] investigated the incremental value of including boys in the immunisation program. Kim and Goldie utilised the empirical model [90] introduced above, in which both girls and boys aged 12 years were vaccinated as inputs in the context of cervical cancer screening. Vaccinating both girls and boys were compared to vaccinating girls only, and the cost of including boys outweighed the benefits it would bring. As a result, the inclusion of boys needed iteratively revising. Brisson et al. [107] simulated vaccinating 12-year-old boys besides girls in the HPV-ADVISE model [20], where the vaccination could provide up to 20-year protection with 99% efficacy and 70% coverage. There was an incremental reduction in HPV incidence over 70 years; however, as the vaccination coverage increased in girls, the incremental benefit of vaccinating boys became limited.

Van de Velde et al. [20] proposed the HPV-ADVISE model to explore factors influencing predicted effectiveness of HPV vaccination, where screening was implemented as potential behaviour in the biological process. Three different vaccination strategies were compared: 1) girls only, 2) girls and boys, and 3) girls with a 1-year catch-up program up to age 18 year, in which characteristics of vaccination were varied, including efficacy, duration, coverage and uptake age. The results implied that the prediction of HPV vaccine effectiveness was sensitive to vaccine efficacy, duration and coverage. Both adding boys and catch-up vaccination could help with HPV infection elimination.

Laprise et al. [110] and Drolet et al. [111] both employed the HPV-ADVISE model [20]. Laprise et al. [110] performed non-inferiority trial on multiple vaccination strategies:

• •

  No vaccination vs. two-dose girls-only;

• •

  Two-dose girls-only vs. three-dose girls-only;

• •

  Two-dose girls-only vs. two-dose girls and boys;

• •

  Three-dose girls and boys vs. two-dose girls and boys or three-dose girls-only.

In the simulation, the efficacy duration of two doses was varied between 10 years and lifelong, and that of three doses was varied between 20 years and lifelong. Under most scenarios, vaccinating boys was not cost-effective unless the cost for boys was substantially lower than that for girl; if the protection duration of two doses reached 10 years, two-dose routing was likely to be cost-effective. Drolet et al. [111] estimated incremental effectiveness of administrating catch-up quadrivalent vaccine to females and males in Australia. Simulations were performed from no vaccination to routine vaccination, routine plus catch-up 18 years and routine plus 26 years with vaccinating girls only or girls and boys. The simulation results exemplified that catch-up vaccination produced rapid declines in prevalence and strong herd immunity effects along with high coverage of routine vaccination. The added benefits brought by catch-up vaccination could last for 20-70 years.

The following publications [86, 89, 93, 97, 72] explored the combination of vaccination and screening as prevention interventions in HPV transmission models.

Goldhaber-Fiebert et al. [86], Diaz et al. [89] and Campos et al. [93] used the same individual-based stochastic model developed by Kim et al. [21]. Goldhaber-Fiebert et al. [86] included screening and vaccination implementation in the primary model [21]. Simulated prevention interventions contained vaccination only, screening only annually or biennially, and combination of screening and vaccination. Practising screening or vaccination individually could reduce lifetime risk of cancer; vaccination combined with screening every five years achieved the highest reduction in the results.

Diaz et al. [89] and Campos et al. [93] examined health and economic impact of cervical cancer screening and bivalent HPV vaccination in India and Eastern Africa. Preadolescent girls under age 12 were vaccinated, women over age 20 got screened, and the combination of vaccination and screening were simulated in the model. The screening types included cytology, VIA and HPV DNA test, containing 1-3 visits for 1-3 times per lifetime at age range 35-45. Both publications [89, 93] agreed that vaccination combined with HPV DNA testing three times during lifetime was the most cost-effective prevention strategy; however, as LMICs, the coverage of vaccination being high enough to provide long-term population-level protection was only achievable if the cost of HPV vaccine was economic enough.

Jarynowski and Serafimovic [72] constructed an epidemic model and analysed impact of potential epidemiological control strategies in Poland in the future 25 years. They mimicked the actual screening situations in Poland as follows:

• •

  Until 2005: opportunistic screening without regular routine;

• •

  2006-2015: regular screening was introduced;

• •

  2016-2039: every 10 years, 1/3 women would do screening every three years.

As for vaccination, girls who were 14 years old were administered mandatory vaccines, and among girls aged 20 to 24, a subpopulation of 5%, opted to have catch-up vaccination. The results indicated that vaccination along with screening would be effective in controlling cervical cancer risk in Poland.