# Contact patterns and HPV-genotype interactions yield heterogeneous HPV-vaccine impacts depending on sexual behaviors: An individual-based model

Mélanie Bonneault, Chiara Poletto, Maxime Flauder, Didier Guillemot, Elisabeth Delarocque-Astagneau, Anne C.M. Thiébaut, Lulla Opatowski

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mélanie Bonneault, Chiara Poletto, Maxime Flauder, Didier Guillemot, Elisabeth Delarocque-Astagneau, et al.. Contact patterns and HPV-genotype interactions yield heterogeneous HPV-vaccine impacts depending on sexual behaviors: An individual-based model. Epidemics, 2022, 39, pp.100584. 10.1016/j.epidem.2022.100584. hal-03775055

HAL Id: hal-03775055

https://hal.science/hal-03775055v1

Submitted on 27 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Contact patterns and HPV-genotype interactions yield heterogeneous HPV-vaccine impacts depending on sexual behaviors: An individual-based model

Mélanie Bonneault<sup>a,b,c,*</sup>, Chiara Poletto<sup>d</sup>, Maxime Flauder<sup>a,b</sup>, Didier Guillemot<sup>a,b</sup>, Elisabeth Delarocque-Astagneau<sup>b</sup>, Anne C.M. Thiébaut<sup>c,1</sup>, Lulla Opatowski<sup>a,b,1</sup>

## ARTICLE INFO

### Keywords:

- Human papillomavirus

- Coinfections

- Individual-based model

- Sexual behavior

- Vaccination

## ABSTRACT

Human papillomaviruses are common sexually transmitted infections, caused by a large diversity of genotypes. In the context of vaccination against a subgroup of genotypes, better understanding the role of genotype interactions and human sexual behavior on genotype dynamics is essential. Herein, we present an individual-based model that integrates realistic heterosexual partnership behaviors and simulates interactions between vaccine and non-vaccine genotypes. Genotype interactions were considered, assuming a previous vaccine-genotype infection shortened (competition) or extended (synergy) the duration of a secondary non-vaccine-genotype infection. Sexual behavior determined papillomavirus acquisition and transmission: only 19.5% of active individuals at most 1 partner during the year, but > 80% of those with ≥ 2 partners, were infected before vaccine introduction. The pre-vaccination situation was consistent with all genotype interaction scenarios. These genotype interactions, despite being undetectable during the pre-vaccination era, markedly impacted genotype prevalence after vaccination started, with a significant increase/decrease of non-vaccine genotypes prevalence for respectively competitive/synergistic interactions. These prevalence changes were more pronounced in individuals with ≤ 3 partners per year (up to 30% of prevalence modification assuming 65% vaccine coverage) but barely visible for individuals with > 3 partners per year (at most 0.30%). Results suggest the presence of genotype interaction, which is consistent with the pre-vaccine situation, may impact the dynamics of non-vaccine genotypes, particularly in less active individuals.

## 1. Introduction

Human papillomavirus (HPV) genital tract infections are among the most common sexually transmitted infections, especially in younger people (de Sanjosé et al., 2007; Markowitz et al., 2013). HPVs are characterized by a high diversity of genotypes, a few of which have been associated with cancers (cervical, vaginal, vulval, penile and anal) (Walboomers et al., 1999). Prevention has recently relied on vaccination with bivalent, quadrivalent and, more recently, nonavalent (Serrano et al., 2012) vaccines that target the subtypes carrying the highest

carcinogenic risks.

Following the introduction of anti-HPV vaccines, HPV-infection prevalences with targeted genotypes have declined in several countries (Drolet et al., 2019; Tota et al., 2020). However, quantifying the vaccines' impacts on global HPV prevalence remains difficult. In particular, an increased prevalence of genotypes not included in the vaccine, known as genotype replacement, is a potential risk. Therefore, when vaccinating against a subgroup of genotypes, prediction of vaccine impact requires better understanding of whether genotypes interact and, if so, how.

<footnote><sup>a</sup> Epidemiology and Modelling of Antibiotic Evasion Unit, Institut Pasteur, 75475 Paris Cedex 15, France</footnote>

<footnote><sup>b</sup> Université Paris-Saclay, UVSQ, Inserm U1018, CESP, Anti-infective Evasion and Pharmacoepidemiology Team, 78180 Montigny-le-Bretonneux Cedex, France</footnote>

<footnote><sup>c</sup> Université Paris-Saclay, Inserm U1018, CESP, High Dimensional Biostatistics Team, 94807 Villejuif Cedex, France</footnote>

<footnote><sup>d</sup> Inserm & Sorbonne Université, iPLESP, 75646 Paris, France</footnote>

<footnote>Abbreviations: HPV, Human papillomavirus; IBM, individual-based model; V, vaccine; NV, non-vaccine.</footnote>

<footnote>* Correspondence to: Institut Pasteur, Epidemiology and Modelling of Antibiotic Evasion Unit, 25, rue du Docteur-Roux, 75475 Paris Cedex 15, France.
E-mail address: melanie.bonneault@gmail.com (M. Bonneault).</footnote>

<footnote n="1">¹ These authors contributed equally to the research.</footnote>
In fact, simultaneous, within-host co-infection by two genotypes can affect each virus genotype's load, cell-infection ability and/or infection duration (McLaughlin-Drubin and Meyers, 2004; Murall et al., 2014; Xi et al., 2009). On the one hand, some evidence indicates that HPV genotypes may compete (Biryukov and Meyers, 2018; McLaughlin-Drubin and Meyers, 2004; Xi et al., 2009); on the other, multiple infections were commonly observed before vaccine introduction (Chaturvedi et al., 2011; Mejlhede et al., 2010; Mendez et al., 2005; Spinillo et al., 2009). The mechanisms leading to the co-existence of widely diverse HPV genotypes in populations have not yet been definitively elucidated. Some researchers argue co-existence is consistent with genotypes acting independently of each other (Chaturvedi et al., 2011; Vaccarella et al., 2013), while others suggest that the observed clustering of distinct HPV genotypes is due to between-genotype interactions (Mejlhede et al., 2010; Merikukka et al., 2011; Spinillo et al., 2009).

