# An agent-based model to simulate co-infection with human papillomaviruses in a partially vaccinated partnership network

Mélanie Bonneault, Maxime Flauder, Elisabeth Delarocque-Astagneau, Anne

Cm Thiébaut, Lulla Opatowski

► To cite this version:

Mélanie Bonneault, Maxime Flauder, Elisabeth Delarocque-Astagneau, Anne Cm Thiébaut, Lulla Opatowski. An agent-based model to simulate co-infection with human papillomaviruses in a partially vaccinated partnership network. 2021. inserm-03367886

HAL Id: inserm-03367886

https://inserm.hal.science/inserm-03367886

Submitted on 6 Oct 2021

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

An agent-based model to simulate co-infection with human papillo-
maviruses in a partially vaccinated partnership network

2

3

4

Mélanie Bonneault<sup>1,2,3</sup>, Maxime Flauder<sup>1,2</sup>, Elisabeth Delarocque-Astagneau<sup>2</sup>, Anne CM

5

Thiébaut<sup>3</sup> and Lulla Opatowski<sup>1,2</sup>

6

7

<sup>1</sup>Epidemiology and Modelling of Antibiotic Evasion Unit, Institut Pasteur, Paris, France

8

<sup>2</sup>Université Paris-Saclay, UVSQ, Inserm, CESP, 78180, Montigny-Le-Bretonneux, France

9

<sup>3</sup>Université Paris-Saclay, UVSQ, Inserm, CESP, 94807, Villejuif, France

10

11

Author for correspondence:

12

Mélanie Bonneault, Institut Pasteur, Epidemiology and Modelling of Antibiotic Evasion Unit,

13

25, rue du Docteur-Roux, 75475 Paris Cedex 15, France.

14

Email: melanie.bonneault@gmail.com

15

16

17

18
**Abstract.** Human papillomaviruses (HPV) are among the most common sexually transmitted infection and a necessary cause of cervical cancer. In the context of vaccination against a subgroup of genotypes, better understanding the role of biological interactions between HPV genotypes and social interactions between humans is essential to anticipate what the vaccine impact could be at the population level. Existing models that study interactions between genotypes are based on basic homogeneous assumptions about human-to-human transmission, whereas agent-based models that take into account the heterogeneity of human behavior and transmission do not consider possible interactions between genotypes. Here, we present a novel stochastic agent-based model formalizing the co-circulation on a human partnership network of multiple interacting genotypes, some of them being preventable by the vaccine (vaccine types) and others not. The model explicitly formalizes heterogeneity in sexual behaviors and makes it possible to explore distinct genotypic interaction mechanisms during intra-host co-infections. Using model simulations, we investigate infection dynamics in the population after vaccine introduction depending on assumptions regarding vaccine coverage and vaccine and non-vaccine genotype interactions.

34

**Keywords:** agent-based model, simulation, human papillomavirus, vaccination, transmission, genotype interactions.

37

38
39 **1 Introduction**

40 Mathematical models are useful tools to understand the dynamics of infectious disease
41 spread in human populations and predict the impact of control strategies such as screening or
42 immunization programs. Compartmental models are a classical form of models for studying
43 infection transmission, in which individuals are grouped indiscriminately in compartments ac-
44 cording to their infectious status (e.g., susceptible, infected, recovered). In the context of sex-
45 ually transmitted infections, the assumption of homogeneous sexual behaviors underlying
46 compartmental models is unrealistic, as informed by surveys showing high heterogeneity
47 across individuals and ages (Bajos and Bozon 2008). Agent-based models (ABMs) appear to be
48 the most suitable to reproduce the heterogeneity of individual behaviors and simulate their
49 effects on a population scale (Auchincloss and Diez Roux 2008). They also facilitate the mod-
50 elling of co-circulation of several pathogens, which is made complicated by the use of com-
51 partmental models.

52
53 In this paper, we present a model specifically applied to the simulation of between-human
54 transmissions of genital infections with human papillomavirus (HPV). HPV infections are
55 among the most common sexually transmitted infections, especially in younger people, and
56 could lead to the development of cancers around the age of 50 years (Walboomers et al. 1999;
57 de Sanjosé et al. 2007). Vaccination protecting against prolonged HPV infection was intro-
58 duced in 2007. However, available vaccines are not universal. They protect against a fraction
59 of HPV genotypes only (called vaccine types). Anticipating the impact of vaccination on the
global evolution of HPV infection at the population level therefore requires integrating as-
sumptions regarding how circulating genotypes, and in particular vaccine and non-vaccine
types, interact for host infection, which was not done before. We propose here a new ABM
combining a heterogenous contact network and between-genotypes interactions and apply it
to the French context.

The paper is organized as follows. After discussing existing works, we describe the specifica-
tions of the new agent-based model we developed. We then present results from experi-
mental simulations and finally provide a discussion on our model and results.

68 **2 Related works**

