Mélanie Bonneault
email: melanie.bonneault@gmail.com

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

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

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 intrahost 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.

Introduction

Mathematical models are useful tools to understand the dynamics of infectious disease spread in human populations and predict the impact of control strategies such as screening or immunization programs. Compartmental models are a classical form of models for studying infection transmission, in which individuals are grouped indiscriminately in compartments according to their infectious status (e.g., susceptible, infected, recovered). In the context of sexually transmitted infections, the assumption of homogeneous sexual behaviors underlying compartmental models is unrealistic, as informed by surveys showing high heterogeneity across individuals and ages [START_REF] Bajos | Sexualité, genre et santé : les apports de l'enquête « Contexte de la sexualité en France[END_REF]. Agent-based models (ABMs) appear to be the most suitable to reproduce the heterogeneity of individual behaviors and simulate their effects on a population scale [START_REF] Auchincloss | A new tool for epidemiology: the usefulness of dynamicagent models in understanding place effects on health[END_REF]. They also facilitate the modelling of co-circulation of several pathogens, which is made complicated by the use of compartmental models.

In this paper, we present a model specifically applied to the simulation of between-human transmissions of genital infections with human papillomavirus (HPV). HPV infections are among the most common sexually transmitted infections, especially in younger people, and could lead to the development of cancers around the age of 50 years [START_REF] Walboomers | Human papillomavirus types 16 and 18 DNA load in relation to coexistence of other types, particularly those in the same species[END_REF][START_REF] De Sanjosé | Worldwide prevalence and genotype distribution of cervical human papillomavirus DNA in women with normal cytology: a meta-analysis[END_REF]. Vaccination protecting against prolonged HPV infection was introduced in 2007. However, available vaccines are not universal. They protect against a fraction 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 assumptions 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 specifications of the new agent-based model we developed. We then present results from experimental simulations and finally provide a discussion on our model and results.

Related works

Only few models of HPV transmission have formalized potential interactions between genotypes [START_REF] Elbasha | A multi-type HPV transmission model[END_REF][START_REF] Pons-Salort | Exploring individual HPV coinfections is essential to predict HPV-vaccination impact on genotype distribution: a model-based approach[END_REF][START_REF] Murall | Revising ecological assumptions about human papillomavirus interactions and type replacement[END_REF][START_REF] Man | Human papillomavirus genotype replacement: still too early to tell?[END_REF]). Elbasha and Galvani studied the conditions of coexistence of two genotypes at equilibrium before vaccine introduction as a function of the baseline reproduction number (number of secondary infections generated by the introduction of an infected individual) of each genotype [START_REF] Elbasha | A multi-type HPV transmission model[END_REF]. Murall et al. questioned the evolutionary response of genotypes to the selection pressure imposed by the immune response induced by vaccination [START_REF] Murall | Revising ecological assumptions about human papillomavirus interactions and type replacement[END_REF]. Pons-Salort et al. examined the influence of co-infection dynamics (simultaneous or sequential elimination 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 [START_REF] Pons-Salort | Exploring individual HPV coinfections is essential to predict HPV-vaccination impact on genotype distribution: a model-based approach[END_REF]. In particular, 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-vaccination surveillance [START_REF] Man | Human papillomavirus genotype replacement: still too early to tell?[END_REF]). However, these former models are all based on the assumption 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 [START_REF] Shiboski | Population-and individual-based approaches to the design and analysis of epidemiologic studies of sexually transmitted disease transmission[END_REF]. Therefore, faithful replication of sexual behaviors is essential to correctly interpret observations of HPV prevalences and provide more accurate projections of the epidemiological and ecological consequences of vaccination. Among the possible modelling approaches, the ABM seems to be the most appropriate to reproduce the heterogeneity of individual behaviors and to simulate their effects at the population level [START_REF] Auchincloss | A new tool for epidemiology: the usefulness of dynamicagent models in understanding place effects on health[END_REF]. Many ABMs have been developed following the introduction of vaccines to study the effectiveness and cost-effectiveness of vaccination and/or screening [START_REF] Olsen | Human papillomavirus transmission and cost-effectiveness of introducing quadrivalent HPV vaccination in Denmark[END_REF][START_REF] Van De Velde | Population-level impact of the bivalent, quadrivalent, and nonavalent human papillomavirus vaccines: a model-based analysis[END_REF][START_REF] Burger | Health and economic benefits of singledose HPV vaccination in a Gavi-eligible country[END_REF]).

However, existing ABMs have not considered possible interactions between genotypes (Olsen and Jepsen 2010; [START_REF] Van De Velde | Population-level impact of the bivalent, quadrivalent, and nonavalent human papillomavirus vaccines: a model-based analysis[END_REF][START_REF] Matthijsse | The role of acquired immunity in the spread of human papillomavirus (HPV): explorations with a microsimulation model[END_REF][START_REF] Burger | Health and economic benefits of singledose HPV vaccination in a Gavi-eligible country[END_REF][START_REF] Johnson | Effect of HPV vaccination and cervical cancer screening in England by ethnicity: a modelling study[END_REF].We therefore undertook to develop our own model and to adapt it to the French context.

Methods

The model description follows the ODD (Overview, Design concepts, and Details) protocol as defined by [START_REF] Grimm | A standard protocol for describing individual-based and agent-based models[END_REF]. We developed a multi-agent system to support the process of HPV genotype transmission in the population. The global architecture of the multi-agent system has been developed in C++ (version 4.9.0). In addition to classical elements of models developed in object programming, the program includes two other components:

(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.

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).

Overview of the model

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.

State variables and scales. The agents are human individuals characterized by a set of variables listed in ECDC 2020). We therefore assumed here that only females were vaccinated, and that vaccination occurred before their entrance in the model. 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.

Figure 1 summarizes the execution of the three major weekly processes in the model: demographic, partnership and infection processes. 3 in appendix). When the end date of a partnership is reached, the partnership is split. For each of the two agents involved, a duration of inactivity is drawn according to his/her sexual activity class, defining the date for the agent to become again available for a new partnership and be included in the list mentioned above. Fig. 1. Schematic representation of the ordered weekly processes that could occur for each agent under conditions. In green the demographic process, in pink the partnership process and in orange the infection process