Previous HPV-transmission models formalizing interactions between HPV genotypes have not yet been able to reproduce simultaneously the complex aspects of HPV biology and individual behaviors. So far, models accounting for genotype interactions have been based on homogenous-mixing assumptions (Elbasha and Galvani, 2005; Pons-Salort et al., 2013; Poolman et al., 2008), although partnership may vary greatly across individuals and ages. These heterogeneities, in addition to affecting infection risk and HPV prevalence (Shiboski and Padian, 1996), could also markedly impact genotype co-infection, co-circulation and interaction dynamics at the population level. Indeed, theoretical study results suggested that the contact network affects the population ecology of pathogen strains (Eames and Keeling, 2006; Pinotti et al., 2019). Gray et al. consistently observed HPV-genotype-replacement differences between individuals with high-risk and low-risk behaviors (Gray et al., 2019). Therefore, an accurate accounting of sexual behavior might be essential to correctly interpret HPV observations and provide more accurate projections of the epidemiological and ecological consequences of vaccination.

Optimally, only individual-based models (IBMs) enable reproduction of the heterogeneity of individual behaviors and simulation of their effects at population levels (Auchincloss and Diez Roux, 2008). A few IBMs have been developed for HPV (Johnson et al., 2018; Matthijsse et al., 2015; Olsen and Jepsen, 2010; Van de Velde et al., 2012, 2010) but none has considered HPV-genotype interactions. Herein, we investigate the role of heterogeneous contact behaviors and genotype interactions in the dynamics of HPV infection before and after the introduction of vaccination. How are infections and co-infections distributed across sexual activity groups? How does vaccination impact the prevalence of vaccine (V) and non-vaccine (NV) genotypes among sexual activity groups? To answer these questions, we developed and used an HPV-transmission IBM based on a realistic heterosexual network of individuals and accounting for the most prevalent HPV V and NV genotypes.

# 2. Materials and methods

## 2.1. Global overview

We developed a stochastic HPV-transmission IBM accounting for contact heterogeneity according to sex, sexual activity and age. We considered the 14 most prevalent genotypes, including two V (HPV-16 and -18), and 12 NV genotypes (HPV-31, -33, -35, -39, -45, -51, -52, -56, -58, -59, -66 and -68). Here each genotype was modeled individually, enabling either single transmission or simultaneous transmission of multiple genotypes in case of a multi-infected partner. Using this framework, we explored a range of competitive and synergistic genotype interactions. The model assumed natural acquired immunity and included vaccination implementation, assuming the vaccine provided full immunity against V-genotype-HPV infections.

## 2.2. Population and individual characteristics

Our model was restricted to the ages at which HPV transmission and prevalence peak so that individuals of both sexes enter the population at 15 years and leave it at 30 years. Individuals who exit the population are directly replaced by new 15-year-old individuals, keeping the population stable over time. The population is divided equally by sex and age per year. Each individual is explicitly modeled for each week and characterized by his/her age, sex, sexual activity class (which determines the number of partners the individual has during the year), partnership status (three categories: in a relationship, available for a new partnership, or neither—concurrent partnership is not allowed) and infection status for each genotype (four categories: susceptible, infected, naturally immunized or vaccinated in case of a V genotype).

## 2.3. Modeling partnership

We modeled heterosexual partnership behaviors as a stochastic process according to the steps summarized in Fig. 1 and detailed in Supplementary material S1.

Model parameters describing the partnership were estimated to reproduce the reported contact patterns from the *Contexte de la Sexualité en France* (CSF) survey (Bajos and Bozon, 2008). For each combination of parameters to be calibrated, the distribution of individuals by sex, age categories (18–19, 20–24 and 25–29 years) and categories of their number of sexual partners since age 15 years obtained by simulations was compared to survey data using the least square distance method. The calibrated sexual behavior parameters are summarized in Supplementary material S2.2. Combinations of parameter values were simulated. After reaching prevalence equilibrium, we extracted the cumulative number of sexual partners from age 15 years to that at the time of data extraction for each individual. Because women's and men's behaviors could not be reproduced simultaneously, we first selected the combinations consistent with women's data, then chose the best combination among them with respect to men's data.

## 2.4. Modeling HPV-genotype transmission and infection

The transmission of any of the 14 genotypes can occur when an infected individual is in partnership with an individual not infected with that specific type (Fig. 1). Transmission-probability parameters are defined by genotype group for V and NV genotypes, respectively $β_V$ and $β_{NV}$, assuming two contacts per week. If the virus is transmitted, acquisition occurs. If an individual is infected—regardless of the genotype—his/her infection duration is sampled from an exponential distribution of mean 52 weeks, in accordance with the literature (Trottier et al., 2008). After infection clearance, the individual is assumed to have acquired natural immunity, conferring total protection against the same genotype for a duration drawn from an exponential distribution whose mean is identical for all genotypes. The distribution of immunity durations is identical for all genotypes. At the end of the immunity period, the individual becomes again fully susceptible to that genotype. The transmission-probability parameters, $β_V$ and $β_{NV}$, and mean immunity duration were calibrated by minimizing the sum of squared differences between simulated and reported prevalence before vaccine introduction by age category (Markowitz et al., 2013), first assuming genotypes to be independent with respect to transmission and infection (neutral interaction scenario, details in Supplementary material Table S4, in S2.3).

With those assumptions, the 14 viruses spread simultaneously among partners in the population. Fig. 2 illustrates transmission over a simulated network.

## 2.5. Between-genotype interaction

