Stéphane Caron
email: stephane.caron@ens.fr

George Kesidis
email: kesidis@engr.psu.edu

Incentive-based Energy Consumption Scheduling Algorithms for the Smart Grid

In this paper, we study Demand Response (DR) problematics for different levels of information sharing in a smart grid. We propose a dynamic pricing scheme incentivizing consumers to achieve an aggregate load profile suitable for utilities, and study how close they can get to an ideal flat profile depending on how much information they share. When customers can share all their load profiles, we provide a distributed algorithm, set up as a cooperative game between consumers, which significantly reduces the total cost and peak-to-average ratio (PAR) of the system. In the absence of full information sharing (for reasons of privacy), when users have only access to the instantaneous total load on the grid, we provide distributed stochastic strategies that successfully exploit this information to improve the overall load profile. Simulation results confirm that these solutions efficiently benefit from information sharing within the grid and reduce both the total cost and PAR.

I. INTRODUCTION

The current U.S. electrical grid is built up according to a static, centralized structure: remote power plants transmit electrical power through long-distance high-voltage lines to substations (transmission network), which in turn adapt and deliver it to local end users (distribution network). In this model, the local network is often statically tuned to match a given average load profile from its consumers.

Yet this structure is about to undergo a major shift: the progressive integration of smart meters [START_REF] Krishnan | Meters of Tomorrow[END_REF] and communicating appliances will upgrade this "blind" system to a decentralized "smart grid" [START_REF]The Smart Grid: An Introduction[END_REF], which is foreseen as a way to save billions dollars in energy consumption [START_REF] Kannberg | GridWise: The Benefits of a Transformed Energy System[END_REF].

In a smart grid infrastructure, utilities can set up dynamic tariffs incentivizing customers to adjust their loads to the current state of the network. This key feature, known as demand side management, will yield several benefits, including:

• Integration of intermittent energy sources such as wind or solar power at the distribution level [START_REF] Ipakchi | Grid of the future[END_REF]; • Demand Response: customers will be encouraged to shift their heavy loads to off-peak hours; • Resilience to attacks or power outages, with users spurred to turn off their "non-critical" devices in case of heavy load on the grid or upstream outages [START_REF] Pipattanasomporn | Multi-Agent Systems in a Distributed Smart Grid: Design and Implementation[END_REF]; and • Energy savings: studies [START_REF] Darby | The Effectiveness of Feedback on Energy Consumption[END_REF] already suggest that giving customers access to real-time consumption information yields significant savings.

Such load management is becoming even more crucial as plug-in hybrid electric vehicles (PHEVs) are coming to the market. With battery capacities varying from 15 to 50 kWh, these vehicles are expected to double the average household load during charging time [START_REF] Ipakchi | Grid of the future[END_REF].

Therefore, the design of appropriate incentives and efficient energy consumption scheduling (ECS) algorithms is a main issue for the deployment of the upcoming smart grid.

In [START_REF] Mohsenian-Rad | Optimal and Autonomous Incentive-based Energy Consumption Scheduling Algorithm for Smart Grid[END_REF], the authors deal with ECS in the case of increasing strictly convex cost functions. They propose a distributed algorithm and show through a game-theoretic analysis that, for incentives satisfying certain properties, it yields optimal energy consumptions for end users. However, they implicitly assume that the daily load on the network is proportional to the daily cost for the utility, with a constant independent of load scheduling (i.e., of game dynamics): this is a strong hypothesis which implies utilities' costs are linearly bounded in any situation.

In this paper, we survey different scheduling problems depending on the DR architecture (centralized or distributed) and the degree of knowledge appliances have on the state of the network. As in [START_REF] Mohsenian-Rad | Optimal and Autonomous Incentive-based Energy Consumption Scheduling Algorithm for Smart Grid[END_REF], we embed consumers in a local distribution network consisting of a single energy source (e.g., a step-down substation) supplying several load subscribers. Customers are incentivized to move their loads to off-peak hours through marginal costs which are linearly increasing with instantaneous value of network overload.

The rest of this paper is organized as follows. We introduce our pricing scheme and notations in section II. In section III, we see that the general scheduling problem is NP-hard, and provide a distributed solution set-up as a cooperative game between consumers. We provide stochastic policies for a decentralized setting in section IV. In section V, we derive the best distributed policy for synchronized users in a power grid, which we use as a reference in section VI where we provide experimental results. We draw conclusions in section VII and suggest future work in section VIII, including considerations of interruptible and non-uniform demand profiles.

II. PROBLEM SET-UP

We consider a T -hour time period, e.g., T = 6 hours from midnight to 6AM, during which N customers need to automatically schedule their electrical jobs. Note that our time horizon is finite since users do not want their jobs to be delayed forever.

A. Loads and Costs

The n th customer has demand profile D n parametrized by flexible start time s n and fixed (d n , τ n ) parameters, where d n denote the instantaneous power consumption of the job and τ n its duration:

D n (t) = d n 1 {sn≤t≤sn+τn} where 0 ≤ s n ≤ T -τ n .
We assume that once the service start time s n is selected, it cannot be interrupted by the user. The total instantaneous load on the network is then

λ(t) := n D n (t) = n d n 1 {sn≤t≤sn+τn}
We denote by C(λ(t)) the cost, in $/kW, experienced by the utility at time t (for the costs charged to consumers, see equation ( 4)). It depends on the instantaneous load λ(t), and the function C itself is likely to depend on additional system parameters. For example, the two-step conservation rate model used by BC Hydro [START_REF] Bc Hydro | Electricity Rates[END_REF] (parametrized by load threshold L) is

C L (λ(t)) = C 0 • 1 {λ(t)<L} + C 1 • 1 {λ(t)≥L} .
In a more general setting, C can be any smooth convex function of λ(t). However, in this article, we will focus on a ramp cost function with load threshold L > 0:

C L (λ(t)) = C 0 + C (λ(t) -L) + (1)
where the base cost C 0 and the overage rate C are positive constants (x + denotes max(0, x)). Threshold L corresponds to the load upon which the utility experiences overages, and therefore raises the cost to dissuade customers from scheduling their jobs (i.e., the second term in sum (1) represents the marginal costs of producing electrical power by an overloaded plant). Otherwise, λ < L corresponds to the plant's nominal operational regime. Finally, we will call "Global Cost" the overall cost (in $) for the utility:

GC := T t=0 n D n (t)C L (λ(t))dt = T t=0 λ(t)C L (λ(t))dt
With a ramp pricing scheme, this global cost becomes:

GC ramp = GC 0 + C T 0 λ(t)(λ(t) -L) + dt. ( 2
)
where GC 0 := C 0 n d n τ n is a schedule-independent incompressible cost.

B. Non-triviality Criterion

We are interested in scenarios where there is too much demand for the system to avoid overages, and so it has to cope with such situations. A simple way to insure this is for jobs to meet the following criterion:

n d n τ n > LT. (3)

III. COMPLETE KNOWLEDGE SETTING

In this section, we survey scheduling when all the jobs' characteristics (d n , τ n ) are known, either to all players or to a single entity who tries to find an optimal schedule for the whole system. This is for example the purpose of the global controller in [START_REF] Molderink | A Three-Step Methodology to Improve Domestic Energy Efficiency[END_REF]. Web portals like Google PowerMeter [START_REF]Google PowerMeter[END_REF], OPOWER [11] or CustomerIQ [START_REF]Silver Sprint CustomerIQ web portal[END_REF] also centralize energy consumption data about their users, which they can use thereafter to derive an efficient schedule and advise consumers to conform to it 1 .

We will first remind that finding an optimal schedule is an NP-hard problem and discuss a greedy approach for it. We will then consider a distributed algorithm set-up as a game between consumers, and derive an optimal strategy for it.

A. NP-hardness

When load profiles (d n , τ n ) are different for different users, the problem of minimizing GC ramp over all start times {s n } N n=1 is NP-hard [START_REF] Ausiello | Complexity and Approximation: Combinatorial Optimization Problems and Their Approximability Properties[END_REF]. Even when all durations τ n are equal (or similarly all d n are equal), finding an optimal schedule is still an NP-hard problem (i.e., the BIN PACKING problem).

B. Discussion : Greedy approach

Since the overall problem is NP-hard, one can consider approximating its optimal solutions, e.g., using well-known metaheuristics such as simulated annealing. Though we won't investigate how these techniques would perform, we will give an incremental greedy solution that may get trapped in suboptimal local extrema.

We consider inserting jobs in a given order i 1 , . . . , i N and denote by λ k the load profile after jobs i 1 , . . . , i k have been scheduled (with λ 0 ≡ 0). Given λ k , we want to schedule the (k + 1)

th job so as to minimize the global cost incurred by

λ k+1 . Let I k+1 (s) := [s; s + τ k+1 ] ∩ {t|λ k (t) > L -d k+1 }, λ + k = (λ k -L) + and λ - k = (L-λ k ) + .
Then, one can show that minimal GC ramp for λ k+1 is achieved when s k+1 minimizes:

I k+1 (s k+1 ) (d k+1 -λ - k ) 2 + 2d k+1 λ + k + L(d k+1 -λ - k )

C. Distributed Action

Now we assume consumers have complete knowledge of each others demands and play a game where they seek to minimize the global cost GC ramp . We provide an effective strategy for players, derived in a pessimistic setting, which turns out to be efficient at peak shaving and yields very good results in practice (see section VI).