Demographic
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 βV or β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 drawn to match the targeted distribution of sexual activity classes by age category. Additionally, in all sexual activity classes, a fraction of agents who will have no partner over the coming year is randomly selected (Table 2).

Design Concept Interaction

Inter-agent contacts. Interactions between agents occur through heterosexual partnership.

The agent's sexual activity class defines the number of partners the agent will have during the year. Partnerships between agents from the same sexual activity class are promoted, although mixing between classes is allowed after several weeks and unsuccessful trials to find a partner within the same sexual activity class (Table 2). When a partnership is set, the two partners involved are assumed to be in contacts at an arbitrary frequency of two sexual intercourses per week over the partnership duration. Simultaneous partnership is not allowed in the model.

Inter-agent transmission. An agent infected with a genotype g can transmit g at any time during partnership with transmission probability V or NV per act, if his/her partner is not already infected with g nor naturally immunized to g.

Within-agent ecological interaction between genotypes. When two genotypes simultaneously infect an agent, interaction between genotypes may occur. This interaction is assumed to affect the duration of infection with one of the genotypes. Interactions are defined unidirectionally: we only consider the impact of the presence of VT on the epidemiological characteristics of NVT. Indeed, we assume that, if an agent already infected with a vaccine genotype gV is subsequently infected with a non-vaccine genotype gNV, the average duration of infection 𝐷 𝐼𝑁𝐹 𝑒𝑓𝑓 with gNV is modified by a multiplicative factor γ, called the strength of interaction (presented in Eqs 1 and in Fig. 2).

𝐷 𝐼𝑁𝐹 𝑒𝑓𝑓 (𝑔 𝑁𝑉 ) = 𝐷 𝐼𝑁𝐹 (𝑔 𝑁𝑉 ) × 𝛾 (1)
Both competitive, whereby γ<1, i.e., the presence of VT reduces duration of infection with NVT, and synergistic interaction, whereby γ>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 gNV.

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 𝛾, 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-tween two partnerships, of infection, of immunization after vaccine introduction) are randomly drawn from probability distributions (Table 2). Finding a partner, changing sexual activity class and getting infected are defined as random processes.

Observations.

A large variety of data can be collected from model simulations. For the purpose of this article, we record for every agent, at each time step, his/her sexual behavior (number of partners during the past year and lifetime) and his/her infection characteristics. Because we aim at monitoring the dynamics of infection at the population level, we aggregate agent data according to human characteristics to obtain numbers of infected agents by age, sex, sexual activity class, etc.

Details on the processes

Initialization. Model initialization is achieved in successive steps. First, all environmental variables are initialized, as detailed in Infection status with respect to VT and NVT is randomly drawn according to the probabilities of infection and co-infection defined for initialization (Table 2). Those initial probabilities are set to arbitrary non null values without any impact on prevalence equilibrium.

Results

Data

Partnership data. We used data reported in a French survey published in 2008 entitled Contexte de la Sexualité en France (CSF) [START_REF] Bajos | Sexualité, genre et santé : les apports de l'enquête « Contexte de la sexualité en France[END_REF] HPV 31,33,35,39,45,51,52,56,58,59,66 and 68) in the pre-vaccine era for females by age group (14-19, 20-24 and 25-29 years old). We calibrated model parameters related to the infection process assuming no genotypic interaction (strength  = 1). These included transmission probability parameters βV and β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 NV in order to properly fit pre-vaccine NVT prevalence.

Input and simulations

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.

Output

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.

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 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. Vaccine impact on NVT prevalence. The model also satisfactorily reproduced NVT prevalence before vaccine introduction, reaching an average 26.5% in females. The same pattern of NVT prevalence was seen across age groups as with VT. Minor differences were visible across calibrated interaction scenarios. Stochasticity was more pronounced for NVT than for VT genotypes, particularly in the 20-24 age group, but standard deviations remained limited to <4% in all scenarios.

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

Sensitivity analyses. With respect to the partnership process, combinations of parameter values that best minimized the calibration criteria after the chosen ones resulted in VT and NVT prevalence curves changing at the same time and with comparable magnitude. For the infec-tion process, modifying one of the parameter values, whether it is the mean duration of infection, the mean duration of immunity or the transmission probabilities of the VT and NVT, affects VT and NVT prevalence equilibrium quantitatively but not the overall qualitative dynamics following the introduction of the vaccine.

Discussion

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

The illustration that we present is an example of how the model can be used to simulate VT and NVT prevalence according to genotypic interaction scenarios and prevention measures such as vaccination. Successively calibrated interaction scenarios help investigate how synergistic interactions between genotypes could contribute to the decrease in overall HPV prevalence aimed by vaccine introduction. Conversely, we show that competitive interactions may lead to substantial increase in NVT prevalence: because VT prevalence decreases following vaccination introduction, so does the prevalence of VT and NVT coinfections, resulting in longer duration of infection with NVT and more opportunities for NVT transmission. This ecological replacement can potentially minimize the decrease in or even increase overall HPV 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 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 [START_REF] Bajos | Sexualité, genre et santé : les apports de l'enquête « Contexte de la sexualité en France[END_REF]. Reporting bias differential by sex has been observed in previous studies [START_REF] Fenton | Measuring sexual behaviour: methodological challenges in survey research[END_REF][START_REF] Mitchell | Why do men report more opposite-sex sexual partners than women? Analysis of the gender discrepancy in a british national probability survey[END_REF]) 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 [START_REF] Markowitz | 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[END_REF]. 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 [START_REF] Mclaughlin-Drubin | Evidence for the coexistence of two genital HPV types within the same host cell in vitro[END_REF][START_REF] Walboomers | Human papillomavirus types 16 and 18 DNA load in relation to coexistence of other types, particularly those in the same species[END_REF][START_REF] Biryukov | Superinfection Exclusion between Two High-Risk Human Papillomavirus Types during a Coinfection[END_REF]. 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.

In conclusion, our model formalizing both the heterogeneity of sexual behaviors at the population level and the co-circulation of distinct groups of genotypes allows us to evaluate pathogen transmission dynamics in realistic sexual networks. In the context of HPV transmission, our simulation results show that interactions between HPV genotypes can significantly impact effectiveness of anti-HPV vaccination at reducing HPV prevalence. This model can be used to further investigate the impact of most recent anti-HPV vaccines targeting a larger number of genotypes, again altering prevalence equilibrium, as well as the recent extension of vaccination to males (starting in 2021 in France). Vaccination of boys is expected to increase the magnitude and speed of change in the prevalence of VT and NVT. Finally, the model could easily be adapted to study other research questions related to sexually transmitted pathogens and their interactions.