We focused only on V–NV-type interactions. When an individual is co-infected with two HPV genotypes, we assumed that their co-infection
<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Schematic representation of partnership and HPV infectious processes. The partnership process is guided by the rules of sexual activity defining active and inactive periods. Simultaneous partnership is not allowed in the model.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Partnership and HPV-transmission network throughout 1 year with a simulated population composed of 1058 individuals and 565 connections. Node colors correspond to age categories. Node diameters correspond to the number of persons infected by the individual over the 1-year study period after reaching prevalence equilibrium. If two individuals were in relationship during the year, a linking line is drawn between them, colored in red if ≥ 1 transmission(s) occurred between them, otherwise black.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Calibration of partnerships and HPV-transmission parameters. (a) Cumulative distribution of the number of sexual partners since the beginning of their sexual life for women 18–30 years old. (b) NV- and V-genotype prevalences with 95% confidence interval (CI) reported in Markowitz et al. (2013) (circles) and results of 10 simulations, each obtained as the average value over 10 years after reaching prevalence equilibrium (connecting lines).</figcaption></figure>

a power law (Liljeros et al., 2001; Schneeberger et al., 2004). Furthermore, the characteristics of the sexual contact network appeared stable over time as shown in Figs. S4–S6 in paragraph S4.2 of Supplementary material.

Under the neutral interaction scenario, the model accurately reproduced realistic distributions of V- and NV-genotype prevalences by age category, with the characteristic bell-shaped curves and prevalence peaks between 20 and 24 years (Fig. 3b). The transmission-probability parameter calibrated for V genotypes was higher than that of NV genotypes ($\beta_V = 0.16$ and $\beta_{NV} = 0.125$, respectively, with 12 weeks of calibrated immunity duration). This difference was not significant according to a profile likelihood approach (paragraph S4.3). In simulations, the average percentage of ever-infected individuals increased sharply with age, reaching 71.4% at 29 years. This simulation result is

consistent with previous reports that 80% of people are infected at least once in their lifetime (Santé publique France, 2019) considering that a limited number of new infections can be acquired after the age of 30 years and that only 14 high-risk HPVs among the 37 detectable sexually transmitted genotypes have been modeled here.

## 3.2. Infections and co-infections did not spread homogeneously throughout the network under the neutral interaction scenario

Individuals infected with a single or multiple genotype(s) were distributed among the population as a function of the number of partners during the past year (Fig. 4a). After reaching prevalence equilibrium, 33.3% of the active population were infected (27.7% of the whole population) and 20.3% were co-infected (61.5% of the infected). This

<figure><img src="image_8.png" /><figcaption>**Fig. 4.** NV- and V-genotype infections, co-infections and transmissions according to the annual number of partners. (a) Mean percentages of individuals with single and multiple infections among the active population over 20 simulations according to their number of partners during the past year at prevalence equilibrium $t = 640$ weeks. (b) Mean percentages of individuals infected with a single NV or V genotype, with the two V types but no NV type, with $t > 1$ distinct NV types but no V type, and with V and NV types among the infected population over 20 simulations according to their annual number of partners at prevalence equilibrium time. (c) Percentages of virus transmissions according to sexual activity group of the transmitter during the year. For each transmission group (x-axis), the sexual activity group distribution of the acquirers is represented on the y-axis.</figcaption></figure>
result, although at the upper limit of the range, is consistent with values described in observational studies, which range from 20% to 70% co-infection among infected individuals (Chaturvedi et al., 2011; Mejlhede et al., 2010; Mendez et al., 2005; Spinillo et al., 2009). Among the infected population (Figs. 4b), 45.2% of infected individuals had 0–1 partner during the past year (15.5% prevalence for this group), 43.2% had 2–3 partners (81.7%) and 11.6% had > 3 partners (99.5%). The co-infection percentage increased with annual sexual activity, with only 8.2% of active individuals with 0–1 partner being co-infected, as opposed to > 60% of individuals with ≥ 2 partners. Co-infections with multiple V genotypes but no NV genotype were extremely rare (0.35% of infected individuals). Indeed, 33.1% of the infected population were V–NV-genotype co-infected.

As shown in Fig. 4c, HPVs were mostly transmitted by individuals with ≥ 2 partners (41.2% by individuals with 2–3 partners and 47.5% by individuals with > 3 partners in the year); 95.5% of individuals with > 3 partners and 48.3% of individuals with 2–3 partners transmitted at least once during the year. Moreover, transmission mostly occurred within the sexual activity group (comparison with partnership matrix in Supplementary material S4.4).

### 3.3. Alternative genotype-interaction scenarios did not modify pre-vaccine infection and co-infection patterns

Before vaccination was introduced, all simulations were fitted to HPV-genotype prevalences within the range of assessed competitive, neutral and synergistic interactions by calibrating the transmission probabilities $β_{NV}$ with higher values in competition and smaller values in synergy (Fig. S9). Only small variations in prevalence were observed across interaction scenarios. The strength of interaction did not markedly affect co-infection patterns either (Fig. S15a-c). We further verified that this lack of influence did not simply reflect differences in transmission probability $β_{NV}$, resulting from calibrating the model for varying interaction strengths (Supplementary material S7.2).

### 3.4. V and NV prevalences were impacted differently according to sexual activity group post-vaccine introduction

After introducing vaccination, hypotheses regarding genotype-interaction impacted NV-prevalence trends (Figs. 5 and S10 for a uniform 65% vaccine coverage among women; results for 25% coverage are given in Supplementary material S4.7). In all genotype-interaction scenarios, V-prevalence range decreased (difference relative to pre-vaccine era) by 83.7–85.9% for individuals with 0–1 and 2–3 partners, and only 51.3–53.2% for individuals with > 3 partners (Fig. 5a). Conversely, NV-genotype prevalences changed markedly after vaccine introduction for all interaction strengths compared to the neutral scenario. Competitive or synergistic interaction scenarios, respectively, led to higher or lower NV prevalence compared to the neutral scenario (Figs. S10, and S15a and d). As expected, the stronger the interaction strength (the more it deviates from 1), the greater the NV-prevalence percentage difference (for additional information, see the Supplementary material on sensitivity analysis S7.3). The NV-prevalence effect was stronger in individuals with 0–1 partner (median 25.9% for $α = 0.67$ and –32.0% for $α = 1.5$) than in individuals with 2–3 partners (median 14.2% for $α = 0.67$ and –23.5% for $α = 1.5$) (Fig. 5b). For individuals with > 3 partners during the year, prevalence did not change appreciably (median 0.08% for $α = 0.67$ and –0.25% for $α = 1.5$).