Only few models of HPV transmission have formalized potential interactions between geno-
types (Elbasha et al. 2008; Pons-Salort et al. 2013; Murall et al. 2014; Man et al. 2020). Elbasha
and Galvani studied the conditions of coexistence of two genotypes at equilibrium before vac-
cine introduction as a function of the baseline reproduction number (number of secondary
infections generated by the introduction of an infected individual) of each genotype (Elbasha
et al. 2008). Murall et al. questioned the evolutionary response of genotypes to the selection
pressure imposed by the immune response induced by vaccination (Murall et al. 2014). Pons-
Salort et al. examined the influence of co-infection dynamics (simultaneous or sequential elim-
ination of co-infections) on the evolution of non-vaccine prevalences after vaccination and
also addressed the issue of genotype replacement, i.e. an increased prevalence of genotypes
not included in the vaccine following vaccine introduction (Pons-Salort et al. 2013). In partic-
ular, they find that in the case of simultaneous clearance, genotypic replacement is always
observed regardless of the nature of the interaction between the genotypes. Finally, Man et
al. sought to determine whether type replacement could be observed in early post-vaccina-
tion surveillance (Man et al. 2020). However, these former models are all based on the as-
sumption that sexual contacts are homogeneous. It is likely that the diversity of behaviors has
a marked impact on the dynamics of co-circulation of genotypes at the population level since
it affects the risk of infection (Shiboski and Padian 1996). Therefore, faithful replication of
sexual behaviors is essential to correctly interpret observations of HPV prevalences and pro-
vide more accurate projections of the epidemiological and ecological consequences of vac-
cination. Among the possible modelling approaches, the ABM seems to be the most appropri-
ate to reproduce the heterogeneity of individual behaviors and to simulate their effects at the
population level (Auchincloss and Diez Roux 2008). Many ABMs have been developed follow-
ing the introduction of vaccines to study the effectiveness and cost-effectiveness of vaccina-
tion and/or screening (Olsen and Jepsen 2010; Van de Velde et al. 2012; Burger et al. 2018).
However, existing ABMs have not considered possible interactions between genotypes (Olsen
and Jepsen 2010; Van de Velde et al. 2012; Matthijsse et al. 2015; Burger et al. 2018; Johnson
et al. 2018). We therefore undertook to develop our own model and to adapt it to the French
context.

98 **3 Methods**

The model description follows the ODD (Overview, Design concepts, and Details) protocol as
99
defined by Grimm et al. (Grimm et al. 2006). We developed a multi-agent system to support
100
the process of HPV genotype transmission in the population. The global architecture of the
101
multi-agent system has been developed in C++ (version 4.9.0). In addition to classical elements
102
of models developed in object programming, the program includes two other components:
103
104 (1) a component for the agent (as detailed later), in which the agent's characteristics and related processes (to become infected, to start or to finish a partnership, to eliminate the infection, etc.) are defined; and (2) a component for the population (as detailed later), in which the variables of the environment of the agents are defined and which allows to initiate the processes (demography, partnership, infectious process, etc.) for all agents at each time step.
105
106
107
108
109 Simulations were run on the computational and storage services (TARS cluster) provided by the IT Department at Institut Pasteur, Paris. Statistical analyses and graphics were computed using R (version 3.5.2).
110
111

112 **3.1 Overview of the model**

113 **Purpose.** The ABM simulates a virtual population of individuals in heterosexual partnership and the co-circulation of distinct HPV genotypes in interaction. A vaccination campaign is also simulated, assuming that vaccinated individuals are protected against a portion of genotypes only (called "vaccine types", VT) whereas their infectious risk to the rest of the circulating types (called "non vaccine types", NVT) is not altered. The main goal is to understand how vaccination affects the global infection dynamics in the presence of between-genotype interactions on a heterogeneous partnership network.
114
115
116
117
118
119

120 **State variables and scales.** The agents are human individuals characterized by a set of varia- bles listed in Table 1 and described below. Each agent is explicitly modeled and characterized 121 by his/her sex, age, sexual activity class (number of partners he or she will have during the 122 year), partnership status (in three categories - virgin, with a partner, or single), vaccination 123
status, and infection status for each genotype *g* (in four categories - susceptible to (infection with) *g*, infected with *g*, naturally immune to *g*, or vaccinated if *g* is a VT).

<table><caption>Table 1. Individual variables characterizing each agent</caption><thead><tr><th>Variable</th><th>Type</th><th>Update time</th></tr></thead><tbody><tr><td>ID</td><td>Code number</td><td>At entrance</td></tr><tr><td>Sex</td><td>Boolean</td><td>At entrance</td></tr><tr><td>Age</td><td>Number weeks</td><td>At each time step</td></tr><tr><td>Vaccination status</td><td>Boolean</td><td>At entrance</td></tr><tr><td colspan="3"><strong>Variables related to partnership process</strong></td></tr><tr><td>Sexual activity class</td><td>Integer</td><td>At 15, 17, 19 and 24 years</td></tr><tr><td>No partner during the current year</td><td>Boolean</td><td>On the first day of each year</td></tr><tr><td>Partnership status</td><td>String</td><td>At the beginning and at the end of partnership</td></tr><tr><td>Partner ID</td><td>Pointer</td><td>At the beginning and at the end of partnership</td></tr><tr><td>End date of partnership or inactivity</td><td>Simulation week number</td><td>At the beginning and at the end of partnership</td></tr><tr><td>Counter of partners during the current year</td><td>Integer</td><td>At the onset of new partnership; and at reset on first day of each year</td></tr><tr><td>Counter of partners over the simulation</td><td>Integer</td><td>At the onset of new partnership</td></tr><tr><td colspan="3"><strong>Variables related to the infection process – for each genotype g</strong></td></tr><tr><td>Infection status</td><td>String</td><td>At the beginning and at the end of infection and immunity</td></tr><tr><td>End date for infection</td><td>Simulation week number</td><td>At the onset of infection</td></tr><tr><td>Natural immunity status</td><td>String</td><td>At the end of infection</td></tr><tr><td>End date for natural immunity</td><td>Simulation week number</td><td>At the end of infection</td></tr><tr><td>Counter of infections over the simulation</td><td>Integer</td><td>At the onset of infection</td></tr></tbody></table>

127