Fig. 3 .

 3 Fig. 3. Temporal evolution of the prevalence of NVT and VT genotypes according to interaction scenarios 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.

  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.

Table 1

 1

	and described below. Each agent is explicitly modeled and characterized
	by his/her sex, age, sexual activity class (number of partners he or she will have during the
	year), partnership status (in three categories -virgin, with a partner, or single), vaccination

Table 1 .

 1 Individual variables characterizing each agent

	Variable	Type	Update time
	ID	Code number	At entrance
	Sex	Boolean	At entrance
	Age	Number weeks	At each time step
	Vaccination status	Boolean	At entrance
		Variables related to partnership process
	Sexual activity class	Integer	At 15, 17, 19 and 24 years
	No partner during the current year	Boolean	On the first day of each year
	Partnership status	String	At the beginning and at the end of partnership
	Partner ID	Pointer	At the beginning and at the end of partnership
	End date of partnership or inac-tivity	Simulation week number	At the beginning and at the end of partnership
	Counter of partners during the current year	Integer	At the onset of new partnership; and at reset on first day of each year
	Counter of partners over the sim-ulation	Integer	At the onset of new partnership

Variables related to the infection process -for each genotype g

	Infection status	String	At the beginning and at the end of infection and im-munity
	End date for infection	Simulation week number	At the onset of infection
	Natural immunity status	String	At the end of infection
	End date for natural immunity	Simulation week number	At the end of infection
	Counter of infections over the simulation	Integer	At the onset of infection

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;

Table 2 .

 2 Environmental variables characterizing the population and their values

	Description	Default value	Rationale
	Population size (N)	800,000 individuals	For computational purposes
	Population age range	From 15 to 30 years of age	Fixed, corresponding to HPV in-
			fection peak
		Variables related to partnership process
	Age at first partnership	Normal distribution with mean 17.5 years in fe-	Mean, calibrated on (Bajos and
		males, 16.8 in male	Bozon 2008)
	Distribution of individuals by sexual	According to sex and age (Table 3 in appendix)	From (Bajos and Bozon 2008)
	activity class		and calibrated
	Proportion of individuals without a	According to sex and age (Table 3 in appendix)	From (Bajos and Bozon 2008)
	partner during the year		and fixed
	Duration of partnership	According to sexual activity class and age (Table 3	Mean, from (Bajos and Bozon
		in appendix)	2008) and calibrated
	Duration of inactivity between two	According to sexual activity class (Table 3 in ap-	Mean, calibrated
	partnerships	pendix)
	Frequency of sexual intercourse	2 per week	Arbitrary
	within a relationship
	Duration of partner search before	According to sexual activity class (Table 3 in ap-	Calibrated
	mixing between sexual activity	pendix)
	groups
	Maximum proportion of individuals	95% for females and 45% for males	Calibrated
	changing from one extreme class to
	another at ages 17, 19 and 24

Variables related to infection process

  Scales. The model time step is the week. A simulation year lasts 52 weeks. Agents are not explicitly distributed in space.

	Number of genotypes	2 vaccine types (VT)	According to genotypes re-
		12 non-vaccine types (NVT)	ported in (Markowitz et al.
			2013)
	Transmission probability parameter	βV = 0.16 per sexual act (VT)	Calibrated
	for VT and NVT genotypes	βNV = 0.125 per sexual act (NVT)
	Duration of infection	Exponential distribution with mean 52 weeks	From (Trottier et al. 2006)
	Duration of immunity	Exponential distribution with mean 12 weeks	Calibrated
	Genotype interaction parameter	1	To be varied
	(𝛾)
	Initial probability of infection	0.2 for NVT and 0.08 for VT	Arbitrary
	Initial probability of co-infection	0.4 for any first genotype g	Arbitrary
	among those initially infected
	Date of vaccine introduction	After 70 years of simulation	To ensure prevalence equilib-
			rium is reached
	Vaccine coverage	60% of females by age cohort	To be varied