Here "complete" knowledge means players will either communicate their demand profiles or make inferences about others demands based on repeated observations (e.g., night after night of [0, T ] = [12 AM-6 AM] activity).

To incent customers to minimize GC ramp (which correspond to the actual cost of supplying their demand, or an upper bound of it), utilities may charge customer i with an amount b i proportional to both the energy he consumed and the global cost, e.g.,

b i := d i τ i j d j τ j × GC ramp = C 0 d i τ i × GC ramp GC 0 , (4)
where C 0 d i τ i is the minimal possible cost for scheduling player i's job.

In what follows, we denote by F i the Cumulative Distribution Function (CDF) of the start time of player i, F i (t) := P[s i ≤ t], and f i the density of dF i . We also define φ i (t) := F i (t) -F i (t -τ i ) which is the probability of job i being active at time t.

Considering equation ( 2), we can upper bound GC ramp as follows:

GC ramp ≤ GC 0 + T 0 λ(t) 2 dt =: GC bound .
Let us consider a game where users seek to minimize the expected value of GC bound , which is the same as minimizing E[ λ2 ]. Since λ 2 = (λ -µ) 2 plus a constant, where µ denotes the temporal mean of λ, this goal is closely related to peak shaving.

For any user i, let us define:

H i (t) := j =i d j (φ j (t) -φ j (t + τ i )) ,
which does not depend on F i . We have:

E T 0 λ(t) 2 dt = i,j d i d j T 0 φ i (t)φ j (t)dt = i d i T 0 F i (t)H i (t)dt =: i γ i .
The game between customers goes like this: users play asynchronously, and at his/her turn, player i updates F i in order to minimize

γ i ∝ F i H i . Claim 1. The optimal CDF F * i minimizing F i H i for any given (right-continuous) function H i is an indicator F * i (t) = 1 {t≥si} for some s i ∈ [0, T -τ i ].
To show this property, let us remark a few facts.

Lemma 1. For any CDF F i , there exist a staircase CDF F i such that