*Environmental variables.* Agents share a number of environmental variables characterizing the population as presented in Table 2 and detailed in the following sections. In most countries, anti-HPV vaccination was first recommended only for young females before the age of 15 years with possible catch-up generally until 18, with no vaccination of males (Baylor 2006;
132 ECDC 2020). We therefore assumed here that only females were vaccinated, and that vaccina-

133 tion occurred before their entrance in the model.

<table><caption>134 **Table 2.** Environmental variables characterizing the population and their values</caption><thead><tr><th>Description</th><th>Default value</th><th>Rationale</th></tr></thead><tbody><tr><td>Population size (N)</td><td>800,000 individuals</td><td>For computational purposes</td></tr><tr><td>Population age range</td><td>From 15 to 30 years of age</td><td>Fixed, corresponding to HPV infection peak</td></tr><tr><td colspan="3"><strong>Variables related to partnership process</strong></td></tr><tr><td>Age at first partnership</td><td>Normal distribution with mean 17.5 years in females, 16.8 in male</td><td>Mean, calibrated on (Bajos and Bozon 2008)</td></tr><tr><td>Distribution of individuals by sexual activity class</td><td>According to sex and age (Table 3 in appendix)</td><td>From (Bajos and Bozon 2008) and calibrated</td></tr><tr><td>Proportion of individuals without a partner during the year</td><td>According to sex and age (Table 3 in appendix)</td><td>From (Bajos and Bozon 2008) and fixed</td></tr><tr><td>Duration of partnership</td><td>According to sexual activity class and age (Table 3 in appendix)</td><td>Mean, from (Bajos and Bozon 2008) and calibrated</td></tr><tr><td>Duration of inactivity between two partnerships</td><td>According to sexual activity class (Table 3 in appendix)</td><td>Mean, calibrated</td></tr><tr><td>Frequency of sexual intercourse within a relationship</td><td>2 per week</td><td>Arbitrary</td></tr><tr><td>Duration of partner search before mixing between sexual activity groups</td><td>According to sexual activity class (Table 3 in appendix)</td><td>Calibrated</td></tr><tr><td>Maximum proportion of individuals changing from one extreme class to another at ages 17, 19 and 24</td><td>95% for females and 45% for males</td><td>Calibrated</td></tr><tr><td colspan="3"><strong>Variables related to infection process</strong></td></tr><tr><td>Number of genotypes</td><td>2 vaccine types (VT)<br/>12 non-vaccine types (NVT)</td><td>According to genotypes reported in (Markowitz et al. 2013)</td></tr><tr><td>Transmission probability parameter for VT and NVT genotypes</td><td>θ<sub>V</sub> = 0.16 per sexual act (VT)<br/>θ<sub>NV</sub> = 0.125 per sexual act (NVT)</td><td>Calibrated</td></tr><tr><td>Duration of infection</td><td>Exponential distribution with mean 52 weeks</td><td>From (Trottier et al. 2006)</td></tr><tr><td>Duration of immunity</td><td>Exponential distribution with mean 12 weeks</td><td>Calibrated</td></tr><tr><td>Genotype interaction parameter (γ)</td><td>1</td><td>To be varied</td></tr><tr><td>Initial probability of infection</td><td>0.2 for NVT and 0.08 for VT</td><td>Arbitrary</td></tr><tr><td>Initial probability of co-infection among those initially infected</td><td>0.4 for any first genotype g</td><td>Arbitrary</td></tr><tr><td>Date of vaccine introduction</td><td>After 70 years of simulation</td><td>To ensure prevalence equilibrium is reached</td></tr><tr><td>Vaccine coverage</td><td>60% of females by age cohort</td><td>To be varied</td></tr></tbody></table>
135 *Scales*. The model time step is the week. A simulation year lasts 52 weeks. Agents are not explicitly distributed in space.

137

138 **Process overview and scheduling.** The model proceeds in several processes which are sequentially updated on a weekly or yearly basis, depending on the processes, as detailed below.
139
140 Figure 1 summarizes the execution of the three major weekly processes in the model: demographic, partnership and infection processes.
141

142 *Demographic process*. Agents age one week with every time step. On the day they turn 30
143 years, agents leave the model. When an agent exits the population, he/she is directly replaced
144 by a new 15-year-old individual with an equal probability of being male or female, to ensure
145 the population remains stable over time. At inclusion the agent's states and counters are ini-
146 tialized. Sexual activity class and date of first partnership are drawn from probabilistic distri-
147 butions (see Table 2). In addition, if a vaccination campaign is in progress at a female's entry
148 date, she will be vaccinated according to a probability defined by vaccine coverage, protecting
149 her against VT genotypes with 100% efficacy for at least 15 years.

150 *Partnership process*. Every week, the list of agents available for a new relationship is updated;
151 each agent from this list then searches for a new partner according to his/her sexual activity
152 class and age. If the agent has been waiting on the list for a predefined period of time without
153 finding a partner within his/her sexual activity class, he/she is allowed to search within other
154 sexual activity classes. If a partner is found, a duration of partnership is drawn according to
155 the agent's sexual activity class and age (Table 3 in appendix). When the end date of a part-
156 nership is reached, the partnership is split. For each of the two agents involved, a duration of
157 inactivity is drawn according to his/her sexual activity class, defining the date for the agent to
158 become again available for a new partnership and be included in the list mentioned above.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic representation of the ordered weekly processes that could occur for each agent under conditions.</figcaption></figure>