Table 2 .

 2 Then all agents are created, with their partnership and infection statuses and counters initialized. At initialization, sex and age of each agent are randomly drawn to ensure homogeneous distribution of population across sex and age. Sexual activity class and age at first partnership are drawn from probabilistic distributions (if the agent is older than his/her age at first partnership, he/she can directly look for a partner).

  to define relevant parameters characterizing the partnership process and calibrated the model to reproduce realistic contact patterns in the country. This population-based cross-sectional survey provides information on sexual behaviors of a representative sample of 12 364 French males and females. We notably used data on the distribution of individuals in sexual activity classes by age categories (18-19, 20-24 and 25-29 years old), the distribution of individuals without partner during the year by

	2013), a country in which HPV epidemiology has been suggested to be similar to that in France
	(de Sanjosé et al. 2007; Bruni et al. 2010). From (Markowitz et al. 2013) we extracted HPV
	prevalence (i.e., the proportion of individuals infected in the population) of VT (gathering 2
	genotypes: HPV 16 and 18) and NVT (gathering 12 genotypes:

age categories, age at first partnership, age difference between partners, duration of partnership, and duration of inactivity between two partnerships.

Infection data. Because no data in the general population were available on the proportions of individuals infected with HPV by age and genotype in the pre-vaccine era in France, we used distributions reported from an epidemiological study carried out in the US

(Markowitz et al.

Acknowledgements

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

Funding 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).

Availability of data and materials

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

Abbreviations

ABM(s), Agent-based model(s); HPV, human papillomavirus; ODD, Overview, Design concepts, and Details; VT, vaccine type; NVT, non-vaccine type.

Competing interests

LO reports research grants from Pfizer through her institution, unrelated to the submitted work. All other authors report no competing interests. Not relevant a Individuals without a partner are randomly selected each year from the population; b When not specified fixed or calibrated, values are extracted from (Bajos and Bozon 2008); c One category of partnership duration is drawn from a multinomial distribution whose parameters are the proportions reported in [START_REF] Bajos | Sexualité, genre et santé : les apports de l'enquête « Contexte de la sexualité en France[END_REF], then a precise value is drawn from a uniform distribution within the duration category; d The duration of partnership and between two partnerships are identical in most sexually active individuals and drawn only once a year .

Authors' contributions

MB

Appendix


Marguerite Léang

Frédéric Ott

Frédérique Giorgiutti-Dauphiné

Ludovic Pauchard

Lay-Theng Lee

Imbibition and structure of silica nanoporous media characterized by neutron imaging

teaching and research institutions in France or abroad, or from public or private research centers.


Abdallah Nassereddine

Qing Wang

David Loffreda
email: david.loffreda@ens-lyon.frjaysen.nelayah@u-paris.frhazar.guesmi@enscm.fr

Christian Ricolleau

Damien Alloyeau

Catherine Louis

Laurent Delannoy

Jaysen Nelayah

Hazar Guesmi

Revealing Size Dependent Structural Transitions in Supported Gold Nanoparticles in Hydrogen at Atmospheric Pressure

Keywords: size dependent restructuration, gold nanoparticles, hydrogen, ETEM, AIMD

come   L'archive ouverte pluridisciplinaire

Introduction

Many scenarios have been proposed to explain the origin of the catalytic performances of gold nanoparticles (Au NPs) in catalytic and electro-catalytic reactions. [1,2,3,4] In particular, the higher activity of Au NPs with size smaller than 5 nm has been often attributed to sizedependent physical features such as the fraction of surface atoms. [5,6,[START_REF] Astruc | Nanoparticles and Catalysis Ch[END_REF] For instance, in reactions of hydrogenation, it is widely accepted that surface atoms on small particles are much more active than those on large particles since the fraction of low coordinated sites, i.e., at edges and corners, which are the hydrogen activation sites, increases when size decreases. [START_REF] Gatin | [END_REF]9,10,11] This concept has been largely adopted for studying the reactivity of Au NPs, especially by theoreticians and surface science researchers. [12] The common practice is to consider Au nanoparticles as perfect model crystals with well-defined symmetry and static surface facets interacting with low coverage of gas molecules. However, the effect of adsorbates, which induces surface changes is generally neglected in these studies, [13,14,15] which is a serious drawback that may prevent reliable description of the catalyst reactivity that mainly depends on the configuration of the surface. Aside from structural effects, quantum size effects are also frequently invoked to explain the exceptional activity of gold nanoparticles. [16,17] For instance, using Density Functional Theory, Illas et al. have reported that the minimum energy pathway of O 2 dissociation for different Au NPs is chemically similar which indicates that the larger activity of smaller ones originates from electronic and not from geometric effects only. [17] Furthermore, several other works strongly suggest that Au particles, below ∼2.5 nm in size, have non-metallic character which make them especially active in H 2 chemisorption and dissociation favoring H-rich environment promoting hydrogenation reactions. [START_REF] Bond | Catalysis by Gold[END_REF] The concept of fluxionality was also reported to be one of the keys to the understanding of the peculiar reactivity of small gold clusters. [START_REF] Arenz | [END_REF]20,21] Nevertheless, all these concepts and interpretations lack a tangible proof, as mechanistic understanding of the working catalyst requires access to the surface atomic evolution under relevant conditions of temperature and pressure.

The understanding of catalytic properties of Au NPs requires access with atomic resolution to the evolution of their structure under relevant conditions of temperature and pressure. Until very recently, this possibility was difficult to reach explaining the lack of a wide consensus on the mechanistic origin of the catalytic activity of Au NPs. Twenty years ago, using environmental transmission electron microscopy (ETEM), Hansen et al. have reported a landmark work evidencing shape changes in supported copper nanocrystals when exposed to low gas pressure under oxidizing and reducing conditions. [22] This technique has since emerged as a powerful tool for real time visualization of structural evolution of materials under dynamic gas reaction conditions. With the advent of high-pressure gas cells and TEM aberration correction, ETEM nowadays enables structural dynamics to be followed in gas conditions up to atmospheric pressure with resolution down to the atomic scale. [23,24] Additionally, ab initio molecular dynamics (AIMD) has been shown to be an efficient method for exploring the conformational three-dimensional complexity of metal NPs in presence of gas adsorbates, [25] providing dynamic trajectories that links different possible structural changes. In this work, the structural dynamics of Au NPs with sizes ranging from 7 nm to 1.5 nm supported on anatase-type titanium oxide (Figure S1, Supporting Information) was examined for the first time, by ETEM at the atomic scale upon exposure to H 2 at atmospheric pressure and at various temperatures. In situ atomic scale observations of individual Au NPs complemented with ex situ TEM observations, show supported fcc Au NPs undergoing dramatic structural changes in H 2 when their size is smaller than  4 nm. AIMD simulations disclose the atomic behavior driving the structural changes under H 2 observed in such NPs and provide a precise picture of the hydrogenated gold surface at the molecular level.

Results and Discussion