F i H i ≤ F i H i .
Proof: Since H i is right continuous, one can take a subdivision 0 = r 0 < r 1 < . . . < r n = T of [0, T ] such that H i is of constant sign on subintervals I k := [r k , r k+1 [, but changes sign between consecutive subintervals. Now define F i on I k as max I k F i if H i is negative on I k , and min I k F i otherwise (see Figure 1). This definition yields a new CDF such that

F i (t)H i (t) ≤ F i (t)H i (t) for all t ∈ [0, T ].
Lemma 2. For any staircase CDF F i , there exists a "one step" CDF F * i such that F * i H i ≤ F i H i . Proof: From the previous lemma, we can suppose without loss of generality that F i is a staircase CDF, so that F i (t) =

F i H i = k p k A k with A k := T r k H i (t)dt.
This is just a convex combination of real constants: if we denote by m the index of the minimum A k , s i := r m and

F * i (t) = 1 {t≥si} , then F * i H i = A m ≤ k p k A k ≤ F i H i .
Hence, given H i , there is an optimum F * i which is an indicator F * i (t) := 1 {t≥si} , where we know how to compute s i from H i . Furthermore,

F * i H i = si+τi si j =i d j φ j (t)dt,
which means the best move for player i is to schedule his job deterministically at a time minimizing the (weighted) sum of the probabilities of other jobs being active during his span

[s i , s i + τ i ].
This game seeks to minimize i γ i by optimizing each γ i iteratively. It does not necessarily lead to the optimal solution since re-scheduling job i may increase any γ j for j = i, yet we will see in section VI that it achieves its goal pretty well in practice.

IV. PARTIAL KNOWLEDGE SETTING

In this section, we suppose players do not share information about each others demands (for privacy reasons), but can still make inferences through the instantaneous total load λ(t) which is assumed actively communicated by the network. 2 We consider an iterative decision process where, at time t, user i decides (stochastically) whether to schedule his job or not according to:

• his own parameters (d i , τ i ),

• the past load profile {λ(t ), t < t}. Concerning the load profile, we will focus on protocols where the decision at time t only depends on the last know value of the load λ(t -).

In what follows, we suppose that all jobs' durations τ n are integer multiples of a unit time slot duration τ 0 dividing T , so that we can without loss of generality schedule jobs at times multiples of τ 0 .

A. ALOHA Strategy

The first strategy we propose is inspired by the slotted ALOHA protocol [START_REF] Abramson | The ALOHA System: another alternative for computer communications[END_REF]. At each time step, if his job has not been scheduled yet, player i applies the following decision procedure, which is parametrized by 0 < q i < p i < 1:

Algorithm IV.1 ALOHA decision procedure for player i if t = T -τ i (last possible scheduling slot) then s i ← t else if λ(t -) + d i ≤ L then s i ← t with probability p i else s i ← t with probability q i end if
Parameters p i should be low enough to avoid customers synchronization, but high enough to allow most of the jobs to induce no overage (keep λ < L).

When L is far below the mean load 1 T λ and all q i = 0, the policy may keep too many jobs for the end, resulting in peak loads at times close to T . Suitable values q i > 0 help deal with this unwanted behavior.

B. Decision Density

A way to generalize this approach is to set up a scheduling decision function for player i, g i (t) := g(λ(t), t, d i , τ i ) ∈ [0, 1], where we assume the form of g is the same for all players. Player i will therefore start at time t with probability g i (t), decisions being independently made by all players.

For example, under this formulation, the decision density for the ALOHA strategy is

g(λ, t, d, τ ) := mux(t = T -τ, 1, mux(λ + d ≤ L, p i , q i )),
where mux is the multiplexer function (mux(c, a 1 , a 2 ) := a 1 is c is true, and a 2 otherwise). Reasonable assumptions about g include:

• g increases with t;

• g → 1 when t → T -τ i ; • g < 1 when λ L and t T ; and • g decreases with λ when λ > L. With this in mind, we devised a new stochastic strategy improving the ALOHA one.

C. Time/Slackness Strategy

One of the issue of the ALOHA strategy lies in the way it discriminates jobs, since it focuses on the instantaneous load d i and only takes τ i into account as t → T -τ i . To remediate this shortcoming, we instead use the slackness σ defined when λ < L and t < T -τ as:

σ(λ, t, d, τ ) := dτ (L -λ)(T -τ -t) ,
i.e., the ratio of the job's overall energy consumption dτ and the residual energy (L -λ)(T -τ -t) which corresponds to the energy available with no overage under the assumption that λ stays constant. . We propose to use the simple density: [START_REF] Pipattanasomporn | Multi-Agent Systems in a Distributed Smart Grid: Design and Implementation[END_REF] so that just three parameters α, β and γ are in play. We call the associated policy Time/Slackness, since it consists of a BERNOULLI trial over g 1 (t) (ensuring the task is scheduled in time), followed by another trial based on slackness, giving a boost to the tasks for which there is enough residual energy.

g(t, σ) = g 1 (t) + (1 -g 1 (t))(β + γ • 1 {0<σ<1} ),
Experimental results (see section VI) confirm this new strategy yields better results than the ALOHA one, suggesting energy is a better discrimination criterion than power.

V. BLIND SETTING

In this section, we survey a power-grid setting where there is no communication layer between users. We also assume all customers have the same demand profile (d, τ ) and decide to schedule their jobs at times multiples of τ (where T = Kτ , K ∈ N) in a discrete time setting. We show that, in this simplified setting, the best strategy for customers is to choose their time slot uniformly at random.

Note that, here, broadcasting λ(kτ -) to the users would be useless since this value is independent from λ(kτ ).

Claim 2. The expected overall cost E[GC] is minimized when (independent) start times are chosen uniformly distributed on {kτ, k ∈ 0, K -1 }.

Proof: Let p k be the Probability Mass Function (PMF) of start-time s, common to all customers by symmetry. Limited information implies independent scheduling decisions. Therefore, the number of customers that select a given service epoch is binomially distributed, i.e.,

P[λ(kτ

) = nd] = N n p n k (1 -p k ) N -n .
So, the overall expected cost

N n=1 K-1 k=0 τ ndC L (nd) N n p n k (1 -p k ) N -n =: K k=1 G(p k )
is to be minimized subject to the PMF p in the K-dimensional simplex

K-1 k=0 p k = 1.
Note that G, defined by swapping order of summation, does not depend on the time-index k. The Lagrangian for this problem is

K-1 k=0 G(p k ) + c 1 - K-1 k=0 p k ,
with Lagrange multipier c, leading to the first-order necessary conditions whose solution is

∀k ∈ 0, K -1 , p k = (G ) -1 (c), i.e., p k is constant in k. (One can check that G is indeed bijective.) Condition K k=1 p k = 1 therefore yields p k = 1/K, so p is the PMF of a uniform distribution.
Recall that the conditions in which this uniform policy is optimal are different from the other settings we studied: we are here in a power-grid with no communication layer. We use this policy only as a reference in our numerical experiments.

VI. NUMERICAL EXPERIMENTS

Experiments on DR scenarios can involve embedding users in one of the IEEE test systems used in [START_REF] Liu | False Data Injection Attacks against State Estimation in Electric Power Grids[END_REF] (which can be found in MATPOWER, a MATLAB package) with a sharedresources game between them, taking into account on the characteristics of the buses. For our experiments we chose the simpler model, used in [START_REF] Mohsenian-Rad | Optimal and Autonomous Incentive-based Energy Consumption Scheduling Algorithm for Smart Grid[END_REF], of a local distribution network with one energy source and several load subscribers (see Figure 2). We implemented a simulator in PYTHON working on a six hours time frame divided into a customizable number of time slots. It implements the different policies we encountered:

• Uniform: the best solution in the blind setting;

• ALOHA I: the ALOHA strategy where all users share the same probabilities p i = p and q i = 0; • ALOHA II: same strategy with ∀i, p i = p > q i = q > 0;

• Time/Slackness: the policy from section IV with decision density (5) parametrized by α, β and γ; • Game: the game from section III.

Our simulator is open-source and available online at [START_REF] Dr Stratcomp | a Demand-Response Strategies Comparator[END_REF].

We set-up different test settings and ensured criterion (3) was met in each of them. For the Game policy, optimal behavior was reached for an average of 3 moves per player, which suggests this strategy converges quickly.

For the ALOHA and Time/Slackness policies, we manually chose good values of the parameters for each setting. In fact, all settings turned out to share approximately the same efficient values of the parameters, i.e.,

• ALOHA I: p ≈ 0.2 • ALOHA II: p ≈ 0.145 and q ≈ 0.0175 • Time/Slackness: α ≈ 45, β ≈ 0.006 and γ ≈ 0.12. This value of α implies time considerations are neglected while t < 90% T . In the last decile however, g 1 (t) yields more balanced schedules than a simple time-over check.

A. Residential Setting

The first scenario we considered is the case where all jobs have the same duration τ and instantaneous cost d, i.e., a residential area where houses have the same first-order load profile. For 1,000 users with a demand profile of 20 kW for 1 hour, the system's nominal load was set to L = 3, 000 kW, while we chose C 0 = 2.8 × 10 -6 $/kW/s (which is the first step in the model used at BC Hydro [START_REF] Bc Hydro | Electricity Rates[END_REF]) and C 1 = 2.8 × 10 -8 $/kW 2 /s. We focused on the global cost GC ramp experienced by the utility for each policy. Results averaged over several runs are shown in Figure 3 with confidence bars.

The best load profile for this configuration is a flat one. The Game strategy achieves a nearly optimal result, which comes from the fact that it is the only policy with enough information to actively seek a flat profile. Stochastic heuristics just try to approximate it (again with only limited information) while the uniform one tends to underload the borders (times close to 0 and T ). We also see that, in this setting where all customers are identical, the uniform policy yields better schedules than the heuristics from section IV.

B. Heterogeneous Setting

We also investigated the case where a lot of different profiles coexist on the network, including:

• a few "big" users demanding 100-400 kW for 2-5 hours,

• about 100 users demanding 10-50 kW for 1-3 hours, • about 100 users demanding 10 kW for about 1 hour, • a few "peak" users demanding > 800 kW for < 30 min. System-wide parameters were set to L = 1000 kW, C 0 = 2.8 × 10 -6 $/kW/s and C 1 = 2.8 × 10 -7 $/kW 2 /s. Results are show in Figure 4.

Again, the Game policy achieves the best behavior, but this time our heuristics perform better than the uniform strategy. Sample load profiles (which we won't produce here but are available online at [START_REF] Dr Stratcomp | a Demand-Response Strategies Comparator[END_REF]) indicate that:

• the Uniform strategy tends to make expensive mistakes, scheduling "big" players when the grid is already stressed and unloading the borders; • ALOHA I achieves a rather flat load, but is likely to keep big players for the end, yielding a final peak; • ALOHA II partially avoids this behaviour when q is high enough, but does not discriminate players in case of overage; • Time/Slackness is the best of the three heuristics and achieves a good compromise in scheduling both small and big users at each time step.

C. Peak-to-average ratio

We mentioned that the Game strategy is effective at peak shaving. We may illustrate this with figures, comparing its peak-to-average ratio (PAR) with the one of the Uniform strategy, averaged over several runs.

Game Uniform

VII. CONCLUSIONS

In this paper, we studied Demand Response problematics on multiple architectures for the dynamic pricing scheme [START_REF] Krishnan | Meters of Tomorrow[END_REF]. We saw that the general problem of finding an optimal schedule under this cost is NP-hard. We then surveyed different strategies depending on the degree of information sharing in the network.

When all demand information is shared, we proposed a game played by customers yielding good results in practice. When only the instantaneous load is known, we provided distributed strategies using the instantaneous load to reduce their costs. To experimentally evaluate all these policies, we developed our own open-source simulator which we released at [START_REF] Dr Stratcomp | a Demand-Response Strategies Comparator[END_REF]. Simulation results confirm that all these strategies perform better than when consumers do not communicate, especially the distributed game which significantly reduces the global cost and PAR, given the required information is available.

VIII. FUTURE WORK

In our jobs model, we chose a first-order approximation of the aggregate profile. Further study could take into account several devices per user with interruptible, non-constant load profiles and power consumptions ranging from 0.01 kW (light bulb) to 1 kW (dishwasher, cloths dryer). Furthermore, multiple hierarchy levels can be considered: devices, users, substations (trying to optimize their own costs, providing incentives to the customers).

Also, in cost (1) we took a threshold L constant within the timeframe of study, which does not encompass the integration of renewable energy sources that may deliver additional power during short time intervals. It would thus be interesting to survey how heuristics perform with a time-varying threshold L(t).

Finally, though the smart grid marketplace information flow is small in volume (and can therefore be "strongly" authenticated with low cost), security problems such as false data injection attacks [START_REF] Liu | False Data Injection Attacks against State Estimation in Electric Power Grids[END_REF] arise, which warrant additional consideration.

Fig. 1 .

 1 Fig. 1. Staircase CDF optimization.

Now define g 1 (

 1

Fig. 2 .

 2 Fig. 2. Setup of the test system.

Fig. 3 .

 3 Fig. 3. Average global costs in the residential setting.

Fig. 4 .

 4 Fig. 4. Average global costs in the heterogeneous setting.

TABLE I PAR

 I BENEFITS OF THE GAME POLICY.

Note that the data rates associated with this framework are so small (at most a few kilobits/second per user) that communication and security overhead is negligible.

We hence suppose that the utility is able to measure the effective state of the grid and compute its load, which is not a minor hypothesis since recent work[START_REF] Liu | False Data Injection Attacks against State Estimation in Electric Power Grids[END_REF] highlighted flaws in the state estimation techniques currently in use.


Ville N Pimenoff
email: vnpimenoff@iconcologia.net

Sara Tous

Yolanda Benavente

Laia Alemany

Francesc Wim Quint

Xavier Bosch

Ignacio G Bravo

Silvia De

Francesc Xavier Bosch

Silvia De Sanjosé

Distinct geographic clustering of oncogenic human papillomaviruses multiple infections in cervical cancers: Results from a worldwide cross-sectional study

Keywords: multiple HPVs infections, invasive cervical cancer, geographic clustering Additional

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Distinct geographic clustering of oncogenic human papillomaviruses multiple infections in cervical cancers: Results from a worldwide cross-sectional study Background Among the more than 350 different Papillomaviruses (PVs) identified (https://pave.niaid.nih.gov), above 200 are Human Papillomaviruses (HPVs). According to their carcinogenic potential, twelve HPVs are acknowledged as carcinogenic to humans (IARC Group I: HPV16/18/31/33/35/39/45/51/52/56/58/59), while 13 additional HPVs are classified as probably or possibly carcinogenic to humans (IARC Groups 2a and 2b: HPV26/30/34/53/66/67/68/69/70/73/82/85/97). All of them are usually referred to as high-risk (HR) HPVs. Two HPVs cannot be classified as to their carcinogenicity to humans (IARC Group 3: HPV6 and 11) and are referred to as low-risk (LR) HPVs. [START_REF] Bouvard | A review of human carcinogens-part B: biological agents[END_REF] Other viruses such as HPV40, 42, 43, 44, 54, and 74 have not been included in any of the IARC Groups but classified by IARC [START_REF] Iarc | Human papillomaviruses[END_REF] as LR-HPVs. All of these above mentioned HPV types belong to Alphapapillomaviruses within the Papillomaviridae family.

HPVs are the most common infectious agents in sexually transmitted infections: above 80% of sexually active women are infected by any HPVs during their lifetime. [START_REF] Einstein | Clinician's guide to human papillomavirus immunology: knowns and unknowns[END_REF] Fortunately, most of these HPVs infections are transient and cleared by the immune system, while persistent infection with HR-HPVs is a necessary cause for invasive cervical cancer (ICC) development. [START_REF] Walboomers | Human papillomavirus is a necessary cause of invasive cervical cancer worldwide[END_REF] Coinfections by multiple HPVs are rare in women with normal cytology (between 1% and 3%), common in women with precancerous lesions (between 15% and 41%) and display intermediate levels in HPV-induced ICCs with coinfections being observed between 6% and 8% of the cases. [START_REF] De Sanjose | Human papillomavirus genotype attribution in invasive cervical cancer: a retrospective crosssectional worldwide study[END_REF] The simultaneous infection by different HPVs seems to have an impact on the natural history of the viral infection. Thus, when compared to single HPV infections, cancerous lesions with multiple HR-HPVs infections contain more oncogenic E6/E7 HPV mRNA, often from more than one HR-HPV, thus suggesting that HPVs coinfections in precancerous lesions might be more susceptible to carcinogenesis than single infections. [START_REF] Andersson | Typedependent E6/E7 mRNA expression of single and multiple high-risk human papillomavirus infections in cervical neoplasia[END_REF] Additionally, women coinfected with at least two different HPVs in the cervix have a higher risk of high-grade cervical lesions, of ICC or of acquiring another HPV infection, compared to those infected with a single type. [START_REF] Wentzensen | New Mexico HPV pap registry steering committee. No evidence for synergy between human papillomavirus genotypes for the risk of high-grade squamous intraepithelial lesions in a large population-based study[END_REF][START_REF] Carrillo-García | Impact of human papillomavirus coinfections on the risk of high-grade squamous intraepithelial lesion and cervical cancer[END_REF][START_REF] Dickson | Cervical cytology and multiple type HPV infection: a study of 8182 women ages 31-65[END_REF][START_REF] Spinillo | Multiple human papillomavirus infection with or without type 16 and risk of cervical intraepithelial neoplasia among women with cervical cytological abnormalities[END_REF][START_REF] Adler | Cumulative impact of HIV and multiple concurrent human papillomavirus infections on the risk of cervical dysplasia[END_REF] Current consensus on the etiology of HPVs-related cancers is that individual lesions are associated with individual HPV infections, meaning that lesions, and thus HPVs-induced cancers, are clonal. [START_REF] Murray | Possible clonal origin of common warts (Verruca vulgaris)[END_REF][START_REF] Quint | One virus, one lesion-individual components of CIN lesions contain a specific HPV type[END_REF] Accordingly, given the large differences in prevalence of HPVs coinfections in lesions of different grade and grade combinations, the etiological role of multiple HPVs coinfections in the development of virus-induced cancers, if any, still needs to be evaluated. This is especially true given the increasing number of large-scale HPVs vaccination programs currently implemented using bivalent vaccine protective against HPV16/18/31/33/35/45 types, and further inducing herd effects against HPV18/31/33/35 types. [START_REF] Lehtinen | Gender-neutral vaccination provides improved control of human papillomavirus types 18/31/33/35 through herd immunity: results of a community randomized trial (III)[END_REF] In parallel, the 9-valent HPV vaccine is protective against seven HR-HPV types (HPV16/18/31/33/45/52/58) and when these vaccines are widely implemented they will thoroughly change the prevalence pattern of the targeted HR-HPVs infections in humans.

The aim of our study was to explore the prevalence and patterns of HPVs coinfections in a large worldwide crosssectional study of ICC cases using standardized HPV genotyping methods and applying different statistical models.

Methods

Study design

Our study is based on an international retrospective crosssectional study of presence of HPVs in ICC cases as previously described. [START_REF] De Sanjose | Human papillomavirus genotype attribution in invasive cervical cancer: a retrospective crosssectional worldwide study[END_REF] The original study included information from 10,575 ICCs and HPV DNA was detected in 85% of the cases. Briefly, paraffin-embedded specimens from consecutive cases (aged 16-97 years) diagnosed with ICC were obtained from histological archives from 38 countries worldwide. Our study and all protocols were approved by both local and ICO ethics committees. After pathological confirmation of primary ICC, HPV DNA detection was performed by PCR with SPF-10 broad-spectrum primers followed by DNA enzyme immunoassay (DEIA) and genotyping with a reverse hybridization line probe assay (LiPA25 SPF10-PCR)-LIPA25 that detects 25 HPV types (HPV6, 11, 16, 18, 31, 33, 34, 35, 39, 40, 42, 43, 44, 45, 51, 52, 53, 54, 56, 58, 66, 68, 70, and 74) (see details in Ref. [START_REF] De Sanjose | Human papillomavirus genotype attribution in invasive cervical cancer: a retrospective crosssectional worldwide study[END_REF] ). Cases showing DEIA positive but LIPA25 negative results (undetermined HPV) were further subjected to Sanger sequencing for HPV identification, as previously described. [START_REF] De Sanjose | Human papillomavirus genotype attribution in invasive cervical cancer: a retrospective crosssectional worldwide study[END_REF] Special care was given to maintain stringent conditions for HPV type-level detection using the SPF10-PCR LiPA25, a method particularly suitable for HPV genotyping from FFPE specimens. [START_REF] De Sanjose | Human papillomavirus genotype attribution in invasive cervical cancer: a retrospective crosssectional worldwide study[END_REF][START_REF] Pérez | Comparison of 2 different PCR-based technologies for the detection of human papilloma virus from paraffin-embedded tissue: genómica clinical arrays versus SPF(10)-LiPA(25)[END_REF] To further prevent contaminations the paraffin blocks in our study were processed under strict conditions: four paraffin sections were systematically obtained from each block. First and fourth sections were used for histopathological assessment after hematoxylin and eosin staining, and the second and

Results

Demographic characteristics of the ICCs

From the 10,575 ICCs previously reported, [START_REF] De Sanjose | Human papillomavirus genotype attribution in invasive cervical cancer: a retrospective crosssectional worldwide study[END_REF] we excluded the ICC cases with a negative HPV DNA test (n = 1598) and cases containing DNA from undetermined HPV further tested with Sanger sequencing (n = 197). Thus, the dataset used in our study contained 8780 ICC HPV-positive cases, most of them (91.8%; N = 8064) histologically diagnosed as squamous cell carcinomas (SCCs, Table 1). The final dataset of HPVpositive ICC cases consisted of women with a mean age of 51 years (range 17-97) at diagnosis and 23.1% recruited from Europe, 39.6% from the Americas, 31.3% from Asia including Oceania, and 6% from Africa (Table 1). Since infection with other pathogens, notably Human Immunodeficiency Virus, could have an impact on age at ICC onset, we verified that the mean age at ICC diagnosis exclusively for African women did not differ from that for women outside Africa (p = 0.240, data not shown). Most (n = 8193) HPV-positive ICC cases showed the presence of a single HPV, and in 83.8% (N = 492) of the cases with multiple infections only two distinct HPVs were detected (Supporting Information Fig. S1). Logistic regression analysis among all ICCs showed prevalence differences between single and multiple HPVs infected ICCs by geographic origin and by period of diagnosis (Table 1). By geographic region, prevalence of multiple HPVs infections was consistently about three-fold less frequent (p < 0.001) among women from Europe (0.36 [95% CI: 0.27-0.49]), Asia (0.31 [95% CI: 0.23-0.42]) or the Americas (0.36 [95% CI: 0.27-0.49]) compared to women from Africa (Table 1). Stratification of the African ICC cases revealed that 41.7% were from Mozambique, 30.1% from Nigeria, 19.6% from Uganda and 8.5% from Algeria (data not shown). Cases diagnosed between 1990 and 2009 showed 7.3% prevalence of multiple infections, 1.32-fold (95%CI: 1.08-1.60) higher than ICCs diagnosed between 1940 and 1989 (5.8% prevalence) (Table 1). However the difference in prevalence for multiple HPVs infections was not significant when stratifying ICCs into cases collected before and after 2000 (results not shown).

HPV type-specific distribution

HPV type-specific distribution in ICCs is presented in Figure 1. With 61.9% (N = 5432) prevalence, HPV16 was the most frequent HPV, followed by HPV18 (11.1% prevalence; N = 977) and HPV45 (6.7% prevalence; N = 585). Their respective contributions in coinfected ICCs were also the highest: HPV16, 58.4% (N = 343), HPV18, 21.0% (N = 123); and HPV45, 18.9% (N = 111) (Supporting Information Table S2). According to viral taxonomy, the Alpha-9 HPVs (which include HPV16) were involved in 86.5% (N = 506) of the coinfected ICCs, while Alpha-7 HPVs (which include HPV18) were involved in 44.6% (N = 262) of the coinfected ICC cases (see also Fig. 1). third sections were used for HPV DNA detection and genotyping as described. [START_REF] De Sanjose | Human papillomavirus genotype attribution in invasive cervical cancer: a retrospective crosssectional worldwide study[END_REF] Importantly, a blank paraffin section was cut and processed in-between specimens to control for any carryover contamination in addition to the standard controls.

Statistical analyses

ICC cases were classified into single or multiple HPV-infected cases based on the above described genotyping results. Appropriate statistical tests were used to compare the proportion of single and multiple HPVs infections by the design variables (age at diagnosis, geographical origin, period of diagnosis and histological type). We estimated prevalence ratio and 95% bootstrapped-confidence interval of single/multiple (S/M) infections for each individual HPV. Unconditional logistic regression models were used to estimate association between demographic and clinical characteristic of the cases, probabil-ity of multiple infections, and to identify pairwise coinfections of HPVs that were more likely to be detected together com-pared to being in combination with any other HPV in multi-ple infections. To simultaneously estimate the possible associations beyond pairwise HPVs comparison among multi-ple HPV infected ICCs collected worldwide, we fitted the his-tological type, period of diagnosis and geographical origin of the sample into our generalized linear model (GLM) of the multivariate HPVs multiple infections dataset. [START_REF] Wang | Mvabund-an R package for model-based analysis of multivariate abundance data[END_REF] Further, a univariate tests for each HPV in coinfected ICCs were imple-mented within multivariate analysis to estimate whether there were typespecific differences in HPV coinfection distribution by histology, period of diagnosis or geographic origin. The observed pairwise and composition level HPVs coinfection pattern differences were visualized using network plots constructed with igraph R package (R Core Team, http://www.Rproject.org/). Classification and regression trees (CART) analyses [START_REF] Breiman | Classification and regression trees[END_REF] were performed using the rpart R package to estimate and visualize the possible HPVs coinfection patterns beyond pairwise HPVs comparison, accounting for age at diagnosis, histology, period of diagnosis and geographic origin of the cases. In practice, CART analyses use recursive partitioning to structure the data into a series of dichotomous splits, e.g. presence or absence of each HPV type and other demo-graphic variables, to create a decision tree with the goal of correctly classifying clusters of the above mentioned interact-ing variables, for example, HPV16 or HPV18 coinfected ICCs. Each independent variable is examined and a split is intro-duced to maximize the sensitivity and specificity, calculated using the predicted probability from the multivariable logistic regression of the classification, and resulting in a decision tree. [START_REF] Breiman | Classification and regression trees[END_REF] All regression models were adjusted by the design vari-ables and Bonferroni's correction was implemented in pair-wise analyses to control for multiple comparisons. All statistical tests were two-sided and rejecting the null hypothe-sis was set at a p value of less than 0.05. Analyses were carried out with Stata software (Release 13, StataCorp, USA) except otherwise specified.

Prevalence of single and multiple HPVs infections in ICC cases along with the patient characteristics, stratified by HR-and LR-HPVs, is given in Supporting Information Table S1. Among ICCs with multiple infections (n = 587), 76.0% were exclusively infected with HR-HPVs while only two cases were exclusively infected with LR-HPVs (HPV54|74 and HPV6|53|66). The remaining 23.7% of the multiple infected ICC cases showed a combination of HR-and LR-HPVs. The corresponding values for ICCs with single infections were 99.5% (N = 8149) for only HR-HPVs and 0.5% (N = 44) for only LR-HPVs (Supporting Information Table S1). Altogether, six (HPV16/18/31/33/35/45) and seven (HPV16/18/31/33/45/52/58) HR-HPVs covered/crossprotected by the 2-and 9-valent HPV vaccines, accounted for respectively 87.9% and 90.9% of the single-infected ICC cases and were present in 92.0% and 95.9% of the multiple HPV coinfected ICCs. 2). A univariate test applied to the multivariate analysis showed a significant virus-specific level difference by histology for HPV31 and by period of diagnosis for HPV74 (Table 2). Furthermore, seven HR-HPVs (HPV18, HPV35, HPV51, HPV52, HPV56, HPV58 and HPV74; Table 2) showed significant or nearly significant coinfection prevalence differences by geography. Accordingly, crude prevalence estimates of the above listed seven HR-HPVs in coinfected ICCs from Africa were either about two-fold higher, or at least five-fold lower than in any other world region (Table 2). In agreement, stratifying by histology and restricting the multivariate analysis on the SCC cases showed significant HPV coinfection pattern differences by period of diagnosis (LRT 41.60, p = 0.012) and by geographic origin (LRT 193.59, p = 0.001). Univariate test revealed also significant prevalence differences by period of diagnosis for HPV74 and by geography for five HR-HPVs (Supporting Information Table S3). HPV type-specific networks shown in Figure 3 stratified by geographic origin visualize the multiple HPVs coinfection combinations observed and more importantly, the significant higher frequency of HPV18, HPV35, HPV51 and S2. Virtually all HR-HPVs, except HPV66 and HPV70, showed S/M prevalence ratios above one indicating a higher prevalence in single than in multiple infections. In contrast, all LR-HPVs (except HPV42) and particularly Alpha-10 LR-HPVs showed S/M prevalence ratios below one, indicating a higher prevalence in multiple than in single infections.

HPV coinfection pattern analysis

Among all HPVs pairwise combinations detected in ICCs with multiple infections, the prevalence of eight HPV pairs was significantly different from expected by chance: HPV16|74, HPV31|33, HPV31|44, HPV33|44 and HPV45|70 pairs were more likely to be detected together than with any other HPV type, while HPV16|31, HPV16|45 and HPV16|51 pairs were less likely to be detected together compared to other HPVs (Fig. 2a). Restricting the analysis to ICCs infected only with two HPVs revealed that the HPV45|70 combination was more likely and HPV16|45 and HPV18|45 less likely to be detected together compared to any other combination of HPV types (Fig. 2b).

To explore individual viral prevalence differences in ICCs with multiple HPVs infections beyond pairwise analysis, we HPV52 present in multiple infections in samples from Africa compared to samples outside Africa. The HPV pairs identified with significantly different pairwise coinfections by logistic regression are also visualized in Figure 3.

Finally, two different recursive partitioning analyses were performed to explore the appearance of HPVs coinfections in ICCs with either HPV16 (Fig. 4) or HPV18 (Supporting Information Fig. S2) and accounting for age at diagnosis, geographic origin and period of diagnosis of the cases. First CART analysis searched for a homogenous split of the data on the chosen variable (i.e. HPV16) and the subsequent best splits among all variables that partitioned the data in each step into further most homogenous subgroups, accounting also for age at diagnosis, period of diagnosis and geographic origin of the cases. The predicted probabilities (PP) for each split for other HPVs occurring with HPV16 (PP = 58.4% [95%CI: 54.3-62.5], Fig. 4) or HPV18 (PP = 21.0%, Supporting Information Fig. S2) were presented. Second and third most likely splits suggest respectively a major exclusion of HPV31 (PP = 62.6% [95%CI: 58.3-66.9]) and of HPV45 (PP = 67.9% [95% CI: 63.1-72.5]) from HPV16-related coinfections. In parallel, the two major exclusion variables from HPV18 coinfections and their predicted probabilities were non-African cases (PP = 39.6%) and HPV45 coinfection exclusion (PP = 47.4%) (Supporting Information Fig. S2).

Discussion

In our study, we describe distinct clustering patterns of different HPVs detected in multiple HPVs infections in the largest available global collection of ICCs in unvaccinated women. Our results confirm that the simultaneous detection of different HPVs in ICCs is rare albeit far from negligible (6.7% in our dataset, N = 8780). The association patterns identified suggest further that the clustering of multiple HPVs in ICCs is not random. Our most relevant findings were a significant causal/coincidential infection of LR-HPV types 44 or 74 in ICCs in combination with HR-HPV types 16, 31 or 33, and a marked geographic clustering of particular oncogenic HPV types 51, 35, 18 and 52 (in descending significance) in multiple HPVs infected cervical cancers in African women compared to women outside Africa. A number of studies have previously communicated that clustering patterns observed among HPVs do not always occur at random in multiple infections detected in normal cytology, in cervical precancerous lesions or in cancer lesions. [START_REF] Wentzensen | New Mexico HPV pap registry steering committee. No evidence for synergy between human papillomavirus genotypes for the risk of high-grade squamous intraepithelial lesions in a large population-based study[END_REF][START_REF] Vaccarella | Patterns of human papillomavirus types in multiple infections: an analysis in women and men of the high throughput human papillomavirus monitoring study[END_REF][START_REF] Chaturvedi | Human papillomavirus infection with multiple types: pattern of coinfection and risk of cervical disease[END_REF][START_REF] Soto-De Leon | Distribution patterns of infection with multiple types of human papillomaviruses and their association with risk factors[END_REF][START_REF] Goldman | Patterns of cervical coinfection with multiple human papilloma virus types in a screening population in Denmark[END_REF][START_REF] Mollers | Population-and type-specific clustering of multiple HPV types across diverse risk populations in The Netherlands[END_REF] Nevertheless, our results disagree with the previous studies, which suggest that the observed nonrandom HPVs clustering in multiple HPVs infections is mainly due to technical biases originating from the HPV genotyping methods. [START_REF] Vaccarella | Concurrent infection with multiple human papillomavirus types: pooled analysis of the IARC HPV prevalence surveys[END_REF] The existence of genuine nonrandom association pattern of HPVs in multiple infections is elusive, and identifying and excluding technology-driven biases is not a simple task. Nevertheless, a public health concern is that HPVs coinfections may have a role in ICCs development not only among unvaccinated [START_REF] Schmitt | VALGENT study group. Multiple human papillomavirus infections with high viral loads are associated with cervical lesions but do not differentiate grades of cervical abnormalities[END_REF] but also among HPV vaccinated women. [START_REF] Gray | Evaluation of HPV type-replacement in unvaccinated and vaccinated adolescent femalespost-hoc analysis of a community-randomized clinical trial (II)[END_REF] The rational is that due to ecological competition between HPV types, the HPV vaccine implementation may further clear the ecological niche for the nonvaccine-targeted oncogenic HPVs. [START_REF] Gray | Evaluation of HPV type-replacement in unvaccinated and vaccinated adolescent femalespost-hoc analysis of a community-randomized clinical trial (II)[END_REF][START_REF] Mesher | Population-level effects of human papillomavirus vaccination programs on infections with nonvaccine genotypes[END_REF][START_REF] Tota | Evaluation of type replacement following HPV16/18 vaccination: pooled analysis of two randomized trials[END_REF] Particularly, some HPVs such as HPV39 and HPV51 not included in the vaccine formulations have previously shown increase in prevalence after vaccination suggesting the possibility of type replacement. [START_REF] Gray | Evaluation of HPV type-replacement in unvaccinated and vaccinated adolescent femalespost-hoc analysis of a community-randomized clinical trial (II)[END_REF][START_REF] Mesher | Population-level effects of human papillomavirus vaccination programs on infections with nonvaccine genotypes[END_REF] However, thus far, competition between HPV types, which is the prerequisite for HPV type replacement, has not been persistently observed among unvaccinated [START_REF] Wentzensen | New Mexico HPV pap registry steering committee. No evidence for synergy between human papillomavirus genotypes for the risk of high-grade squamous intraepithelial lesions in a large population-based study[END_REF][START_REF] Vaccarella | Patterns of human papillomavirus types in multiple infections: an analysis in women and men of the high throughput human papillomavirus monitoring study[END_REF][START_REF] Mejlhede | Multiple human papilloma virus types in cervical infections: competition or synergy?[END_REF][START_REF] Merikukka | Differences in incidence and co-occurrence of vaccine and nonvaccine human papillomavirus types in Finnish population before human papillomavirus mass vaccination suggest competitive advantage for HPV33[END_REF] or vaccinated populations. [START_REF] Gray | Evaluation of HPV type-replacement in unvaccinated and vaccinated adolescent femalespost-hoc analysis of a community-randomized clinical trial (II)[END_REF][START_REF] Mesher | Population-level effects of human papillomavirus vaccination programs on infections with nonvaccine genotypes[END_REF][START_REF] Tota | Evaluation of type replacement following HPV16/18 vaccination: pooled analysis of two randomized trials[END_REF]

HPVs coinfection prevalence

Here in our study we describe known as well as novel patterns of nonrandom HPVs clustering in global ICC cases, identified using a single well established HPV detection and genotyping technique. For clarity, it is important to remind that the current consensus of ICC etiology follows the logic of superinfection exclusion, [START_REF] Biryukov | Superinfection exclusion between two high-risk human papillomavirus (HPV) types during a co-infection[END_REF] with a subsequent clonal origin of HPVinduced lesions indicating that a single HPV virus contributes to the development of a single ICC lesion. [START_REF] Murray | Possible clonal origin of common warts (Verruca vulgaris)[END_REF][START_REF] Quint | One virus, one lesion-individual components of CIN lesions contain a specific HPV type[END_REF] This is to say, that the presence of DNA from two or more HPV types in an ICC sample does not necessarily mean that the different HPVs detected, actually contribute to the cancer development. Previously, techniques of laser capture microdissection (LCM) have been successfully used to affirm the likely involvement of causal or coincidential HPVs with particular lesion component in multiple HPVs infected cases, further supporting the hypothesis of clonal origin of HPV-induced cervical cancers. [START_REF] Quint | One virus, one lesion-individual components of CIN lesions contain a specific HPV type[END_REF][START_REF] Guimerà | Laser capture microdissection shows HPV11 as both a causal and a coincidental infection in cervical cancer specimens with multiple HPV types[END_REF] In ICC cases worldwide, we did not observe significant differences in HPVs coinfection prevalence between the age groups at diagnosis. These results seem to be in contrast to the studies reporting higher prevalence of multiple HPVs infections in sexually active adolescents and in young women, [START_REF] Chaturvedi | Human papillomavirus infection with multiple types: pattern of coinfection and risk of cervical disease[END_REF][START_REF] Schmitt | VALGENT study group. Multiple human papillomavirus infections with high viral loads are associated with cervical lesions but do not differentiate grades of cervical abnormalities[END_REF][START_REF] Nielsen | Type-specific HPV infection and multiple HPV types: prevalence and risk factor profile in nearly 12,000 younger and older Danish women[END_REF] and thus, our results might actually indicate that regardless of the sexual risk-taking behavior coinfection by multiple HPVs may be a risk factor for malignant cervical transformations. Interestingly however, a recent study suggests reduced HPV-driven carcinogenesis in cases of coinfection with oncogenic HPV types from the same papillomavirus species among HIV-infected patients. [START_REF] Sobota | Co-infection with HPV types from the same species provides natural cross-protection from progression to cervical cancer[END_REF] This was interpreted as an evidence of natural cross-protection against HPV-induced carcinogenesis by multiple HPV infections from the same HPV species group. [START_REF] Mollers | Population-and type-specific clustering of multiple HPV types across diverse risk populations in The Netherlands[END_REF][START_REF] Sobota | Co-infection with HPV types from the same species provides natural cross-protection from progression to cervical cancer[END_REF][START_REF] Verdenius | Cross protection against HPV might prevent type replacement[END_REF] Nonetheless, we did observe that HPVs coinfection prevalence was consistently about three-fold higher in ICCs from Africa compared to the rest of the world. One explanation could be concomitant HIV infection which may increase the prevalence of multiple HPVs coinfections. [START_REF] Chaturvedi | Prevalence and clustering patterns of human papillomavirus genotypes in multiple infections[END_REF] Information on the sexual risk-taking behavior or HIV status was, unfortunately, not available for our ICC cases but the corresponding population estimates by Bruni and colleagues [START_REF] Bruni | Summary report: human papillomavirus and related diseases in Africa. ICO|IARC Information Cent HPV Cancer (HPV Inf Centre) Hum Papillomavirus Relat[END_REF] showed a relatively high prevalence of HIV infections among women in Mozambique, Uganda and Nigeria (11.3%, 7.2 and 3.7%, respectively). Thus, the larger fraction of women infected by HIV in these African populations compared to other populations worldwide may partly play a role in the higher presence of multiple HPVs infections also in our African ICC cases.

DNA degradation

Despite using a highly sensitive PCR based method for HPV genotyping (based on a 65 base pairs amplicon), we detected a minor 1.3-fold higher prevalence of coinfections in the more recently diagnosed ICC cases (between 1990 and 2009) than in older cases (between 1940 and 1989), but the difference was not significant when ICCs were stratified for cases collected before and after 2000 (results not shown), and these results were further supported by our multivariate analysis showing significance only for the 1990 stratification. Essentially, most paraffin embedded samples are initially fixed in formalin-containing solution, which reacts with the lateral amino groups in proteins and DNA, and thus, complicates the DNA isolation from FFPE samples. [START_REF] Gilbert | The isolation of nucleic acids from fixed, paraffinembedded tissues-which methods are useful when?[END_REF] More importantly, the aldehyde in the formalin is unstable if unbuffered (used in most historical FFPE samples at least before the 1990s [START_REF] Gilbert | The isolation of nucleic acids from fixed, paraffinembedded tissues-which methods are useful when?[END_REF] ), and disproportionate rendering methanol and formic acid, causing extra fragmentation of the sample DNA. In short, we cannot completely rule out the possible effect of additional DNA degradation stemming from older tissue-fixation protocols, as described by Gilbert and colleagues. [START_REF] Gilbert | The isolation of nucleic acids from fixed, paraffinembedded tissues-which methods are useful when?[END_REF] Notwithstanding, we estimate that older tissue-fixation protocol related DNA degradation plays a minor role (e.g. lowest LRT estimated for period of diagnosis) in explaining the higher prevalence of multiple HPVs infections in cases collected since 1990 compared to the older ICC cases.

HPV type-specific clustering

At pairwise association level, we found a significant causal or coincidental infection of LR-HPV types 44 or 74 in combination with HR-HPV types 16, 31 or 33 in global ICCs. Moreover, all Alpha-10 LR-HPVs (HPV6/11/44/74) showed a markedly higher prevalence in multiple HPV type infections compared to single HPV type infected ICCs.

Altogether, we identified five positive and four negative significant associations between HPV types among the multiple infected ICCs. Three of the five positive associations were novel (i.e. HPV16|74, HPV31|44 and HPV33|44) while association between Alpha-9 species HPV31 and HPV33 or between Alpha-7 species HPV45 and HPV70 had also been previously estimated. [START_REF] Mendez | Instituto Nacional de Cancerologia human papillomavirus study group. Cervical coinfection with human papillomavirus (HPV) types and possible implications for the prevention of cervical cancer by HPV vaccines[END_REF] Current epidemiological research supports the double etiological origin for HPV-induced cervical lesions: LR-HPVs, mainly HPV6 and HPV11, mostly leading to benign low-grade squamous intraepithelial lesions (LSILs), and oncogenic HR-HPVs, if persistent, mostly leading to high-grade SILs, and considered cervical cancer precursor lesions. Nonetheless, low-risk HPVs such as HPV6, 11 and 44 from the Alpha-10 species have previously been occasionally associated with causal and coincidental infections in anogenital cancers. [START_REF] Guimerà | Laser capture microdissection shows HPV11 as both a causal and a coincidental infection in cervical cancer specimens with multiple HPV types[END_REF][START_REF] Guimerà | The occasional role of low-risk human papillomaviruses 6, 11, 42, 44, and 70 in anogenital carcinoma defined by laser capture microdissection/PCR methodology: results from a global study[END_REF] In addition, previous studies have also shown that a fraction of low-grade cervical SILs may also harbor components of HSIL or carcinoma often associated with independent HPVs infections, thus following the logic of "one-virus onelesion". [START_REF] Quint | One virus, one lesion-individual components of CIN lesions contain a specific HPV type[END_REF][START_REF] Herfs | A discrete population of squamocolumnar junction cells implicated in the pathogenesis of cervical cancer[END_REF] In agreement, our results match well the logic of an essentially clonal origin of HPV-related lesions and support the occasional causal and/or coincidental role of low-risk Alpha-10 HPVs in cervical cancer development.

Theoretical models have predicted natural competition between HPV types, [START_REF] Mejlhede | Multiple human papilloma virus types in cervical infections: competition or synergy?[END_REF][START_REF] Murall | Revising ecological assumptions about human papillomavirus interactions and type replacement[END_REF] such as between types from the same HPV species group, and here we observed four significant negative pairwise HPV type associations: HPV16|HPV31 and HPV18|HPV45 pairs within the same species, and HPV16| HPV45 and HPV16|HPV51 pairs with longer pairwise phylogenetic distance between the types. Furthermore, our multivariate CART analysis beyond pairwise HPVs estimation revealed the same patterns of negative association between three of the HPV pairs, and thus, sustained the idea of possible competition or antagonistic interaction between these types. Novel negative association between HPV16|HPV45 types and two also previously reported negative associations, namely HPV16|HPV31 and HPV18|HPV45 pairs within the same species group. [START_REF] Mejlhede | Multiple human papilloma virus types in cervical infections: competition or synergy?[END_REF][START_REF] Xi | Human papillomavirus types 16 and 18 DNA load in relation to coexistence of other types, particularly those in the same species[END_REF] Notably, antagonistic interactions between HPV6/11 and HPV16, and between HPV6 and HPV31 types, earlier infections with the LR-HPVs protecting against cervical cancer risk associated with the latter HR-HPV adenocarcinomas (ADC) or mixed adenosquamous carcinomas (ADSC). Previous studies have already reported these HPV type distribution differences between single infected cervical SCCs and ADCs/ADSCs, [START_REF] De Sanjose | Human papillomavirus genotype attribution in invasive cervical cancer: a retrospective crosssectional worldwide study[END_REF][START_REF] Li | Human papillomavirus type distribution in 30,848 invasive cervical cancers worldwide: variation by geographical region, histological type and year of publication[END_REF] and our results among multiple HPV coinfected ICCs further sustain the results.

Nevertheless, and despite the large worldwide crosssectional sampling of ICC cases and the standardized HPV genotyping methodology used, our study suffers from certain limitations. Firstly, detection of more than one HPV type or variant in a sample is not enough to assign causality of the multiple HPVs with the particular lesion, in this case, the ICC. Secondly, further cohort studies and use of other biomarkers providing biological activity of particular HPVs such as mRNA or p16INK4a expression, and with new statistical models are needed to better disentangle the role of multiple HPVs infections in the etiology of cervical cancer, especially in the post-HPV-vaccinated environment. Thirdly, a better characterization of the patients' genetic background, possible HIV infections, heterogeneous ICCs biopsy material and detection method cross-validation is needed to confirm the potential role of multiple HPVs infections in ICCs. Lastly, the detection and ultimately the likely role of multiple HPVs infections at HPV variant level has only recently been verified thanks to the next-generation sequencing applied to HPV genotyping, [START_REF] De Oliveira | High-level of viral genomic diversity in cervical cancers: a Brazilian study on human papillomavirus type 16[END_REF] but within host HPV type-level variation is beyond the scope of our study.

Our genotyping strategy was based on the differential hybridization of very short PCR products, and as such, we cannot completely rule out the possibility that some of the observed clustering of multiple HPV type coinfections in our study is not partly due to PCR and hybridization biases and directional cross-detection of multiple HPVs. However, several lines of evidence argue against that the possible technology-based cross-detection biases are the main explanation for the HPVs clustering patterns observed. First, it has been previously shown that the SPF10-PCR LiPA25 methodology is more sensitive and less prone to overemphasize multiple HPVs infections than other PCR-based methods. [START_REF] Pérez | Comparison of 2 different PCR-based technologies for the detection of human papilloma virus from paraffin-embedded tissue: genómica clinical arrays versus SPF(10)-LiPA(25)[END_REF] Second, it is unlikely that the observed differential HPVs clustering after histological or geographical stratification exclusively originates from unspecific genotyping, as it would imply that such hybridization biases are both histology and geography dependent. Third, we found the worldwide nonrandom clustering of different HR-HPVs in coinfected ICCs at two different analysis depths, at pairwise and at multivariate level as described above.

Taken together, HPVs coinfection patterns are rare in ICCs, and thus difficult to study at population-level. In our study, however, we provide robust and worldwide estimations on the proportion of single and multiple HPV type infected ICCs and reveal HPVs pairwise combinations significantly deviating from the null expectation. Further, we identified for the first time beyond pairwise analysis a distinct coinfection pattern in African ICCs compared to any other world region.

infections have been documented, [START_REF] Luostarinen | No excess risk of cervical carcinoma among women seropositive for both HPV16 and HPV6/11[END_REF][START_REF] Luostarinen | Order of HPV/chlamydia infections and cervical high-grade precancer risk: a case-cohort study[END_REF] and indirectly fit with our observations on less likely coinfections with HPV16 and 31, and HPV18 and 45. Taken together, the assumption that HPV types within the same species group, such as Alpha-7 (HPV18, 39, 45, 59, 68, 70 and 85) or Alpha-9 (HPV16, 31, 33, 35, 52, 58 and 67), compete with one another [START_REF] Mejlhede | Multiple human papilloma virus types in cervical infections: competition or synergy?[END_REF][START_REF] Xi | Human papillomavirus types 16 and 18 DNA load in relation to coexistence of other types, particularly those in the same species[END_REF] is not clear as we also observed positive pairwise HPV type associations for HPV31|33 and HPV45|70 pairs.

Our most relevant finding is the distinct geographic clustering of HPVs multiple infections in ICCs among women from Africa compared to women outside Africa. Particularly, a significant geographic clustering of oncogenic HPV types 51, 35, 18 and 52 in multiple infections was observed in African women. That is, beyond pairwise level association and exclusively among the HPVs coinfected ICCs, we identify for the first time seven (i.e. HPV51, HPV35, HPV18, HPV74, HPV58, HPV52, HPV56) geographic clustering patterns for specific HPVs, in descending significance. Earlier studies regarding worldwide HPV genotype attribution in ICCs have described particular type-specific geographic structure, namely higher prevalence for HPV16, HPV18 and HPV45 in African populations, and a higher prevalence for HPV52 and HPV58 types in Asian populations compared to other world regions. [START_REF] De Sanjose | Human papillomavirus genotype attribution in invasive cervical cancer: a retrospective crosssectional worldwide study[END_REF][START_REF] Vuyst | The burden of human papillomavirus infections and related diseases in sub-saharan Africa[END_REF][START_REF] Chan | Metaanalysis on prevalence and attribution of human papillomavirus types 52 and 58 in cervical neoplasia worldwide[END_REF] In agreement, our multivariate analysis among the coinfected ICCs or only among the coinfected SCC cases demonstrated the significantly higher prevalence clustering of HPV18 in Africa and HPV58 in Asia. In contrast, HPV52 showed marginally significant higher geographic clustering among African women with intermediate prevalence in Asia while HPV16 or HPV45 coinfection prevalence values did not reveal any geographical structure. Most interestingly the observed patterns of geographic clustering of HPVs among the coinfected ICCs practically originate from differences between Africa and other world regions. As such, the prevalence of HPV51, HPV35, HPV18 and HPV52 types were about two-fold higher in African multiple infected ICCs, while prevalence of HPV56, HPV58 and HPV74 was either absent or at least two-fold lower in African ICCs compared to any other world region. Most importantly, the highest LRT values observed in our multivariate analysis were both the geographic clustering of multiple infections in general and at type level the geographic clustering of oncogenic HPV51, HPV35, HPV18 and HPV52 types in African multiple infected ICCs. Simultaneously, a recent community-randomized clinical trial data from Finnish population showed an increase of HPV51 and HPV52 in multiple HPVs infected women post-vaccinated. [START_REF] Gray | Evaluation of HPV type-replacement in unvaccinated and vaccinated adolescent femalespost-hoc analysis of a community-randomized clinical trial (II)[END_REF] However, within a public health perspective earlier studies have shown that although HPV51 is classified as possibly/probably oncogenic by the IARC, it presents a low oncogenic potential to cause ICCs and is commonly more prevalent in low-grade cervical lesions. [START_REF] Bzhalava | A systematic review of the prevalence of mucosal and cutaneous human papillomavirus types[END_REF] In addition, a minor clustering of HPVs by histological type was also observed among the HPV coinfected ICCs. A significantly higher contribution of HPV31 was observed in SCCs compared to other type of ICCs, mostly Our results suggest that this differential clustering of HPV type combinations in ICCs after stratifying by histological presentation and geographical origin might be compatible with a specific modulation of the natural history of the infection and/or of the oncogenic potential of particular HPVs multiple infections. Finally, although our results may be valuable in estimating competition between nonvaccine HPV types and warrant monitoring for post-vaccinated, it should be emphasized that clustering of multiple HPVs infections are marginal in ICCs with likely only a minor public health importance, and mostly among high sexual risk-taking behavior population groups.

2 Column percent. 3 Row

 3 percent of multiple infections among HPV/DNA-positive ICC cases. 4 Asia including Oceania. performed a multivariate analysis. A significant difference in HPVs composition by geographic origin [likelihood ratio test (LRT) 201.24, p = 0.001], by histology (LRT 58.604, p = 0.001) and by period of diagnosis (LRT 43.35, p = 0.008) was observed (Table

Figure 1 .5

 1 Figure 1. HPV type-specific prevalence among all single HPV type infected ICCs (n = 8193), and exclusively among multiple HPVs infected ICCs (N = 587). Abbreviations: HPV, Human papillomavirus; Taxonomy grouping based on ICTV (http://www.ictvonline.org) and visualized with gray color. [Color figure can be viewed at wileyonlinelibrary.com]

Figure 2 .

 2 Figure 2. Logistic regression analysis p values of the pairwise comparisons (a) among all HPV coinfections and (b) among cases only with two HPV type coinfections, and adjusted for age, period of diagnosis and geography. Abbreviations, OR, Odds ratio, lines indicate -log10(p value) = 2 and Bonferroni's corrected threshold for 2 HPV types co-occurrence -log10(0.01/146) = 4. [Color figure can be viewed at wileyonlinelibrary.com]

Figure 3 .

 3 Figure 3. Visualization of the coinfection patterns among all 587 multiple HPVs infected ICCs, and stratified by geographical origin of the samples. Nodes are numbered by the corresponding HPV type and node size is proportional to the observed HPV type prevalence in coinfected ICCs. Node colors represent the following: HPV types with significantly higher prevalence in coinfected ICCs in Africa compared to other world regions are shown in orange. HPV types with significantly lower prevalence in coinfected ICCs in Africa compared to other world regions are shown in blue. Each line represents the coinfection observed between the two HPV types with a width proportional to the total number of these coinfections observed in this study. Orange color of the coinfection line indicated that the two HPV types were more likely to be detected as coinfections and blue color indicated that the two HPV types were less likely to be detected as coinfections compared with any other HPVs coinfected ICCs. [Color figure can be viewed at wileyonlinelibrary.com]

Figure 4 .

 4 Figure 4. Classification and regression tree (CART) of the HPV16 coinfections. On top of the tree is shown the basal node containing the entire data set with the predicted probability of HPV16 coinfections, and also visualized in gray color. Two subsequent splits between HPV16 -HPV31, and between HPV16 -HPV45 suggest a likely competition between the types among the HPV coinfected ICCs. Number of cases (N) in each node and the corresponding predicted probabilities (PP) in percentages [95% confidence interval].

Table 1 .

 1 Characteristics of the ICC cases included in the study

		HPV positive 1	Multiple HPVs	Crude		Adjusted
	ICCs	N	(%) 2	N	(%) 3	OR	[95%CI]	OR	[95%CI]
	Age at diagnosis
	Mean (SD)	51.0	(13.17)	51.6	(13.8)	1.00	[0.996-1.01]
	Median (Min-Max)	50	(17-97)	51	(24-93)
	ANOVA p value			0.318
	LRT p value					0.315
	Age at diagnosis
	65+	1226	(14.0)	73	(6.0)	Ref.		Ref.
	55-64	1501	(17.1)	102	(6.8)	1.15	[0.84-1.57]	1.17	[0.86-1.60]
	45-54	2123	(24.2)	119	(5.6)	0.94	[0.69-1.27]	0.96	[0.71-1.30]
	35-44	1806	(20.6)	88	(4.9)	0.81	[0.59-1.11]	0.83	[0.61-1.15]
	25-34	635	(7.2)	40	(6.3)	1.06	[0.71-1.58]	1.14	[0.76-1.70]
	24	48	(0.5)	5	(10.4)	1.84	[0.71-4.78]	1.96	[0.75-5.13]
	Chi-squared test p value			0.162
	Trend test p value			0.365		0.365		0.435
	LRT p value					0.180		<0.001
	Missing	1441		160
	Geographic origin
	Africa	525	(6.0)	101	(19.2)	Ref.		Ref.
	Europe	2031	(23.1)	134	(6.6)	0.30	[0.22-0.39]	0.36	[0.27-0.49]
	Asia 4	2751	(31.3)	147	(5.3)	0.24	[0.18-0.31]	0.31	[0.23-0.42]
	America	3473	(39.6)	205	(5.9)	0.26	[0.20-0.34]	0.36	[0.27-0.49]
	Chi-squared test p value			<0.001
	LRT p value					<0.001		<0.001
	Period of diagnosis
	1940-1989	3513	(40.0)	202	(5.8)	Ref.		Ref.
	1990-2009	5267	(60.0)	385	(7.3)	1.29	[1.08-1.54]	1.32	[1.08-1.60]
	Chi-squared test p value			0.296
	LRT p value					0.004		0.005
	Histology
	Squamous cell carcinoma	8064	(91.8)	530	(6.6)	Ref.
	Other 5	716	(8.2)	57	(8.0)	1.22	[0.93-1.63]
	Chi-squared test p value			0.154
	LRT p value					0.165
	TOTAL	8780	(100.0)	587	(6.7)
	1 HPV positive by PCR with SPF-10 primers followed by DEIA and LIPA25 genotyping of 25 HPV types (HPV6, 11, 16, 18, 31, 33, 34, 35, 39, 40, 42,
	43, 44, 45, 51, 52, 53, 54, 56, 58, 66, 68, 70, and 74).

Abbreviations: ICC, Invasive cervical cancer; HPV, Human Papillomavirus; Crude, crude odds ratio (OR); adjusted, ORs adjusted by age at diagnosis, period of diagnosis and geographic origin; N, Number of HPV/DNA-positive cases, SD: standard deviation; LRT: log-likelihood ratio test.

Table 2 .

 2 Multivariate analysis of all 25 HPV types available among coinfected ICCs, and HPV type-specific prevalence distribution for eight HPV types showing significant clustering by histology, period of diagnosis or geographic origin of the samples

								Geographic origin
		Histology		Period of diagnosis	Geographic origin	Prevalence among HPV coinfected ICCs
		LRT	p value	LRT	p value	LRT	p value	Africa	Europe	Asia 1	America
	ICCs	58.60	0.001	43.35	0.008	201.24	0.001	(n = 101)	(n = 134)	(n = 147)	(n = 205)
	HPV51	0.02	1.000	1.6	0.976	30.81	0.001	29.7	6.0	6.1	8.3
	HPV35	3.60	0.673	0.76	0.998	22.17	0.003	18.8	10.5	2.0	5.9
	HPV18	7.40	0.087	2.2	0.943	22.15	0.003	39.6	20.2	21.1	12.2
	HPV74	0.35	1.000	12.10	0.010	18.99	0.006	1.0	4.5	5.4	13.2
	HPV58	0.08	1.000	1.21	0.998	15.48	0.026	3.0	6.0	17.0	8.3
	HPV52	0.25	1.000	3.69	0.671	12.61	0.065	30.7	12.7	18.4	13.2
	HPV56	0.22	1.000	0.28	1.000	12.50	0.065	0.0	8.2	5.4	5.4
	HPV31	17.08	0.001	0.15	1.000	8.70	0.347	5.9	11.9	15.0	18.0

Diagnosis of adenocarcinoma, adenosquamous carcinoma and to a lesser extent undifferentiated and neuroendocrine tumors.

Abbreviations: ICC, Invasive cervical cancer; Asia 1, Asia including Oceania; LRT: log-likelihood ratio test; Italic values indicate significant or marginally significant p values; Bolded values indicate the highest geographic prevalence values observed. HPV types in descending order regarding the LRT values by geographic origin.

Acknowledgements

We thank Jo Ellen Klaustermeier for excellent technical assistance. The study is part of the ICO SURVEY RIS HPVTT and has been partially supported by grants from the Instituto de Salud Carlos III-ISCIII (Spanish Government), cofunded by FEDER funds / European Regional Development Fund (ERDF) -a way to build Europe (References: RD12/0036/ 0056, CIBERESP, CIBERONC CB16/12/00401, PI030240, PI061246, PI17/ 00092, PERIS SLT002/16/00496) and from the Agència de Gestió d'Ajuts Universitaris i de Recerca (2014SGR1077, 2014SGR756, 2017SGR1718, 2017SGR1085), Marató de TV3 Foundation (051530), The European Society of Clinical Microbiology and Infectious Diseases (ESCMID Research Grant 2015), GSK, MSD, and Merck.

Research Program has received sponsorship and grants from Merck. Xavier F. Bosch has received occasional personal and travel grants from Merck, Hologic, Qiagen, MSD and Seegene. Grant sponsor: Merck; Grant sponsor: MSD; Grant sponsor: GSK; Grant sponsor: European Society of Clinical Microbiology and Infectious Diseases; Grant numbers: ESCMID Research Grant 2015; Grant sponsor: Marato de TV3 Foundation; Grant numbers: 051530; Grant sponsor: Agencia de Gestio d'Ajuts Universitaris i de Recerca; Grant numbers: 2014SGR1077, 2014SGR756, 2017SGR1718, 2017SGR1085; Grant sponsor: Instituto de Salud Carlos III (Government of Spain) and FEDER/ ERDF; Grant numbers: RD12/0036/0056, CIBERESP, CIBERONC CB16/12/00401, PI030240, PI061246, PI17/00092, PERIS SLT002/16/00496.

Author contributions

Conceived the project and wrote the study: VNP and ST. Performed the analysis: VNP, ST, YB. Contributed in data generation, study design and data interpretation: VNP, ST,YB, LA, NG, WQ, XFB, IGB, SS. Supervised the project: SS. All authors assisted in the writing of the study and had final approval of the submitted and published versions.