In green the demographic process, in pink the partnership process and in orange the infec-
tion process
Infection process and natural disease history. At each time step, the model proceeds with infection transmission. For each genotype $g$, all partnerships characterized by one partner infected with $g$ and one partner susceptible to $g$ are evaluated: transmission of $g$ from the infected to the susceptible partner occurs with a transmission probability parameter $\beta_V$ or $\beta_{NV}$ depending on whether genotype $g$ belongs to VT or NVT group, respectively. If an infection occurs, the duration of infection is drawn from an exponential distribution (Table 2). In a second step, the model proceeds with the evolution of infection status for all infected agents. For any genotype $g$, all agents infected with $g$ who reach their date of recovery change their infection status: they become immune to genotype $g$ and cannot transmit genotype $g$ anymore.

Acquired natural immunity is assumed to confer total protection against genotype $g$ for a limited duration which is drawn from an exponential distribution (Table 2). All immunized agents reaching the end date of their immunized status become fully susceptible to genotype $g$ again.

The processes scheduled annually are executed on the first day of the year in the following order:

Reset of counters of annual number of partners. For each agent of the model, the counter of his/her cumulative number of partners over the current year is set to 0.

Change of sexual activity class. A portion of agents aged 17, 19 and 24 years and not engaged in long partnerships (i.e., >1 year) are randomly selected to change sexual activity class to conform to age groups available in the survey data and changing sexual behavior with age. Changing from one extreme class to the other is limited to a maximum proportion of individuals by sex (Table 2) to favor evolution between adjacent sexual activity classes. Numbers are
187 drawn to match the targeted distribution of sexual activity classes by age category. Addi-
188 ally, in all sexual activity classes, a fraction of agents who will have no partner over the coming
189 year is randomly selected (Table 2).

190 **3.2 Design Concept**

191 **Interaction**

192 *Inter-agent contacts.* Interactions between agents occur through heterosexual partnership.
193 The agent’s sexual activity class defines the number of partners the agent will have during the
194 year. Partnerships between agents from the same sexual activity class are promoted, although
195 mixing between classes is allowed after several weeks and unsuccessful trials to find a partner
196 within the same sexual activity class (Table 2). When a partnership is set, the two partners
197 involved are assumed to be in contacts at an arbitrary frequency of two sexual intercourses
198 per week over the partnership duration. Simultaneous partnership is not allowed in the
199 model.

200 *Inter-agent transmission.* An agent infected with a genotype *g* can transmit *g* at any time dur-
ing partnership with transmission probability *β*<sub>V</sub> or *β*<sub>NV</sub> per act, if his/her partner is not already
201 infected with *g* nor naturally immunized to *g*.

203 *Within-agent ecological interaction between genotypes.* When two genotypes simultaneously
204 infect an agent, interaction between genotypes may occur. This interaction is assumed to af-
205 fect the duration of infection with one of the genotypes. Interactions are defined unidirection-
206 ally: we only consider the impact of the presence of VT on the epidemiological characteristics
207 of NVT. Indeed, we assume that, if an agent already infected with a vaccine genotype *g*<sub>V</sub> is
subsequently infected with a non-vaccine genotype $g_{NV}$, the average duration of infection $D_{INF}^{eff}$ with $g_{NV}$ is modified by a multiplicative factor $\gamma$, called the strength of interaction (presented in Eqs 1 and in Fig. 2).

$$D_{INF}^{eff}(g_{NV}) = D_{INF}(g_{NV}) \times \gamma \quad (1)$$

Both competitive, whereby $\gamma<1$, i.e., the presence of VT reduces duration of infection with NVT, and synergistic interaction, whereby $\gamma>1$, i.e., the presence of VT enhances the duration of infection with NVT, are allowed. Since the period of infection coincides with the period of possible transmission, this also has an impact on the length of time the infected agent can transmit $g_{NV}$.

**Fig. 2.** Schematic representation of the interaction mechanism between VT and NVT during an agent co-infection.

Here the agent is first infected with VT and subsequently acquires NVT. The interaction alters the duration of infection with NVT by a multiplicative factor $\gamma$, the strength of interaction. The vertical arrows indicate varying recovery times depending on the interaction. C: competitive, N: neutral, S: synergistic

**Stochasticity.** The model is stochastic, stochasticity being included both through the definition of random parameters and through random processes. Sexual activity class, age at first partnership, age difference between partners, and all durations (of partnership, of inactivity be-
227 tween two partnerships, of infection, of immunization after vaccine introduction) are ran-
228 domly drawn from probability distributions (Table 2). Finding a partner, changing sexual activ-
229 ity class and getting infected are defined as random processes.

230 **Observations.** A large variety of data can be collected from model simulations. For the pur-
231 pose of this article, we record for every agent, at each time step, his/her sexual behavior (num-
232 ber of partners during the past year and lifetime) and his/her infection characteristics. Be-
233 cause we aim at monitoring the dynamics of infection at the population level, we aggregate
234 agent data according to human characteristics to obtain numbers of infected agents by age,
235 sex, sexual activity class, etc.

236 **3.3 Details on the processes**

237 **Initialization.** Model initialization is achieved in successive steps. First, all environmental var-
238 iables are initialized, as detailed in Table 2. Then all agents are created, with their partnership
239 and infection statuses and counters initialized. At initialization, sex and age of each agent are
240 randomly drawn to ensure homogeneous distribution of population across sex and age. Sexual
241 activity class and age at first partnership are drawn from probabilistic distributions (if the
242 agent is older than his/her age at first partnership, he/she can directly look for a partner).
243 Infection status with respect to VT and NVT is randomly drawn according to the probabilities
244 of infection and co-infection defined for initialization (Table 2). Those initial probabilities are
245 set to arbitrary non null values without any impact on prevalence equilibrium.
246 **4 Results**

247 **4.1 Data**