Figure 1 shows series of in situ high angle annular dark field scanning ETEM (HAADF STEM) side-view images of three Au NPs with size ranging from 4.5 nm to 3 nm under 10 5 Pa H 2 . Figure 1a and 1b compare the structural evolution of a 4.5 nm (a) and a 3 nm (b) truncated octahedron (TOh) Au NP when they are cooled from 400 to 25 °C in H 2 . The full series of in situ HAADF and bright-field (BF) STEM images acquired on these two NPs are provided in (Figure S2 and S3 respectively, Supporting Information) and all experimental details are described in the Methods section. Under H 2 exposure, size-dependent structural changes of Au NPs are revealed. As reported in Figure 1a and (Figure S2, Supporting Information), the gold NP of 4.5 nm is found to roughly retain its initial shape and its fcc structure down to 25 °C. Similar results are found for a larger NP presented in (Figure S4, Supporting Information). In contrast, both the shape and the structure symmetry of the Au NP with 3 nm size drastically evolve upon exposure to H 2 at 400°C (Figure 1b and Figure S3, Supporting Information). Indeed, comparison of STEM image simulations with experimental in situ data shows that the initial fcc structure transform into icosahedron (Figure S5, Supporting Information) with the five-fold icosahedral symmetry remaining stable as temperature is reduced to 200 °C. Additional in situ ETEM observations of another 3-nm Au NP in 10 5 Pa H 2 at 175 °C show instead a non-fcc structure with unidentified symmetry evolving in the presence of hydrogen (Figure 1c). As illustrated by the HAADF STEM images extracted from a time-lapse video recording (Video S1, Supporting Information), the structure of this Au NP fluctuates continuously between different structural isomers with shape in projection remaining roughly spherical during the course of the video recording (144 seconds duration). S2c and S3c, Supporting Information). In H 2 , the 4.5-nm NP remains stable down to 25 °C while the 3-nm NP restructures at 400 °C, leading to changes in both its structure and shape with the observation of an icosahedral structure under H 2 as from 400 °C. c) Sequential in situ STEM images of another 3 nm Au NP in 10 5 Pa H 2 at 175 °C, showing a round-shape particle undergoing continuous restructuring under H 2 (captured from Video S1, Supporting Information). In order to further investigate the size-dependent structural evolution of gold nanoparticles when exposed to H 2 , we perform complementary ex situ STEM imaging under vacuum and at room temperature to analyze the local structures of Au NPs previously exposed to H 2 in a fixed-bed reactor under pressure and heating conditions close to those of in situ TEM observations (see Experimental section). Figure 2 shows atomic scale ex situ HAADF STEM images of typical NPs exposed ex situ to H 2 . The analysis of the structure of 61 NPs clearly identifies a critical particle size around 3-4 nm. While particles with size greater than 4 nm display perfect fcc structure bounded by well-defined facets (Figure 2a-c), ex situ TEM identifies various non-fcc structures for sizes smaller than 4 nm. We observe particles with a core reminiscent of an fcc order (core showing lattice deformation) and a highly-disordered surface (Figure 2e-f), NPs showing five-fold symmetries (Figure 2g) and disordered structure NPs (Figure 2h-l). Ex situ TEM observations thus provide evidences of the important effect of adsorbed hydrogen on the structural transition of gold NPs being size-dependent. Based on these both ex situ and in situ TEM, we can firmly state that NPs below 4-3 nm undergo substantial structure evolution when exposed to H 2 , resulting in the loss of their fcc arrangement. In contrast, H 2 exposure has roughly no effect on the structure of Au NPs with size above 4 nm, which is in line with theoretical results predicting that low index gold surfaces as those forming the surface of TOh NPs are unable to instantly adsorb H 2 . [15,26,27] These experimental observations provide the first evidence of particle size dependent structural transformation of gold NP under hydrogen. Moreover, while hydrogen adsorption on Au NPs has already been demonstrated by DFT calculations and spectroscopic measurements, [28,29] the adsorption inducing structural transformation had never been directly observed.