Sensitivity analyses showed that targeting high-risk sexual activity groups instead of vaccinating all women < 15 years old led to greater reduction of V-type prevalence for individuals with > 3 partners but did not modify NV-prevalence-variation patterns across sexual activity groups (Supplementary material S4.8). Furthermore, among the other hypotheses tested for genotype interactions (Supplementary material S8), it appeared that only interactions by HPV-genotype groups could be detected within the considered range of interaction strengths. Moreover, when notable NV-prevalence differences were observed, compared to the neutral interaction scenario, they appeared only among the less active individuals.

## 4. Discussion and conclusion

We found that heterogeneous contact behaviors and natural

<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Percentage differences between post- and pre-vaccine prevalences of (a) V and (b) NV genotypes under 65% vaccine coverage by sexual activity group and interaction scenarios. Vaccine was introduced after reaching prevalence equilibrium. Average prevalences were assessed at three distinct times each before and after vaccine introduction. Boxplots show the variability over 20 simulations for each interaction scenario (bold horizontal line inside the box is the median, lower and upper box limits are the 1st and 3rd quartiles, antennae correspond to 1.5 times the interquartile range and circles mark the values outside this range).</figcaption></figure>
immunity acquired after infection were essential to reproduce the bell-shaped curve of HPV-infection prevalence by age category. Calibrated interaction scenarios revealed no notable differences in terms of single infection and co-infection patterns before vaccine introduction. However, thereafter, interactions mattered: competitive interactions led to significantly increased NV prevalence, while the opposite was seen for synergistic interactions. Of interest, our results suggested that the interaction effect on prevalence was stronger for less active individuals.

While it is expected that the introduction of the nonavalent vaccine will again modify the prevalence equilibrium, solid comprehension of the potential impacts of the initial bi- and quadrivalent vaccines is required to improve understanding of changing HPV ecology. Our IBM allowed us specifically to examine the impact of vaccination scenarios of V-, NV- or V-NV genotype interactions. Two former studies based on compartmental models analyzed the impact of genotype interactions on acquisition or clearance, without explicitly modeling heterogeneous sexual behaviors as IBMs do (Elbasha and Galvani, 2005; Pons-Salort et al., 2013). The results of both showed that, when considering between-genotype interactions, NV-genotype prevalences could be modified by vaccine introduction.

Based on our simulations for the prevaccine era, 71.4% of individuals were infected with at least one HPV genotype by 29 years of age. This simulation finding agrees with previous reports of 80% infected at least once during their lifetime (Santé publique France, 2019), acknowledging that a limited number of new infections can be acquired after 30 years of age and that only 14 high-risk HPV among 37 detectable sexually transmitted genotypes were modeled here. Furthermore, co-infections were frequent in our results, representing about two-thirds of HPV-positive women for all interaction strengths tested. That outcome, although at the upper limit of the interval, is consistent with 20-70% co-infections among all infections described in observational studies (Chaturvedi et al., 2011; Mejlhede et al., 2010; Mendez et al., 2005; Spinillo et al., 2009). Pertinently, herein, all co-infections were assumed to be detected and counted, unlike epidemiological studies in which, depending on the techniques used, detection of multiple infections may be more-or-less sensitive (Qu et al., 1997).

Concerning transmission probabilities under the neutral interaction scenario, our estimated $\beta_{NV}$- and $\beta_V$-transmission probabilities were 0.125 and 0.16 per contact, respectively. Considering two sexual acts per week on average, our results are of the same order of magnitude as values from other IBMs, range 0.048-0.95 per sex act (Johnson et al., 2018; Matthijsse et al., 2015; Olsen and Jepsen, 2010; Van de Velde et al., 2012, 2010). Nevertheless, our estimates strongly depend on infection duration and the definition of immunity. Notably, assuming a probability < 1 to acquire type-specific, natural, lifelong immunity yielded higher transmission probabilities, range 0.5-0.95 (Johnson et al., 2018; Van de Velde et al., 2012, 2010). In addition, either with acquired immunity during a limited period, as assumed herein, or without immunity, modeling showed lower transmission probabilities, range 0.048-0.3 (Matthijsse et al., 2015; Olsen and Jepsen, 2010).

Pertinently, the pre-vaccine situation could be explained by totally distinct interaction hypotheses. This is an important result since many epidemiological studies have claimed that pre-vaccination patterns are consistent with the absence of genotypic interactions only, although potential for biases cannot be ruled out (Malagón et al., 2016). Our finding is indeed consistent with other transmission models investigating biological mechanisms to explain pre-vaccine epidemiological data (Durham et al., 2012; Murall et al., 2014). Moreover, vaccine introduction more strongly impacted NV prevalence in less active individuals, whereas any prevalence modification would be expected to be seen in high-risk groups (Gray et al., 2019; Tota et al., 2013). Considering transmission throughout the network, the acquisition risk for any individual depends on two factors: the virus-circulation level in the population (defined as the transmission rate, and durations of infection and immunity), and the number of partners. As those numbers remain unchanged on average, the vaccination impact depends on

virus-transmission potential before vaccination. That potential was substantially higher for the > 3-partner group than the low- or intermediate-activity groups. Interestingly, sensitivity analysis on the impact of the proportions of individuals with a medium or high number of contacts confirmed the robustness of our results that prevalence variations of NV genotypes were stronger in low activity than highly connected individuals (Supplementary material S5).