248 **Partnership data.** We used data reported in a French survey published in 2008 entitled *Contexte de la Sexualité en France* (CSF) (Bajos and Bozon 2008) to define relevant parameters
249 characterizing the partnership process and calibrated the model to reproduce realistic contact
250 patterns in the country. This population-based cross-sectional survey provides information on
251 sexual behaviors of a representative sample of 12 364 French males and females. We notably
252 used data on the distribution of individuals in sexual activity classes by age categories (18-19,
253 20-24 and 25-29 years old), the distribution of individuals without partner during the year by
254 age categories, age at first partnership, age difference between partners, duration of partner-
255 ship, and duration of inactivity between two partnerships.
256

257

258 **Infection data.** Because no data in the general population were available on the proportions
259 of individuals infected with HPV by age and genotype in the pre-vaccine era in France, we used
260 distributions reported from an epidemiological study carried out in the US (Markowitz et al.
261 2013), a country in which HPV epidemiology has been suggested to be similar to that in France
262 (de Sanjosé et al. 2007; Bruni et al. 2010). From (Markowitz et al. 2013) we extracted HPV
263 prevalence (i.e., the proportion of individuals infected in the population) of VT (gathering 2
264 genotypes: HPV 16 and 18) and NVT (gathering 12 genotypes: HPV 31, 33, 35, 39, 45, 51, 52,
265 56, 58, 59, 66 and 68) in the pre-vaccine era for females by age group (14-19, 20-24 and 25-
266 29 years old). We calibrated model parameters related to the infection process assuming no
267 genotypic interaction (strength $\gamma = 1$). These included transmission probability parameters $\delta_v$
268 and $\beta_{NV}$ and the mean duration of immunity for all genotypes (default values are presented in Table 2). For each interaction scenario, calibration was repeated for $\beta_{NV}$ in order to properly fit pre-vaccine NVT prevalence.
269
270

271 **4.2 Input and simulations**

272 The model was initialized with 800,000 individual agents interconnected through a contact network. We experiment 7 interaction scenarios (0.5, 0.7, 0.9, 1, 1.1, 1.3, 1.5) with 10 simulations each. The explored genotypic interaction values and the corresponding transmission rates are given as input to the model. First the simulations are run for 70 years to reach prevalence equilibrium in the pre-vaccination period. Vaccination is introduced at date 70 years, post-vaccination time is then defined as t=0 year, assuming 60% vaccine coverage. The model is simulated over 40 years post-vaccine introduction. The execution of a single simulation over 110 years takes about 2 hours. The 10 independent stochastic simulations of each scenario are generated thanks to a bash code taking the explored parameter values as input and run in parallel on the cluster nodes. Sensitivity analyses were also carried out to evaluate the impact of our model assumptions on observed results.
273
274
275
276
277
278
279
280
281
282

283 **4.3 Output**

284 Fig. 3 displays the dynamics of VT and NVT prevalence before and after vaccine introduction for 10 independent repetitions of the model in seven distinct scenarios of between-genotype interaction.
285
286

287 **Vaccine impact on VT prevalence.** Before vaccine introduction, the model was able to reproduce realistic VT prevalence patterns, with an average 10.7% in females. The prevalence was
288
slightly higher in the age group 20-24 years old than in the other groups and lowest in 15-19 year olds. After vaccine introduction, an expected decrease of VT prevalence was observed in females, which stabilized to a new equilibrium (3.2%) once the age cohort had had a chance to be immunized (after 5 years among 15-19 year olds up to 15 years among 25-29 year olds).
All calibrated interaction scenarios yielded undistinguishable patterns of VT prevalence and stochasticity was weak, as demonstrated by the limited variability across the 10 simulations displayed for VT prevalence in Fig. 3.

<figure><img src="image_3.png" /><figcaption>Fig. 3. Temporal evolution of the prevalence of NVT and VT genotypes according to interaction scenarios</figcaption></figure>

Proportion of infected of NVT (upper plain line curves) and VT (lower dashed curves) genotypes according to interaction scenarios (colors) (A) All females; (B) Females aged 15 to 19 years; (C) Females aged 20 to 24 years; (D) Females aged 25 to 29 years. The vaccine is offered from t = 0 to all females <15 years of age with 60% vaccine coverage. For each interaction scenario, results of 10 simulations are displayed.
303 **Vaccine impact on NVT prevalence.** The model also satisfactorily reproduced NVT prevalence 304 before vaccine introduction, reaching an average 26.5% in females. The same pattern of NVT 305 prevalence was seen across age groups as with VT. Minor differences were visible across cali- 306 brated interaction scenarios. Stochasticity was more pronounced for NVT than for VT geno- 306 types, particularly in the 20-24 age group, but standard deviations remained limited to <4% in 307 all scenarios. 308

309 After vaccine introduction, the evolution of NVT prevalence strongly depended on the inter- 310 action scenario (Fig. 2). In the absence of genotypic interaction (neutral scenario), NVT preva- 311 lence remained constant over time, while in non-neutral interaction scenarios some dynamics 312 was observed at the population level. Under competitive interaction (y<1), NVT prevalence 313 increased while, conversely, under synergistic interaction (y>1), NVT prevalence decreased 314 over time until reaching a new equilibrium. The magnitude of prevalence variations depended 315 on the strength of interaction: the more it deviated from 1, the more NVT prevalence deviated 316 from pre-vaccine equilibrium. The magnitude of the increase or decrease also appeared to be 317 proportional to NVT prevalence before vaccine introduction. Indeed, variations in NVT preva- 318 lence appeared to be most pronounced in 20-24 year old females who also had the highest 319 prevalence before vaccine introduction. Times to reach the new prevalence equilibrium was 320 the same for NVT and VT. 321