In order to understand the atomic origin of the structural changes observed experimentally in small sized Au NPs under hydrogen pressure, we performed AIMD simulations of model TOh hydrogen-covered Au NPs at 500 K (227°C) and at 300 K (27°C). We follow the structural evolution dynamics of Au NPs of 201-atom (2 nm) in vacuum (denoted as TOh_Au 201 ) and with a coverage of 1ML atomic hydrogen (denoted as TOh_Au 201 H 122 and corresponding to the adsorption of 122 H atoms, see Methods). The AIMD simulations show a strong evolution of the structure of the Au NP induced by the adsorption of atomic hydrogen on its surface (Video S2, Supporting Information). Within few ps of dynamics performed at 500K, the TOh_Au 201 H 122 NP (Figure 3a, configuration A) loses its TOh equilibrium shape and transforms into a rounded shape NP with a highly distorted non-fcc structure (Figure 3a In addition to the evidenced geometric ensemble effect where all gold surface atoms are henceforth lower coordinated sites, the electronic structure of the gold NP surface sites is unique. This is primary illustrated by the distribution of charge density at the NP surface depicted in (Figure 4a) showing a strong accumulation of electron density over the Au-H-Au crown lines and a strong depletion in between (The inner shell atoms). The Bader charge calculations show that the Au atoms of the crown lines have larger negative charges (>0.04e, with e the elementary charge) than those from the inner shell which are almost positive (see enlarged part in Figure 4a). This suggests that the surface restructuration may offer two types of accessible Au sites for reactants and intermediates of hydrogenation catalytic reactions.

Finally, the analysis of the structure of the NP core shows important atomic ordering changes with surface segregation of bulk Au atoms resulting in a strongly disordered inner shell and a highly distorted icosahedral-like center (Figure 3d). The latter recalls the icosahedral structures observed experimentally (Figure 1b). These results elucidate how hydrogen adsorption induces fcc to icosahedral-like structural transformations in Au NPs observed by in situ TEM. They are in excellent agreement with the reported classical MD results of H.

Ogawa predicting the possible transformation from fcc to icosahedral arrangements for a given hydrogenated fcc metal nanoparticle smaller than 2 nm. [30] This author also reported that this strong disorder and lattice deformations might occur if the M-H interactions are strong enough. In our work, the adsorption of atomic hydrogen is found to stabilize at 0.8 ML coverage with adsorption energy of -0.025 eV/H. This value is higher than the calculated adsorption energies on flat gold surfaces or perfect crystal NPs, [31,32,33,34] but is fully consistent with experimental experimental works predicting the dissociative chemisorption of H 2 on gold catalyst [10,28,35,36,37] for temperature in the range of 298 -523 K [10] providing by the way, a rational explanation of the observed structural transformations. The predicted strong hydrogen adsorption governing the structural changes is also illustrated by the drastic changes in the Au 5d-band density of states of Au atoms from both the surface and the center in configuration D' compared to bare TOh_Au 201 NP. Figure 4b and4c show that, while the occupied orbitals of the hydrogenated surface gold atoms downshift away from the Fermi level as expected upon hydrogen chemisorption, those of the center atoms upshift toward lower energies, signature of charge transfers for the benefit of the upper surface atoms. In addition, the narrowing of the d-band density of state of the center Au atoms (19 atoms) found for configuration D' compared to TOh_Au 201 is in complete agreement with the theoretically predicted and experimentally observed transition from fcc to icosahedral-like symmetry. [38] All these results open the way for a deep investigations of the reactivity of these revealed surface active sites to understand the "real" reaction mechanisms occurring over hydrogenated gold nanocatalyst.

Conclusion

In conclusion, using environmental transmission electron microscopy, we provide the first real space evidence of the size dependent structural evolution on Au NPs under hydrogen. We reveal symmetry transformations from fcc to icosahedral-like or other undefined structures occurring in Au NP with size below 4 nm. Performing AIMD simulations, we demonstrate that these structural evolutions resulting from the high mobility of gold atoms are induced by strong Au-H interactions and we bring insights into the unique atomic and electronic structure of the interfacial sites of Au NPs and hydrogen in this critical size regime. This work recalls the necessity of considering the adsorption-inducing restructuration to explain the sizedependent catalytic activity of gold nanocatalysts. It shows that the idea of low coordinated surface atoms at particle edges and corners as the main active sites is incomplete and that in addition, ensemble and electronic effects induced by the surface and core restructuration have to be considered. This work also opens the question of the mechanistic origins of the structural flexibility of Au NPs under reactive gas as a key point for the understanding and the optimization of active gold-based heterogeneous nanocatalysts.

Experimental Section

Sample preparation: A titania-supported gold catalysts (Au/TiO 2 with 6 wt % Au) sample was prepared by deposition-precipitation technique in the presence of urea. This preparation method produces gold NPs with average size around 3-4 nm and catalytically active in selective hydrogenation reactions. [39,40] The range of Au particle sizes studied in this work is typically the one encountered in the literature for active gold catalysts in various oxidation or hydrogenation catalytic reactions. [START_REF] Bond | Catalysis by Gold[END_REF] Moreover, experimental limitations resulting from electron scattering by the silicon nitride windows of the high-pressure gas cell, and STEM contrast from the TiO 2 support, it was not possible to resolve the atomic structure of Au NPs smaller than ~3 nm. As support, a commercial TiO 2 (100 % anatase, US research nanomaterials, surface area 50 m 2 .g -1 , purity = 99.98%) was used. First, a suspension of TiO 2

(2 g) in distilled water (150 ml) was prepared in a double wall-reactor. The solution was stirred and heated in order to well disperse the support up to 80 °C. Then, an appropriate quantity (27.5 mL) of a solution of gold precursor HAuCl 4 .3H 2 O (Sigma-Aldrich, 4.5 g.L -1 ) to achieve Au loading of 6 wt%, was added to the TiO 2 suspension. Instantly, urea CO(NH 2 ) 2 (Sigma-Aldrich) was also added in excess ([urea]/ [metal] = 100). The mixture at 80 °C was then stirred 20 hours in the dark to avoid uncontrolled reduction of gold. The obtained solid was then separated from the liquid phase by centrifugation. After centrifugation, the solid was washed three times with distilled water to remove impurities, in particular residual chlorides.

The sample was then dried under vacuum, in dark and at room temperature for 24 hours. The Au gold loading determined by X-ray fluorescence was 5.8 wt %. The catalyst was then activated in a fixed bed microreactor in Pyrex (home-built, 4 mm of internal diameter), leading to the so-called as-prepared Au/TiO 2 ., i.e. gold was reduced into metal particles, by temperature-programmed reduction (ramp 2 °C.min -1 ) under flow of H 2 (Air Liquide, purity 99.999 %), (100 ml.min -1 ) from RT to 300 °C then maintained at 300 °C for 2 hours of pure hydrogen. The microreactor was then purged with N 2 (Air Liquide, purity 99.999 %) before exposure to air. The same sample batch of as-prepared Au/TiO 2 was used for all the experiments presented in this work.

In situ environmental TEM: All electron microscopy observations were carried out in double aberration-corrected JEOL ARM 200F microscope equipped with a cold field emission gun operated at 200 kV. For in situ STEM observations in gas phase, a Protochips Atmosphere TM TEM environmental high-pressure gas cell (HPGC) was used. It consists of a pair of silicon microchips (called E-chips) aligned and assembled in closed-cell configuration at the tip of a dedicated TEM sample holder. The small E-chip is a 2 mm × 2 mm, 300 μm-thick silicon chip with a 300 μm × 300 μm 50 nm-thick amorphous SiN window at its center. As for the large Echip, it is a 6.0 mm × 4.5 mm and 300 μm-thick Si microchip with a 300 μm × 300 μm SiC membrane at its center. The SiC film is patterned with six (2 lines of 3) 7 μm in diameter circular windows made up of 30 nm thick amorphous SiN. The as-prepared Au/TiO 2 sample were deposited on the large E-chip by drop casting 6 µl of a suspension of Au NPs ultrasonically dispersed in distilled water. TEM observations were carried out through these electron-transparent windows. The Au NPs were subjected to pure argon (Messer, purity > 99.999 %) and pure hydrogen (Messer, purity > 99.999 %) in static gas conditions using a manually-operated gas manifold from Protochips. Gas pressure in the HPGC cell was controlled using a MKS PDR 2000 dual capacitance manometer. The nominal gas layer thickness was 5 µm fixed by gold spacers patterned on the E-chips HPGC. The sample was heated through resistive heating of the SiC membrane via current injection using an electronic control unit interfaced with Protochips Aduro 300 software with closed-loop temperature control. For the latter, the software uses chip-specific ohmic resistance of the SiC film versus temperature calibration files as input. Gold NPs were imaged in STEM mode, using both the HAADF and BF detectors. The probe setting used for image acquisition was 8C and the convergence semi-angle was set to 20.5 mrad (by using a 40 μm condenser lens). The resulting probe size was 0.7-0.8 Å and the corresponding probe current was close to 27 pA (both quantities being measured in vacuum). The camera lens was set to 8 cm corresponding respectively to inner and outer collection semi-angles of 68 and 280 mrad for HAADF imaging. For BF imaging, the collection semi-angle was limited to 17 mrad by using a 3 mm BF aperture. High-resolution STEM images (2048 × 2048 pixels) were acquired at magnifications of × 10 M and × 20 M (figure 1b and c only) with pixel area equal to 8.97×10 - 3 and 2.24×10 -3 Å 2 respectively. To limit the effect of the electron beam during in situ STEM imaging, all images were captured with a limited acquisition time of 2 or 3 µs per pixel and the electron illumination was blanked between image acquisition. Under these acquisition conditions, the dose per frame, as defined in [41] was equal to 3.76×10 10 or 5.64×10 10 e Å -2 at × conditions. Another important simulation parameter is the exponent power (also denoted as the Exponent number in the software) which describes the dependency of the HAADF intensity with the atomic number. Its value varies with the inner acceptance angle of the HAADF detector as shown in. [43] Based on the selected collection conditions, its value was fixed to 1.7 in all simulations. Atomic structure models used as inputs of STEM simulations were either taken from the simulation software database (bare 931-atom icosahedron in Figure S5, Supporting information) or extracted from AIMD simulations of the 201-atom Au NP with 1 ML coverage of hydrogen (simulated atomic structures in Figure 3).

Cluster model: Regular truncated octahedron of 1.8 nm size composed by 201 Au atoms (TOh_Au 201 ) was considered. [44] This size is representative of gold nanoparticles undergoing structural transformations under hydrogen in our TEM observations. To account for the presence of atmospheric hydrogen pressure, the hydrogen coverage was fixed to 1ML, which corresponds to the adsorption of 122 H atoms (TOh-Au201-H122). The gold NP model structures were modeled using a cubic cell of dimensions The NP models for MD simulations (TOh_Au 201 and TOh_Au 201 H 122 noted configuration A in Figure 3a)

were initially optimized by DFT calculations, leading to fcc local minimum energy structures with slightly distorted surface facets, for the case of hydrogen adsorption with H atoms mainly occupying twofold bridge and threefold hollow sites.

AIMD simulations and DFT calculations: Ab Initio Molecular Dynamics simulations

involving Born-Oppenheimer approximation were performed in the canonical (NVT) ensemble employing Nosé-Hoover thermostat with a time step of 1.5 fs at 300K and 500K with a total integration time of 60 ps. [45,46] All ab initio MD simulations and DFT calculations were conducted using Vienna ab initio Simulation Package (VASP 541). [47,48] All local minima associated with specific atomic arrangements selected from MD simulations were isolated and their geometrical and electronic properties were calculated and analyzed by DFT method. The Generalized Gradient Approximation functional of the exchange-correlation energy was calculated within the Perdew, Burke, and Ernzerhof formulation of the generalized-gradient approximation (GGA-PBE). [49,50] Van der Waals interactions have been considered for all the models throughout the DFT study by selecting the dispersion corrected PBE-D3 functional proposed by Grimme. [51 ] The cut off energy was fixed to 400 eV. Fine geometry optimizations were performed by choosing a convergence threshold of 10 -6 eV for the total electronic energy and a threshold of ±10 -2 eV•Å -1 for minimizing residual forces acting on the nuclei. In order to describe the interaction between gold nanoparticle and atomic hydrogen, the adsorption energy normalized by the number of H atoms (N_Hads) was defined as follows:

ΔEads = (E(TOh+adsorbate) -E(TOh) -E(adsorbate))/N_H ads , where E(TOh+adsorbate), E(slab), E(adsorbate) are the calculated electronic energy of H species adsorbed on the gold nanoparticle, of bare TOh nanoparticle, of gas phase H 2 , respectively. The adsorption energy per H atom referres to 1/2H 2 (g).

Table of Contents

Table of Contents Figure

Figure 1 .

 1 Figure 1. Direct atomic scale visualization of the size-dependent structural evolution of Au NPs supported on anatase-type titania in H 2 at atmospheric pressure. a-c) Time-and temperature resolved in situ HAADF STEM images reveal distinct structure transformation pathways for NPs greater than 4 nm (a) and smaller than 4 nm (b,c). Prior to H 2 exposure and as a reference experiment, the two NPs were individually imaged in 40 Pa argon (near-vacuum condition) at 400 °C. In situ HAADF STEM imaging of a  4.5-nm Au NP (a) and a  3-nm Au NP (b) are sequentially acquired in 40 Pa Ar at 400 °C then in 10 5 Pa H 2 at 400, 200 and 25 °C. The broken white lines highlight the projected truncated octahedral equilibrium shapes of both NPs in Ar at 400 °C. The