Our results obtained with this methodology should also be interpreted in light of the following limitations. First, aggregated data were used, impeding a precise description of individual trajectories throughout their sexual lifetimes, such as switching between sexual activity groups, and durations of partnership and between two partnerships. As previously stressed in other studies, our model was unable to perfectly reproduce the distributions of the numbers of partners for both sexes. Indeed, CSF data included more men with high numbers of partners than women (Bajos and Bozon, 2008). Differential sex-dependent reporting bias was observed previously (Fenton et al., 2001; Mitchell et al., 2019) and may, in part, be at the origin of the discrepancy between model and data. To overcome that difficulty, we focused on women's data to fit our model and estimate parameters, while trying to conform to men's data in terms of cumulative partner numbers. Reassuringly, the resulting distribution of total numbers of partners was typical of a partnership network (Liljeros et al., 2001; Schneeberger et al., 2004).

Second, because no data are available on HPV-infection prevalence by age and V genotype for the pre-vaccine era in France, we used distributions reported for the US, whose current HPV epidemiology is similar (Bruni et al., 2010; de Sanjosé et al., 2007). NV-prevalence simulations appeared to be underestimated for 15-19-year-olds, compared with US data (Markowitz et al., 2013). Unfortunately, partnership and infection data are scarce for this age category, making calibration difficult. We think that underestimation for 15-19-year-olds had little impact on our results for all ages combined. The average prevalence of 27.7% for the whole 15-30-year-old population obtained in our simulations before vaccination is comparable to those reported in epidemiological studies in France, range 25.1-28.5% for 15-79-year-old women with normal vaginal smear cytology (Dalstein et al., 2003; Monsonego et al., 2005; Riethmuller et al., 1999).

Third, to keep our model relatively simple, only two transmission-probability parameters were estimated for all genotypes, despite studies having reported that HPV-16 and -18 prevalences may differ from one another and from those of other HR genotypes (Markowitz et al., 2013). Hence, our estimated transmission-probability parameter $\beta_V$ can be considered an average probability for HPV-16 and -18, and that of $\beta_{NV}$ an average probability for the 12 high-risk NV types. Others in previous studies have chosen to simulate only two groups of HPV V and NV equivalent to the circulation of two genotypes only. This choice impedes addressing the questions of multiple infections and interactions properly. In addition, it does not allow to define specific immunity for each genotype, resulting in poor calibration to the data (results not shown). Similarly, because epidemiological study results have not yet demonstrated clear between-genotype differences for infection duration (Ramanakumar et al., 2016; Trottier et al., 2008), a unique parameter was defined for all genotypes.

Fourth, while much is still unknown about immunity following HPV infection (Gravitt and Winer, 2017), we chose to define immunity as total protection against reinfection to the same genotype during a fixed duration, identical for all genotypes, in accordance with findings suggesting some degree of protection against reinfection (Ho et al., 2002; Safaeian et al., 2010). Previous HPV modeling considered either complete lifelong immunity, partial or waning immunity, acquired immunity during a defined period, or protection increasing with the number of past infections (Franceschi and Baussano, 2014). Matthijsse et al. showed that assuming natural immunity in their model was necessary to reproduce the age-specific HPV-infection patterns and that assuming full genotype-specific immunity yielded the best calibration results
(Matthijsse et al., 2015). We made the same observation herein, after testing various immunity hypotheses (data not shown). Moreover, allowing for greater heterogeneity in the distribution of immunity duration as previously suggested (Beachler et al., 2016), we confirmed our main result that NV variations are mostly apparent in low-activity individuals (Supplementary material S6).

Finally, while interaction mechanisms between HPV genotypes also remain unknown, we assumed that genotype interaction would affect the infection duration of a second acquired virus. Some authors suggested that simultaneous within-host co-infection with two genotypes could impact each virus-genotype load and thus the successful infection of cells by each (Biryukov and Meyers, 2018; McLaughlin-Drubin and Meyers, 2004; Xi et al., 2009). We explored alternative interaction hypotheses in sensitivity analyses and obtained similar results for NV genotypes, when genotype interaction impacted acquisition risk instead of infection duration.

To conclude, our results confirmed that HPV infections reached not only the most sexually active but also less active individuals and unexpectedly revealed that the impact of vaccine introduction on genotype ecology could be more detectable in less active individuals. These analyzes suggest that focusing on more sexually active individuals to address questions related to vaccination impact might potentially lower the ability to detect HPV prevalence variations, and highlight that public health information and prevention efforts should target all sexually active individuals, not only the population at-risk for sexually transmitted infections. Moreover, our results showed that genotype-interaction-associated hypotheses are plausible as they were all consistent with the reported pre-vaccine prevalences. In the post-vaccine era, better understanding those interactions and sexual contacts is key to anticipating the long-term impact of anti-HPV vaccines with respect to the prevention of cervical and other cancers.

## Funding sources

MB was funded by the INCa [Grant DOC 2017-123] and MGEN, and her work was supported by internal resources of Institut Pasteur, the French National Institute of Health and Medical Research (Inserm) and the University of Versailles Saint-Quentin-en-Yvelines (UVSQ).

## CRediT authorship contribution statement

MB, CP, MF, EDA, DG, ACMT and LO conceptualized the project. AT and LO supervised the project. MB, ACMT and LO designed the model. MB and MF developed the model. MB, CP, ACMT and LO performed the network and statistical analyses and validated them. All authors approved the latest version of this article.

## Declaration of Competing Interest

The authors declare no conflicts of interest.

## Data availability

No data was used for the research described in the article.

## Acknowledgements

We thank Margarita Pons-Salort for helpful discussions and ideas provided early during the project. We are also grateful to Janet Jacobson for editorial assistance.

## Appendix A. Supporting information

Supplementary data associated with this article can be found in the online version at doi:10.1016/j.epidem.2022.100584.

## References

Auchincloss, A.H., Diez Roux, A.V., 2008. A new tool for epidemiology: the usefulness of dynamic-agent models in understanding place effects on health. Am. J. Epidemiol. 168, 1–8. https://doi.org/10.1093/aje/kwn118.