321 **Sensitivity analyses.** With respect to the partnership process, combinations of parameter val- 322 ues that best minimized the calibration criteria after the chosen ones resulted in VT and NVT 323 prevalence curves changing at the same time and with comparable magnitude. For the infec-
324 tion process, modifying one of the parameter values, whether it is the mean duration of infec-
325 tion, the mean duration of immunity or the transmission probabilities of the VT and NVT, af-
326 fects VT and NVT prevalence equilibrium quantitatively but not the overall qualitative dynam-
327 ics following the introduction of the vaccine.

328 5 Discussion

329 In this paper, we introduce a new agent-based model simulating the co-circulation of vaccine
330 and non-vaccine genotypes in a heterogeneous human population. The main originality of this
331 ABM is to consider genotype interactions. The agent-based approach makes it possible to de-
332 fine precisely natural history of HPV infection as well as co-infections in individuals and to
333 assess vaccine impact on the overall dynamics of HPV infection at the population level.

334 The illustration that we present is an example of how the model can be used to simulate VT
335 and NVT prevalence according to genotypic interaction scenarios and prevention measures
336 such as vaccination. Successively calibrated interaction scenarios help investigate how syner-
337 gistic interactions between genotypes could contribute to the decrease in overall HPV preva-
338 lence aimed by vaccine introduction. Conversely, we show that competitive interactions may
339 lead to substantial increase in NVT prevalence: because VT prevalence decreases following
340 vaccination introduction, so does the prevalence of VT and NVT coinfections, resulting in
341 longer duration of infection with NVT and more opportunities for NVT transmission. This eco-
342 logical replacement can potentially minimize the decrease in or even increase overall HPV
343 prevalence, therefore limiting vaccine benefit at the population level.
Although vaccination coverage is much lower in France (around 20%), we have chosen to present here a scenario with 60% vaccination coverage for a better observation of the replacement phenomenon. With lower vaccination coverage, variations in NVT prevalence are still visible under non-neutral interactions, but with reduced magnitude (results not shown).

We also chose a fixed population of 800,000 individuals, a size typical of a big metropolitan area such as Paris, to avoid additional complexity of modelling spatial constraints on contact mixing. Considering such a large population yields more stable results and allows us to investigate subpopulations such as age groups in Fig. 3. Due to strong heterogeneity in contacts across ages and sexual activity classes, simulating smaller populations may lead to more variability and possible extinctions of the infectious process, making any calibration difficult. Simulating the real French population size (~9,200,000 individuals 18-30 years of age) would require extremely long calculation times and would only slightly reduce the already limited variability we observe in our results.

Despite a detailed agent-based approach is used here, some simplifications and assumptions were made to compensate for the lack of data and gaps in knowledge. Several limitations should be underlined.

First, it was not possible to calibrate the model simultaneously to sexual behaviors reported by males and females, as distributions of numbers of partners did not match. In the CSF survey, more males declared large numbers of partners than females (Bajos and Bozon 2008). Reporting bias differential by sex has been observed in previous studies (Fenton et al. 2001; Mitchell et al. 2019) and may explain part of our calibration difficulties. To overcome those difficulties, we focused on females' data in order to adapt the model and the estimation parameters, while
trying to make it consistent with males' data in terms of cumulative number of partners. Moreover, because our interest was to simulate realistic heterogeneous contact patterns than can reproduce the typical bell-shaped curve of HPV prevalence according to age, we did not consider any temporal variation in sexual behaviors. Although the model was run for decades assuming unrealistically stable behaviors, that step only aimed at reaching prevalence equilibrium before the introduction of the vaccination and further analysis of the post vaccination dynamics was restricted to 40 years.

Second, to keep our model relatively simple, only two transmission probability parameters were estimated for VT and NVT genotypes, thus mimicking average dynamics across the two VT and the 12 NVT genotypes, respectively. However, large differences in prevalence are reported not only between VT and NVT groups, but also within each group (Markowitz et al. 2013). It could therefore be relevant to simulate VT and NVT genotypes more finely by calibrating transmission probability parameters by genotype. It would then be possible to study each genotype individually as well as to simulate interactions more precisely, for example between certain defined genotypes.

Finally, simulation results presented here were restricted to the assumption that genotypic interaction would affect the infection duration of a second acquired virus. Interaction mechanisms between HPV genotypes remain unknown and could as well involve viral load and/or cell-infection ability (McLaughlin-Drubin and Meyers 2004; Xi et al. 2009; Biryukov and Meyers 2018). Alternative mechanism interaction scenarios were therefore incorporated in the model code, affecting probability of acquisition. Such mechanisms could be explored in the future to measure their population-level impact on prevalence and compare their credibility.
388 In conclusion, our model formalizing both the heterogeneity of sexual behaviors at the popu-
389 lation level and the co-circulation of distinct groups of genotypes allows us to evaluate path-
390 ogen transmission dynamics in realistic sexual networks. In the context of HPV transmission,
391 our simulation results show that interactions between HPV genotypes can significantly impact
392 effectiveness of anti-HPV vaccination at reducing HPV prevalence. This model can be used to
393 further investigate the impact of most recent anti-HPV vaccines targeting a larger number of
394 genotypes, again altering prevalence equilibrium, as well as the recent extension of vaccina-
395 tion to males (starting in 2021 in France). Vaccination of boys is expected to increase the mag-
396 nitude and speed of change in the prevalence of VT and NVT. Finally, the model could easily
397 be adapted to study other research questions related to sexually transmitted pathogens and
398 their interactions.
399

399

400

401

402

403

404

405

406

407

408

409

410
411 **Abbreviations**