Figure 2 .

 2 Figure 2. Atomic scale STEM imaging of TiO 2 supported Au NPs after ex situ exposure to H 2 at atmospheric pressure at 400 and 200 °C. Ex situ exposure to hydrogen is conducted on the same starting materials used for in situ experiments in (Figure 1). Gas pressure and temperature conditions for ex situ experiment are consistent with in situ ones (see Experimental Section). TEM pressure during STEM image acquisition,  10 -6 Pa. Particle size decreases from (a) to (l). The scale, 1 nm, in (e) is common to all subsequent Figures. a, c) 4.3-nm and 4.1nm particles with fcc structure (FFT in (b) and (d) respectively). e-l) For particle size below 4 nm, highly distorted non fcc structures are observed.

  , configuration B), which is primarily caused by surface motion and extraction of Au atoms induced by the adsorption of hydrogen. The fcc to non-fcc structure transition is evidenced by examining the distribution of Au atom distances with respect to the mass center of the NP. As depicted in Figure 3b, the well-defined peaks of TOh_Au 201 signature of fcc structure progressively disappear during adsorption of hydrogen, suggesting the formation of a disordered structure. During the first 20 ps of time equilibration (Figure 3a, configuration C), combination of atomic hydrogen and desorption of H 2 molecules occur at the surface of the Au NP and the number of adsorbed atomic hydrogen reduces to 102, i.e., to 0.8 ML. The configuration C shows the initial stage of formation of several H-Au-H-Au chains on the surface of the Au NP, with one hydrogen bridging two gold atoms. After 60 ps time span, well-organized long H-Au-H-Au crown lines are stabilized (Figure 3a, configuration D). They are found to continuously slide all around a highly distorted gold core. Interestingly, HAADF-STEM simulation of the selected hydrogenated MD structures at 0, 20, 40, 60 ps (Figure 3c) gives rise to image contrast similar to the ones observed experimentally (Figure 2e-l). These MD simulation results thus highly support the experimentally-observed structural changes in Au NPs smaller than  4 nm. Further analysis of the dynamic behavior of individual atoms during the simulation time reveals strong structural fluctuations with Au-Au bond breaking and formation of new Au-Au and Au-H bonds (Figure S6, Supporting Information). This Figure shows the relative motion between atoms monitored by following the separation of Au-Au and Au-H atom pairs at the particle surface (i.e. at the crown chains) and of Au-Au atom pairs in its core (i.e. in the inner shell and in the center of 19 atoms) during time simulation. While highly mobile gold and hydrogen atoms at the surface reach the stage of bond retention after 40 ps corresponding to the formation of well-ordered H-Au-H-Au crown chains in configuration D, the Au-Au bond breaking and re-formation still proceed after 60 ps simulation time span for gold atoms in the particle center. These simulated structural fluctuations are in agreement with the strong structural instability observed in the gold nanoparticles presented in (Figure 1c and Video S1, Supporting Information) Similar results are obtained for the MD simulations performed at room temperature.