Bajos, N., Bozon, M., 2008. Sexualité, genre et santé: les apports de l'enquête "Contexte de la sexualité en France". Enquête sur la sexualité en France. La Découverte.

Beachler, D.C., Jenkins, G., Safaeian, M., Kreimer, A.R., Wentzensen, N., 2016. Natural acquired immunity against subsequent genital human papillomavirus infection: a systematic review and meta-analysis. J. Infect. Dis. 213, 1444–1454. https://doi.org/10.1093/infdis/jiv753.

Biryukov, J., Meyers, C., 2018. Superinfection exclusion between two high-risk human papillomavirus types during a coinfection. J. Virol. 92 https://doi.org/10.1128/JVI.01993-17 (e01993-17).

Bruni, L., Diaz, M., Castellsagué, X., Ferrer, E., Bosch, F.X., de Sanjosé, S., 2010. Cervical human papillomavirus prevalence in 5 continents: meta-analysis of 1 million women with normal cytological findings. J. Infect. Dis. 202, 1789–1799. https://doi.org/10.1086/657321.

Chaturvedi, A.K., Katki, H.A., Hildesheim, A., Rodríguez, A.C., Quint, W., Schiffman, M., Van Doorn, L.-J., Porras, C., Wacholder, S., Gonzalez, P., Sherman, M.E., Herrero, R., Group, C.V.T., 2011. Human papillomavirus infection with multiple types: pattern of coinfection and risk of cervical disease. J. Infect. Dis. 203, 910–920. https://doi.org/10.1093/infdis/jiq139.

Dalstein, V., Riethmuller, D., Prétet, J.-L., Le Bail Carval, K., Sautière, J.-L., Carbillet, J.-P., Kantelip, B., Schaal, J.-P., Mougin, C., 2003. Persistence and load of high-risk HPV are predictors for development of high-grade cervical lesions: a longitudinal French cohort study. Int. J. Cancer 106, 396–403. https://doi.org/10.1002/ijc.11222.

de Sanjosé, S., Diaz, M., Castellsagué, X., Clifford, G., Bruni, L., Muñoz, N., Bosch, F.X., 2007. Worldwide prevalence and genotype distribution of cervical human papillomavirus DNA in women with normal cytology: a meta-analysis. Lancet Infect. Dis. 7, 453–459. https://doi.org/10.1016/S1473-3099(07)70158-5.

Drolet, M., Bénard, E., Pérez, N., Brisson, M., HPV Vaccination Impact Study Group, 2019. Population-level impact and herd effects following the introduction of human papillomavirus vaccination programmes: updated systematic review and meta-analysis. Lancet Lond. Engl. 394, 497–509. https://doi.org/10.1016/S0140-6736(19)30298-3.

Durham, D.P., Poolman, E.M., Ibuka, Y., Townsend, J.P., Galvani, A.P., 2012. Reevaluation of epidemiological data demonstrates that it is consistent with cross-immunity among human papillomavirus types. J. Infect. Dis. 206, 1291–1298. https://doi.org/10.1093/infdis/jis494.

Eames, K.T.D., Keeling, M.J., 2006. Coexistence and specialization of pathogen strains on contact networks. Am. Nat. 168, 230–241. https://doi.org/10.1086/505760.