412 ABM(s), Agent-based model(s); HPV, human papillomavirus; ODD, Overview, Design concepts,

413 and Details; VT, vaccine type; NVT, non-vaccine type.

414 **Availability of data and materials**

415 All the source codes of the model and those for the analyses are available on request.

416 **Competing interests**

417 LO reports research grants from Pfizer through her institution, unrelated to the submitted

418 work. All other authors report no competing interests.

419 **Funding**

420 MB was funded by the INCa [grant DOC 2017-123] and MGEN, and her work was supported

421 by internal resources of Institut Pasteur, the French National Institute of Health and Medical

422 Research (Inserm) and the University of Versailles Saint-Quentin-en-Yvelines (UVSQ).

423 **Authors' contributions**

424 MB, MF, EDA, ACMT and LO conceptualized the project. ACMT and LO supervised the project.

425 MB, ACMT and LO designed the model. MB and MF developed the model. MB, ACMT and LO

426 performed the validation of the model and the experimental simulations. All the authors re-

427 viewed the article and approved the submission.

428 **Acknowledgements**

429 We thank Margarita Pons-Salort for helpful discussions and ideas provided early during the

430 project.
431 Appendix

<table><caption>Table 3. Complementary values of variables related to the partnership process according to sexual activity class and possibly age and sex</caption><thead><tr><th>Sex</th><th>Age (years)</th><th colspan="4">Number of partners per year</th></tr><tr><th></th><th></th><th>1</th><th>2-3</th><th>≥4</th><th>0<sup>a</sup></th></tr></thead><tbody><tr><td colspan="6">Distribution of individuals by sexual activity class (%)</td></tr><tr><td rowspan="4">Female</td><td>15-17</td><td>85.00<br>(calibrated)<sup>b</sup></td><td>12.00<br>(calibrated)</td><td>3.00<br>(calibrated)</td><td>3.63<br>(fixed)</td></tr><tr><td>18-19</td><td>68.99</td><td>24.72</td><td>6.29</td><td>3.63</td></tr><tr><td>20-24</td><td>75.84</td><td>21.56</td><td>2.60</td><td>7.24</td></tr><tr><td>25-29</td><td>85.75</td><td>11.07</td><td>3.18</td><td>5.59</td></tr><tr><td rowspan="4">Male</td><td>15-17</td><td>75.00<br>(calibrated)</td><td>19.00<br>(calibrated)</td><td>6.00<br>(calibrated)</td><td>14.38<br>(fixed)</td></tr><tr><td>18-19</td><td>60.42</td><td>27.59</td><td>11.99</td><td>14.38</td></tr><tr><td>20-24</td><td>61.85</td><td>28.15</td><td>10.00</td><td>10.89</td></tr><tr><td>25-29</td><td>76.30</td><td>17.23</td><td>6.47</td><td>6.69</td></tr><tr><td colspan="6">Duration of partnership</td></tr><tr><td rowspan="3">Any</td><td>15-19</td><td>4 categories (&lt;1, 1-4, 5-9 and &gt;10 years) distributed as in (Bajos and Bozon 2008)<sup>c</sup> (mean 224 weeks)</td><td rowspan="3">Gamma distribution with calibrated mean 14 weeks (fixed variance 8)</td><td rowspan="3">Gamma distribution with calibrated mean 4 weeks (fixed variance 5)<sup>d</sup></td><td rowspan="3">Not relevant</td></tr><tr><td>20-24</td><td>4 categories distributed as in (Bajos and Bozon 2008)<sup>c</sup> (mean 264 weeks)</td></tr><tr><td>25-29</td><td>4 categories distributed as in (Bajos and Bozon 2008)<sup>c</sup> (mean 440 weeks)</td></tr><tr><td colspan="6">Duration of inactivity between two partnerships</td></tr><tr><td>Any</td><td>All</td><td>Uniform distribution [1:104] weeks</td><td>Gamma distribution with calibrated mean 14 weeks (fixed variance 8)</td><td>Gamma distribution with calibrated mean 4 weeks (fixed variance 5)<sup>d</sup></td><td>At least 52 weeks</td></tr><tr><td colspan="6">Duration of partner search before mixing between sexual activity classes</td></tr><tr><td>Any</td><td>All</td><td>10 weeks (calibrated)</td><td>5 weeks (calibrated)</td><td>30 weeks (calibrated)</td><td>Not relevant</td></tr></tbody></table>

<footnote><sup>a</sup> Individuals without a partner are randomly selected each year from the population;</footnote>

<footnote><sup>b</sup> When not specified fixed or calibrated, values are extracted from (Bajos and Bozon 2008);</footnote>

<footnote><sup>c</sup> One category of partnership duration is drawn from a multinomial distribution whose parameters are the proportions reported in (Bajos and Bozon 2008), then a precise value is drawn from a uniform distribution within the duration category;</footnote>

<footnote><sup>d</sup> The duration of partnership and between two partnerships are identical in most sexually active individuals and drawn only once a year.</footnote>
441 References

442 Auchincloss AH, Diez Roux AV (2008) A new tool for epidemiology: the usefulness of dynamic-agent models in understanding place effects on health. Am J Epidemiol 168(1):1-8.
443
444 https://doi.org/10.1093/aje/kwn118

445 Bajos N, Bozon M (2008) Sexualité, genre et santé : les apports de l'enquête « Contexte de la sexualité en France ». La Découverte
446

447 Baylor NW (2006) Approved Products - June 8, 2006 Approval Letter - Human Papillomavirus Quadrivalent (Types 6, 11, 16, 18) Vaccine, Recombinant. http://wayback.archive-it.org/7993/20170722145339/https://www.fda.gov/BiologicsBloodVaccines/Vaccines/ApprovedProducts/ucm111283.htm. Accessed 12 Dec 2020