Figure 3 .

 3 Figure 3. Ab initio molecular-dynamics (AIMD) simulations of model gold NP in presence of hydrogen and DFT optimized structure. a) Optimized structure of TOh_Au 201 with fcc structure. Different configurations extracted from the AIMD simulation at 500K of TOh_Au 201 covered by adsorbed hydrogen atoms: Configuration A (t=0) obtained after DFT geometry optimization (at 0K) of TOh_Au 201 H 122 . Configuration B (t=3 ps) corresponding to the earlier stage of structural change process. Configuration C (t=20 ps) obtained after thermal equilibration, with 116 adsorbed H atoms and 3 desorbed H 2 molecules in the gas phase. Configuration D (t=60 ps) with well-formed linear chains on the surface, 102 adsorbed H atoms and 10 desorbed H 2 molecules. Color

Figure 4 .

 4 Figure 4. Electronic structure analyses of the optimized DFT structure (D') of hydrogenated Au NP. a) 3D charge-density difference on configuration D' with an isosurface value of 0.00015 e/bohr 3 . Color code: yellow, gaining electrons; blue, losing electrons. The inset reports the Bader charges of a selected region where red and black values are for hydrogen and gold atoms, respectively. b,c) Au 5 d-band density of states of the surface crown lines (a) and of the center (19 atoms) (b) calculated from D' configuration (red curves) and compared with initial non hydrogenated gold NP, TOh_Au 201 (black curves). The d-band centers of surface atoms from hydrogenated and non-hydrogenated Au NP are calculated to be of -3.33 eV and -3.20 eV, respectively. The dband centers of center atoms (19 atoms) from hydrogenated (icosahedral-like) and non-hydrogenated (fcc symmetry) Au NP are calculated to be of -3.47 eV and -3.72 eV, respectively. E f is the Fermi energy.

  Text Atomic-resolution transmission electron microscopy reveals size-dependent structural transformations in supported gold nanoparticles upon exposure hydrogen at atmospheric pressure. Molecular dynamics simulations calculations highlight the role of mobile gold atoms in the evolution of particle crystal structure and morphology observed below 4 nm and unveil the unique atomic and electronic structures of hydrogenated gold nanoparticles in this critical size regime.

Acknowledgements

This work was supported by the French National Research Agency through the TOTEM project, Grant Number ANR-17-CE07-0031. A. N was supported by a PhD fellowship granted through the same project. The authors are grateful to Region Ile de France for convention SESAME E1845 for the support of the JEOL electron microscope installed at Université de Paris. The authors thank IDRIS in Paris, CINES in Montpellier, TGCC in Grenoble (project 609 and A0090807369, GENCI/CT8) and especially PSMN in Lyon for CPU time and assistance. The authors thank the CPER/SYSPROD 2015-2022 project (N°2019-AURA-P5B) and AXELERA Pôle de Compétitivité (PSMN Data Center). The authors are also grateful to GDR Or-Nano (GDR n°2002) for his work to stimulate effective collaborations. A. N and Q. W contributed equally to this work. Received: XXXX Revised: XXXX Published online: XXXX

Data Availability Statement

The ETEM data and STEM simulations that support the main findings of this study are properties of the Quantum Materials and Phenomena laboratory of Université de Paris. They are available from the corresponding author J. N on request. All DFT and MD simulation results are also available from the corresponding authors H.G and D. L. on request.

10M magnification and was raised to 1.51×10 11 or 2.26×10 11 e Å -2 at × 20M magnification.

The same doses per frame were selected for ex situ STEM observation.

Ex situ STEM: Ex situ observations of the as-prepared Au/TiO 2 catalyst was carried out under static atmospheric pressure of pure hydrogen (Air Liquide, purity 99.999 %) in the home-built fixed bed microreactor in Pyrex. The catalyst (15 mg) was heated from room temperature to 400 °C (ramp 10 °C.min -1 ) under H 2 , held for 1 h at this temperature, then cooled to 200 °C (ramp 10 °C.min -1 ) and maintained for 1 h at 200 °C before cooling to room temperature. After heating, the microreactor was purged with N 2 (Air Liquide, purity 99.999%) before exposure to air. The pressure and temperature conditions under hydrogen were chosen to be consistent with those applied during in situ TEM observations. The H 2treated Au/TiO 2 catalyst was then recovered and transferred in air for high resolution ex situ STEM observations at room temperature and under vacuum ( 10 -6 Pa). The sample for electron microscopy was prepared by dispersing the hydrogen-exposed catalyst in distilled water and depositing 10 drops of 2 µl of the resulting solution onto a holey carbon-coated copper TEM grid (Agar Scientific). For STEM imaging, the illumination and image acquisition conditions were similar to those set for in situ STEM.

STEM simulation: Atomically-resolved HAADF STEM image simulation was performed with the STEM simulation software developed by He, Li and Yuan using Matlab 7.12.0.635 (R2011a). [42] This software assumes kinetic scattering of the electrons (i.e. the incident electrons are only scattered once) by the atoms of the structure under consideration. Within this approximation, which holds for atomic columns with less than 10-50 atoms, the integrated HAADF intensity at each point of scan is proportional to the total number of atoms encountered by the incident electron along its propagation path. The Gaussian width that describes the effective width of the electron probe was set to 0.7 Å in all simulations to match the spatial resolution of the JEOL ARM 200F microscope under the selected illumination

Supporting Information

Supporting Information is available from the Wiley Online Library or from the author.

Supporting Information

Supporting Information is available from the Wiley Online Library or from the author.

Conflict of Interest

The authors declare no conflict of interest.

Author Informations

These authors contribute equally: Abdallah Nassereddine, Qing Wang