ECDC, 2020. Guidance on HPV vaccination in EU countries: focus on boys, people living with HIV and 9-valent HPV vaccine introduction [WWW Document]. Eur. Cent. Dis. Prev. Control. (https://www.ecdc.europa.eu/en/publications-data/guidance-hpv-vaccination-eu-focus-boys-people-living-hiv-9vHPV-vaccine), (Accessed 7 February 2020).

Elbasha, E.H., Galvani, A.P., 2005. Vaccination against multiple HPV types. Math. Biosci. 197, 88–117. https://doi.org/10.1016/j.mbs.2005.05.004.

Fenton, K.A., Johnson, A.M., McManus, S., Erens, B., 2001. Measuring sexual behaviour: methodological challenges in survey research. Sex. Transm. Infect. 77, 84–92. https://doi.org/10.1136/sti.77.2.84.

Franceschi, S., Baussano, I., 2014. Naturally acquired immunity against human papillomavirus (HPV): why it matters in the HPV vaccine era. J. Infect. Dis. 210, 507–509. https://doi.org/10.1093/infdis/jiu143.

Gravitt, P.E., Winer, R.L., 2017. Natural history of HPV infection across the lifespan: role of viral latency. Viruses 9, 267. https://doi.org/10.3390/v9100267.

Gray, P., Luostarinen, T., Vänskä, S., Eriksson, T., Lagheden, C., Man, I., Palmroth, J., Pimenoff, V.N., Söderlund-Strand, A., Dillner, J., Lehtinen, M., 2019. Occurrence of human papillomavirus (HPV) type replacement by sexual risk-taking behaviour group: post-hoc analysis of a community randomized clinical trial up to 9 years after vaccination (IV). Int. J. Cancer 145, 785–796. https://doi.org/10.1002/ijc.32189.

Ho, G.Y.F., Studentsov, Y., Hall, C.B., Bierman, R., Beardsley, L., Lempa, M., Burk, R.D., 2002. Risk factors for subsequent cervicovaginal human papillomavirus (HPV) infection and the protective role of antibodies to HPV-16 virus-like particles. J. Infect. Dis. 186, 737–742. https://doi.org/10.1086/342972.

Johnson, H.C., Lafferty, E.I., Eggo, R.M., Louie, K., Soldan, K., Waller, J., Edmunds, W.J., 2018. Effect of HPV vaccination and cervical cancer screening in England by ethnicity: a modelling study. Lancet Public Health 3, e44–e51. https://doi.org/10.1016/S2468-2667(17)30238-4.

Liljeros, F., Edling, C.R., Amaral, L.A., Stanley, H.E., Aberg, Y., 2001. The web of human sexual contacts. Nature 411, 907–908. https://doi.org/10.1038/35082140.

Malagón, T., Lemieux-Mellouki, P., Laprise, J.-F., Brisson, M., 2016. Bias due to correlation between times-at-risk for infection in epidemiologic studies measuring biological interactions between sexually transmitted infections: a case study using human papillomavirus type interactions. Am. J. Epidemiol. 184, 873–883. https://doi.org/10.1093/aje/kkw152.

Markowitz, L.E., Hariri, S., Lin, C., Dunne, E.F., Steinau, M., McQuillan, G., Unger, E.R., 2013. Reduction in human papillomavirus (HPV) prevalence among young women following HPV vaccine introduction in the United States, National Health and Nutrition Examination Surveys, 2003–2010. J. Infect. Dis. 208, 385–393. https://doi.org/10.1093/infdis/jit192.

Mathijsse, S.M., van Rosmalen, J., Hontelez, J.A.C., Bakker, R., de Kok, I.M.C.M., van Ballegooijen, M., de Vlas, S.J., 2015. The role of acquired immunity in the spread of
human papillomavirus (HPV): explorations with a microsimulation model. PLoS One
10, e0116618. https://doi.org/10.1371/journal.pone.0116618.

McLaughlin-Drubin, M.E., Meyers, C., 2004. Evidence for the coexistence of two genital HPV types within the same host cell in vitro. Virology 321, 173-180. https://doi.org/10.1016/j.viroi.2003.12.019.

Mejlhede, N., Pedersen, B.V., Frisch, M., Fomsgaard, A., 2010. Multiple human papilloma virus types in cervical infections: competition or synergy? APMIS Acta Pathol. Microbiol. Immunol. Scand. 118, 346-352. https://doi.org/10.1111/j.1600-0463.2010.02602.x.

Mendez, F., Munoz, N., Posso, H., Molano, M., Moreno, V., van den Brule, A.J.C., Ronderos, M., Meijer, C., Munoz, A., Instituto Nacional de Cancerologia Human Papillomavirus Study Group, 2005. Cervical coinfection with human papillomavirus (HPV) types and possible implications for the prevention of cervical cancer by HPV vaccines. J. Infect. Dis. 192, 1158-1165. https://doi.org/10.1086/444391.

Merikukka, M., Kaasila, M., Namujju, P.B., Palmroth, J., Kirnbauer, R., Paavonen, J., Surcel, H.-M., Lehtinen, M., 2011. Differences in incidence and co-occurrence of vaccine and nonvaccine human papillomavirus types in Finnish population before human papillomavirus mass vaccination suggest competitive advantage for HPV33. Int. J. Cancer 128, 1114-1119. https://doi.org/10.1002/ijc.25675.

Mitchell, K.R., Mercer, C.H., Prah, P., Clifton, S., Tanton, C., Wellings, K., Copas, A., 2019. Why do men report more opposite-sex sexual partners than women? Analysis of the gender discrepancy in a British national probability survey. J. Sex Res. 56, 1-8. https://doi.org/10.1080/00224499.2018.1481193.

Monsonego, J., Bohbot, J.M., Pollini, G., Krawec, C., Vincent, C., Merignargues, I., Haroun, F., Sednaoui, P., Monfort, L., Dachez, R., Syrjänen, K., 2005. Performance of the Roche AMPLICOR human papillomavirus (HPV) test in prediction of cervical intraepithelial neoplasia (CIN) in women with abnormal PAP smear. Gynecol. Oncol. 99, 160-168. https://doi.org/10.1016/j.ygyno.2005.05.030.

Murall, C.L., McCann, K.S., Bauch, C.T., 2014. Revising ecological assumptions about human papillomavirus interactions and type replacement. J. Theor. Biol. 350, 98-109. https://doi.org/10.1016/j.tbi.2013.12.028.

Olsen, J., Jepsen, M.R., 2010. Human papillomavirus transmission and cost-effectiveness of introducing quadrivalent HPV vaccination in Denmark. Int. J. Technol. Assess. Health Care 26, 183-191. https://doi.org/10.1017/S0266462310000085.

Pinotti, F., Fleury, É., Guillemot, D., Boëlle, P.-Y., Poletto, C., 2019. Host contact dynamics shapes richness and dominance of pathogen strains. PLoS Comput. Biol. 15, e1006530 https://doi.org/10.1371/journal.pcbi.1006530.

Pons-Salort, M., Letort, V., Favre, M., Heard, I., Dervaux, B., Opatowski, L., Guillemot, D., 2013. Exploring individual HPV coinfections is essential to predict HPV-vaccination impact on genotype distribution: a model-based approach. Vaccine 31, 1238-1245. https://doi.org/10.1016/j.vaccine.2012.11.098.

Poolman, E.M., Elbasha, E.H., Galvani, A.P., 2008. Vaccination and the evolutionary ecology of human papillomavirus. Vaccine 26 (Suppl. 3), SC25-SC30. https://doi.org/10.1016/j.vaccine.2008.04.010.

Qu, W., Jiang, G., Cruz, Y., Chang, C.J., Ho, G.Y., Klein, R.S., Burk, R.D., 1997. PCR detection of human papillomavirus: comparison between MY09/MY11 and GP5+/GP6+ primer systems. J. Clin. Microbiol. 35, 1304-1310.

Ramanakumar, A.V., Naud, P., Roteli-Martins, C.M., de Carvalho, N.S., de Borba, P.C., Teixeira, J.C., Blatter, M., Moscicki, A.-B., Harper, D.M., Romanowski, B., Tyring, S. K., Ramjattan, B., Schuind, A., Dubin, G., Franco, E.L., HPV-007 Study Group, 2016. Incidence and duration of type-specific human papillomavirus infection in high-risk HPV-naïve women: results from the control arm of a phase II HPV-16/18 vaccine trial. BMJ Open 6, e011371. https://doi.org/10.1136/bmjopen-2016-011371.

Riethmuller, D., Gay, C., Bertrand, X., Bettinger, D., Schaal, J.P., Carbillet, J.P., Lassabe, C., Arveux, P., Seilles, E., Mougin, C., 1999. Genital human papillomavirus infection among women recruited for routine cervical cancer screening or for colposcopy determined by Hybrid Capture II and polymerase chain reaction. Diagn. Mol. Pathol. Am. J. Surg. Pathol. Part B 8, 157-164. https://doi.org/10.1097/00019606-199909000-00009.

Ryston, P., 2007. Profile likelihood for estimation and confidence intervals. Stata J. 7, 376-387.

Safaiean, M., Porras, C., Schiffman, M., Rodriguez, A.C., Wacholder, S., Gonzalez, P., Quint, W., van Doorn, L.-J., Sherman, M.E., Hxenseval, V., Herrero, R., Hildesheim, A., Costa Rican Vaccine Trial Group, 2010. Epidemiological study of anti-HPV16/18 seropositivity and subsequent risk of HPV16 and 18 infections. J. Natl. Cancer Inst. 102, 1653-1662. https://doi.org/10.1093/jnci/djq384.

Santé publique France, 2019. Cancer du col de l'utérus [WWW Document].
(https://www.santepubliquefrance.fr/maladies-et-traumatismes/cancers/can
cer-col-col-de-l-uterus/la-maladie/#tabs), (Accessed 7 June 2020).

Schneeberger, A., Mercer, C.H., Gregson, S.A.J., Ferguson, N.M., Nyamukapa, C.A., Anderson, R.M., Johnson, A.M., Garnett, G.P., 2004. Scale-free networks and sexually transmitted diseases: a description of observed patterns of sexual contacts in Britain and Zimbabwe. Sex. Transm. Dis. 31, 380-387. https://doi.org/10.1097/00007435-200406000-00012.

Serrano, B., Alemany, L., Tous, S., Bruni, L., Clifford, G.M., Weiss, T., Bosch, F.X., de Sanjosé, S., 2012. Potential impact of a nine-valent vaccine in human papillomavirus related cervical disease. Infect. Agent. Cancer 7, 38. https://doi.org/10.1186/1750-9378-7-38.

Shiboski, S., Padian, N.S., 1996. Population- and individual-based approaches to the design and analysis of epidemiologic studies of sexually transmitted disease transmission. J. Infect. Dis. 174, S188-S200. https://doi.org/10.1093/infdis/174. Supplement_2_S188.

Spinillo, A., Dal Bello, B., Alberizzi, P., Cesari, S., Gardella, B., Roccio, M., Silini, E.M., 2009. Clustering patterns of human papillomavirus genotypes in multiple infections. Virus Res. 142, 154-159. https://doi.org/10.1016/j.virusres.2009.02.004.

Tota, J.E., Ramanakumar, A.V., Jiang, M., Dillner, J., Walter, S.D., Kaufman, J.S., Coutlée, F., Villa, L.L., Franco, E.L., 2013. Epidemiologic approaches to evaluating the potential for human papillomavirus type replacement postvaccination. Am. J. Epidemiol. 178, 625-634. https://doi.org/10.1093/aje/kwt018.

Tota, J.E., Struyf, F., Sampson, J.N., Gonzalez, P., Ryser, M., Herrero, R., Schussler, J., Karkada, N., Rodriguez, A.C., Folschweiller, P., Porras, C., Schiffman, M., Schiller, J. T., Quint, W., Kreimer, A.R., Wheeler, C.M., Hildesheim, A., 2020. Efficacy of the AS04-adjuvanted HPV16/18 vaccine: pooled analysis of the Costa Rica vaccine and PATRICIA randomized controlled trials. JNCI J. Natl. Cancer Inst. 112, 818-828. https://doi.org/10.1093/jnci/djz222.

Trottier, H., Mahmud, S., Prado, J.C.M., Sobrinho, J.S., Costa, M.C., Rohan, T.E., Villa, L. L., Franco, E.L., 2008. Type-specific duration of human papillomavirus infection: implications for human papillomavirus screening and vaccination. J. Infect. Dis. 197, 1436-1447. https://doi.org/10.1086/587698.

Vaccarella, S., Söderlund-Strand, A., Franceschi, S., Plummer, M., Dillner, J., 2013. Patterns of human papillomavirus types in multiple infections: an analysis in women and men of the high throughput human papillomavirus monitoring study. PLoS One 8, e71617. https://doi.org/10.1371/journal.pone.0071617.

Van de Velde, N., Boily, M.-C., Drolet, M., Franco, E.L., Mayrand, M.-H., Kliewer, E.V., Coutlée, F., Laprise, J.-F., Malagón, T., Brisson, M., 2012. Population-level impact of the bivalent, quadrivalent, and nonavalent human papillomavirus vaccines: a model-based analysis. J. Natl. Cancer Inst. 104, 1712-1723. https://doi.org/10.1093/jnci/djs395.

Van de Velde, N., Brisson, M., Boily, M.-C., 2010. Understanding differences in predictions of HPV vaccine effectiveness: a comparative model-based analysis. Vaccine 28, 5473-5484. https://doi.org/10.1016/j.vaccine.2010.05.056.

Walboomers, J.M., Jacobs, M.V., Manos, M.M., Bosch, F.X., Kummer, J.A., Shah, K.V., Snijders, P.J., Peto, J., Meijer, C.J., Muñoz, N., 1999. Human papillomavirus is a necessary cause of invasive cervical cancer worldwide. J. Pathol. 189, 12-19. (https://doi.org/10.1002/(SICI)1096-9896(199909)189:13.00;2-F).

Xi, L.F., Edelstein, Z.R., Meyers, C., Ho, J., Cherne, S.L., Schiffman, M., 2009. Human papillomavirus types 16 and 18 DNA load in relation to coexistence of other types, particularly those in the same species. Cancer Epidemiol. Biomark. Prev. Publ. Am. Assoc. Cancer Res. Cosponsored Am. Soc. Prev. Oncol. 18, 2507-2512. https://doi.org/10.1158/1055-9965.EPI-09-0482.