451 Biryukov J, Meyers C (2018) Superinfection Exclusion between Two High-Risk Human Papillomavirus Types during a Coinfection. J Virol 92(8):e01993-17.
452
453 https://doi.org/10.1128/JVI.01993-17

454 Bruni L, Diaz M, Castellsagué X, Ferrer E, Bosch FX, 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(12):1789-1799. https://doi.org/10.1086/657321

457 Burger EA, Campos NG, Sy S, Regan C, Kim JJ (2018) Health and economic benefits of single-dose HPV vaccination in a Gavi-eligible country. Vaccine 36(32 Pt A):4823-4829.
458
459 https://doi.org/10.1016/j.vaccine.2018.04.061

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

464 ECDC (2020) Guidance on HPV vaccination in EU countries: focus on boys, people living with HIV and 9-valent HPV vaccine introduction. In: European Centre for Disease Prevention and Control. https://www.ecdc.europa.eu/en/publications-data/guidance-hpv-vaccination-eu-focus-boys-people-living-hiv-9vHPV-vaccine. Accessed 2 Jul 2020

468 Elbasha EH, Dasbach EJ, Insinga RP (2008) A multi-type HPV transmission model. Bull Math Biol 70(8):2126-2176. https://doi.org/10.1007/s11538-008-9338-x

470 Fenton KA, Johnson AM, McManus S, Erens B (2001) Measuring sexual behaviour: methodological challenges in survey research. Sex Transm Infect 77(2):84-92.
471
472 https://doi.org/10.1136/sti.77.2.84

473 Grimm V, Berger U, Bastiansen F, Eliassen S, Ginot V, Giske J, Goss-Custard J, Grand T, Heinz SK, Huse G, Huth A, Jepsen JU, Jørgensen C, Mooij WM, Müller B, Pe'er G, Piou C, Railsback SF, Robbins AM, Robbins MM, Rossmanith E, Rüger N, Strand E, Souissi S, Stillman RA, Vabø R, Visser U, DeAngelis DL (2006) A standard protocol for describing individual-
477 based and agent-based models. Ecological Modelling 198(1):115-126.
478 https://doi.org/10.1016/j.ecolmodel.2006.04.023

479 Johnson HC, Lafferty EI, Eggo RM, Louie K, Soldan K, Waller J, Edmunds WJ (2018) Effect of HPV vaccination and cervical cancer screening in England by ethnicity: a modelling study. Lancet Public Health 3(1):e44-e51. https://doi.org/10.1016/S2468-2667(17)30238-4

483 Man I, Vänskä S, Lehtinen M, Bogaards JA (2020) Human papillomavirus genotype replacement: still too early to tell? J Infect Dis. https://doi.org/10.1093/infdis/jiaa032

485 Markowitz LE, Hariri S, Lin C, Dunne EF, Steinau M, McQuillan G, Unger ER (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(3):385-393. https://doi.org/10.1093/infdis/jit192

489 Matthijsse SM, van Rosmalen J, Hontelez JAC, Bakker R, de Kok IMCM, van Ballegooijen M, de Vlas SJ (2015) The role of acquired immunity in the spread of human papillomavirus (HPV): explorations with a microsimulation model. PLoS One 10(2):e0116618.
491
492 https://doi.org/10.1371/journal.pone.0116618

493 McLaughlin-Drubin ME, Meyers C (2004) Evidence for the coexistence of two genital HPV types within the same host cell in vitro. Virology 321(2):173-180.
494
495 https://doi.org/10.1016/j.virol.2003.12.019

496 Mitchell KR, Mercer CH, 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):1-8.
497
498
499 https://doi.org/10.1080/00224499.2018.1481193

500 Murall CL, McCann KS, Bauch CT (2014) Revising ecological assumptions about human papillomavirus interactions and type replacement. J Theor Biol 350:98-109.
501
502 https://doi.org/10.1016/j.jtbi.2013.12.028

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

506 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(8):1238-1245.
507
508
509 https://doi.org/10.1016/j.vaccine.2012.11.098

510 Shiboski S, Padian NS (1996) Population- and individual-based approaches to the design and analysis of epidemiologic studies of sexually transmitted disease transmission. J Infect Dis 174(Supplement 2):S188-S200. https://doi.org/10.1093/infdis/174.Supplement_2.S188
514 Trottier H, Mahmud S, Costa MC, Sobrinho JP, Duarte-Franco E, Rohan TE, Ferenczy A, Villa LL,
515 Franco EL (2006) Human papillomavirus infections with multiple types and risk of cer-
516 vical neoplasia. Cancer Epidemiol Biomarkers Prev 15(7):1274-1280.
517 https://doi.org/10.1158/1055-9965.EPI-06-0129

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

522 Walboomers JM, Jacobs MV, Manos MM, Bosch FX, Kummer JA, Shah KV, Snijders PJ, Peto J,
523 Meijer CJ, Muñoz N (1999) Human papillomavirus is a necessary cause of invasive cer-
524 vical cancer worldwide. J Pathol 189(1):12-19. https://doi.org/10.1002/(SICI)1096-
525 9896(199909)189:1<12::AID-PATH431>3.0.CO;2-F

526 Xi LF, Edelstein ZR, Meyers C, Ho J, Cherne SL, Schiffman M (2009) Human papillomavirus types
527 16 and 18 DNA load in relation to coexistence of other types, particularly those in the
528 same species. Cancer Epidemiol Biomarkers Prev 18(9):2507-2512.
529 https://doi.org/10.1158/1055-9965.EPI-09-0482