key: cord-0021090-8d2144gp
authors: Beneteau, Thomas; Selinger, Christian; Sofonea, Mircea T.; Alizon, Samuel
title: Episome partitioning and symmetric cell divisions: Quantifying the role of random events in the persistence of HPV infections
date: 2021-09-07
journal: PLoS Comput Biol
DOI: 10.1371/journal.pcbi.1009352
sha: c3c321ab43f112b5cc5536de720c1090950b2934
doc_id: 21090
cord_uid: 8d2144gp

Human Papillomaviruses (HPV) are one of the most prevalent sexually transmitted infections (STI) and the most oncogenic viruses known to humans. The vast majority of HPV infections clear in less than 3 years, but the underlying mechanisms, especially the involvement of the immune response, are still poorly known. Building on earlier work stressing the importance of randomness in the type of cell divisions in the clearance of HPV infection, we develop a stochastic mathematical model of HPV dynamics that combines the previous aspect with an explicit description of the intracellular level. We show that the random partitioning of virus episomes upon stem cell division and the occurrence of symmetric divisions dramatically affect viral persistence. These results call for more detailed within-host studies to better understand the relative importance of stochasticity and immunity in HPV infection clearance.

a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Human Papillomaviruses (HPV) are the most oncogenic viruses known to humans [1] . Up to 8.6% of all cancers occurring worldwide in women are caused by HPV infections, versus 0.8% for cancers in men [2] . Among the various types of cancers attributable to HPV, cervical cancers stand out since they are the third-most diagnosed form of cancer in women [3] .

The burden HPVs impose on human populations largely originates from their high prevalence in the general population. For example, it was estimated that the average lifetime probability (from sexual debut to age 70) of acquiring a genital HPV infection in the US was 84.6% for women and 91.3% for men [4] . But since up to 90% of these infections remain subclinical and are cleared in less than 3 years [5] , viral persistence, in particular for specific genotypes such as HPV16, determines potential oncogenicity [6] .

Despite its ubiquity, the mechanisms underlying HPV infection clearance and persistence are still poorly understood. Although the innate immune response (e.g. natural killer cells) is known to be involved in chronic infections [7] , few data exist for non-persistent infections. Adaptive immune responses (e.g. the role of T helper cells) have been studied in the context of vaccine trials, but the results remain inconsistent [8] . It has been noticed that people with recurrent genital warts were more likely to be HPV seropositive than those who do not have a history of genital warts [9] . However, more generally, the breadth and magnitude of adaptive immune responses following HPV infection are believed to be limited [10] due to innate immune evasion mechanisms and a lack of inflammatory responses elicited by infected epithelial cells. Furthermore, in natural infections, there is no viremia (i.e. presence of viruses in the blood), and the production of detectable antibodies against the virus is low [10] . Another longitudinal study investigated HPV16 seroconversion statuses among various groups of infected women but did not reach any conclusive result, potentially by lack of statistical power [11] .

The inherently stochastic nature of the episomal life cycle of HPV, i.e. replication of the viral genome as independent circular DNA copies, could be an additional factor contributing to the clearance of HPV infections. HPVs are non-lytic viruses that persist in the basal stem cells by replicating with host cell chromosomes [12] . In the basal layer, most cell divisions are asymmetric: one basal stem cell gives birth to one daughter stem cell, and one daughter differentiated cell upon cell division. Occasionally, especially in the context of a perturbation such as wounding [13, 14] , symmetric cell divisions occur that either yield two stem cells or two differentiated cells. Infection persistence in the epithelial tissue occurs when viral copies are transmitted to the daughter stem cell during cell divisions. Conversely, transmission to new hosts occurs via the viral copies in the daughter differentiated cell. These copies migrate and maturate with the cell before being released to the surface upon shedding.

Symmetric divisions of infected basal stem cells into two daughter stem cells or two daughters differentiated cells [13, 15] strongly impact the spread and persistence of the virus in the basal layer of the tissue. Since HPV infections begin with few infected cells [16] , randomness in the type of cell division plays a crucial role in the clearance of HPV infections [17] . A second source of stochasticity occurs within the nucleus of infected cells. During stem cell divisions, the number of episomes that are passed on to each daughter cell is random. Since cells are infected with a small number (potentially unique) of viral copies [18] , stochasticity in episome distribution upon division could impact the persistence of HPV infection. The partition of viral plasmids in human cells is uncommon but not specific to HPV. This phenomenon is still not clearly understood despite its importance in the maintenance of infections and progression towards cancers [19] . To date, this aspect of HPV infections has not been taken into account. In a recent review [20] on the use of mechanistic mathematical models applied to HPV infection, the authors globally call for more within-host modeling studies. While healtheconomic models have been extensively used to evaluate the cost-benefit of public interventions in the prevention of HPV infection or cervical cancer, mechanistic models have mainly been used to study some specific topics (e.g. viral evolution induced by vaccine use [21, 22] or impact of cell division in HPV natural infection [17] ).

Building on the work of Ryser [17] that emphasized the contribution of symmetric cell divisions on HPV clearance, we develop an original model of HPV viral copies dynamic in epithelial basal layers.

Taken separately, the intra-and inter-cellular stochastic processes fail to capture part of the life cycle of HPV infections. The intracellular approach does not account for the spread of the infection inside the epithelial tissue, whereas the inter-cellular approach ignores the necessity of viral persistence in proliferating cells. Thus, we hypothesize that combining the two aspects can provide additional insights into the dynamics of HPV infection.

We investigate the role of stochasticity in the clearance of HPV infections and quantify the role of both symmetric divisions and episomal partitioning in this process. We provide further evidence for the crucial effect of randomness in the clearance of HPV infections.

Our model focuses on the basal layer of a stratified squamous epithelium (SSE), but it can easily be adapted to other tissues as the HPV replication pattern is essentially the same in most sites [23, 24] . We only follow the dynamics of episomes in infected stem cells because SSE is in perpetual renewal and the turnover interval of the non-dividing cells lasts 3 to 6 weeks [25] [26] [27] . Therefore, in the absence of cell-to-cell transmission, it is improbable that an infection could persist only in the differentiated layers of SSE. We also neglect potential reinfection of the SSE from virions released by desquamating cells because these are thought to be unlikely given the importance of micro-wounding in the establishment of new infections [12] .

The life-cycle is divided into two steps following the non-lytic properties of HPV: 1) during the division of an infected cell, HPV episomes are distributed randomly between the two daughter cells, and 2) infected stem cells can divide symmetrically, thereby ending or, conversely, further spreading the infection. We graphically summarize each aspect of the model in the S1 Fig. Both the intra-and inter-cellular mechanisms are modelled as discrete-time branching processes [28, 29] .

More formally, let X n ; n 2 N be the total number of HPV episomes in all infected stem cells, n being the number of stem cell generations since the beginning of the infection. Unless stated otherwise, a time step corresponds to one cell generation. X n is a random variable, the distribution of which is non-trivial, for it is constructed following several random events. We assume the random events between cell lineages, i.e. descending cells by asymmetric divisions, to be independent and identically distributed. Within a cell lineage, we assume other stochastic events to be identically distributed (but not independent). Independence is verified under certain specific conditions, for instance if we assume a cell lineage only divides asymmetrically.

HPVs are non-lytic viruses that follow the life cycle of their host cells. In basal (stem) cells, the number of viral copies is usually considered to be limited to several dozens of episomes per cell: 10 to 200 according to Doorbar [16] , and 50 to 100 according to more recent work [12] . We denote by C the maximum number of episomes per infected stem cell. We model intracellular dynamics as the superposition of two distinct random steps:

1. A division, where each episome in the dividing cell is distributed to the two daughter cells according to a Bernoulli distribution Bernoulli(p). For asymmetric divisions, parameter p stands for the probability for an episome to be relocated in the daughter stem cell. Episomes distributed to the differentiated daughter cell are considered withdrawn from the system as they migrate to the surface without contributing to the infection of stem cells.

2. An amplification, which occurs once the division over. For this random phenomenon we either assume a Dirac distribution δ(λ) or a Poisson distribution Poisson(λ), with l 2 N. We refer to λ as the episome replication factor.

Without prior knowledge of the biology of the amplification phase, the Dirac distribution appears to be the most parsimonious as it implicitly assumes that all viral copies exhibit the same average behavior. However, this assumption could be biased at the beginning of the infection as viruses remain in small numbers. In this case, the number of copies amplified could diverge from the average behavior due to individual variance. As little information is available in the literature, we choose a Poisson distribution as it only requires the same parameter as for the Dirac case.

The order of the two steps can be reversed without affecting the intrinsic characteristics of this branching process (see S1 Text for more details). Since HPV genomes are present in low copy numbers and HPV gene expression is set at a minimum level in basal stem cells [24] , we assume that the number of HPV viral copies does not affect the dynamics of the host stem cells.

The infection can spread in the basal layer following symmetric divisions in two stem cells, a phenomenon governed by strong stochastic effects and detailed in the next section. To distinguish between daughter stem cells, we use the notation described in the supporting information S1 Text and define p as the probability for an episome to be distributed in the left (L) daughter stem cell. We assume episome distribution bias towards one of the daughter cells does not apply when the two daughter cells are stem cells, hence we assume distribution to be even (p = 0.5) during symmetric division.

In what follows, we group all the layers apical to the basal layer in a general population of differentiated cells. Let ðG n Þ n2N be the number of infected stem cells at time n. We denote by s resp. r the probability of a symmetric division to yield two stem cells resp. two differentiated cells. We assume s and r to be small compared to the probability of asymmetric cell division. Following the work of Ryser [17] , we model infected stem cells dynamics using the following branching process:

where ϕ i are independent and identically distributed random variables with the following distribution:

where ψ is a random variable accounting for intracellular effects leading to the end of the infection of the infected stem cell S. A histological study from 1970 estimated that r + s ⩽ 0.08 [25] , and a more recent analysis also estimated 8% of cell divisions to be symmetric [15] .

Some analytical insights can be obtained from the previous equations but in-depth investigations require numerical simulations. We estimated the cumulative probability of extinction (p ext ) given a set of parameters (see Table 1 for further details on the parameters used) by carrying out O = 10 3 independent runs parametrized by the same parameter values. We also include the type of intracellular amplification in the analysis (Dirac or Poisson distributed). We displayed a random trajectory of episomal dynamic in the S2 Fig. Using the O different trajectories, we calculate the cumulative probability of extinction p ext (t) as follows:

where X k t jS is the total number of episomes at time t given parameters set S for the k-th trajectories. If Γ 0 > 1, we assume that all infected cells initially harbor N 0 viral copies, hence a total of X 0 = Γ 0 × N 0 episomes. p ext is computed for all t ⩽ 100 cell generations. Based on data from the 1970s [25] , we bound the rate of divisions between [0.03, 0.07] day -1 , which means the upper bound of 100 cell divisions is sufficient to estimate infection clearance after 3 years. For each t 2 N, we estimate p ext using Maximum Likelihood estimation for a Binomial proportion. For each of these proportion we computed the 95% confidence interval using Agresti-Coull method [31] .

To evaluate the role of each parameter on the probability of extinction, we perform a global sensitivity analysis using variance-based methods (also known as Sobol indices). This approach decomposes the variance of the simulation outputs (here p ext ) into fractions that can Table 1 . Parameters notation and description along with their biologically estimated ranges. The parameters p, λ, C and N 0 govern the intracellular process, the three others are involved in the inter-cellular approach. Another parameter (not numerical) can be added to the list: the type of intracellular amplification (Dirac or Poisson distributed) as the two types slightly differ in their construction. Extensive analysis of the impact of each distribution on the probability of extinction is shown in S1 Text. be attributed to each of our parameters or interactions between them [32, 33] . To do so, we store parameters sets in two different matrices A and B of dimension n × 7, where the columns correspond to the six parameters p, λ, s, r, C, and N 0 (see Table 1 ) and the type of intracellular amplification. We denote by A ðkÞ B the matrix A whose k-th column has been replaced by the kth column of matrix B. We generate A and B sets using the Latin Hypercube Sampling (LHS) method implemented in the lhs package [34] in R v4.0. 3 [35] . Each matrix contains n = 500 different parameter sets (more details can be found in S1 Text). We execute the global sensitivity analysis using multisensi package [36]. We compute model outputs (i.e. extinction probabilities) for all sets A, B, A ðkÞ B ; k 2〚1; 7〛to obtain the following Monte-Carlo estimator for the variance of the extinction probability up to time t [37]:

We do not include the initial number of infected stem cells (Γ 0 ) in this analysis, since the sensitivity of p ext with respect to Γ 0 follows from the branching process property:

Our model contains two nested branching processes, which makes it difficult to obtain analytical results. However, analytical expressions can be derived if we neglect either of the two stochastic processes (see S1 Text for detailed calculus and S3 and S4 Figs for numerical verification) and we use these results as boundaries to interpret the results of our numerical simulations and to determine how the cumulative probability of extinction (p ext ) varies for the different parameters. Simulations also allow us to explicitize this direction of variation in p ext with each parameter. Fig 1 displays the difference between probabilities of extinction for two parameter sets that only differ by the value of one parameter, we subtract the probability of extinction for the lowest of the two varying parameters to the probability of extinction of the highest of the two.

If stochasticity only acts at the within-cell level, i.e. symmetric divisions are assumed not to occur, we retrieve classical results from Bienaymé-Galton-Watson (BGW) branching process theory [29] : the extinction probability p ext then decreases with the episome replication factor, λ, and with the probability for an episome to end up in the daughter stem cell, p (see Fig 2 for the Poisson case and S5 Fig for the Dirac distribution) . As λ increases, we observe that p ext converges toward 1 − p. Biologically, this means that the main event governing the probability of extinction is the first cell division and the distribution of the original episome into a differentiated or a stem cell.

Numerical simulations allow us to investigate the effect of intracellular processes while accounting for symmetric divisions. In particular, our sensitivity analysis highlights the effect of λ and p (Fig 3) . At the beginning of the infection, intracellular phenomena explain around 70% of the variance, mostly through the probability of allocation to daughter stem cell p. Indeed, in the early phase of the infection, the total number of episomes is usually low and, since symmetric divisions are rare (r + s ⩽ 0.08 [25] ), the early dynamic of episomes is mostly governed by intracellular processes. We also retrieve the results from the analytical model, although the decrease in p ext with increasing λ is more limited (Fig 1) .

As explained in the Methods, the numerical simulations slightly deviate from the analytical framework because the number of episomes per cell is assumed to be limited to C. However, we expect most of these results to hold in the general framework because we focus on a maximum time scale of 100 cell generations. We observe more deviations on that time interval when C is small and when the intracellular regime is slightly supercritical (further details can be found in S1 Text and S5 Fig) .

When omitting intracellular stochasticity, theory allows us to predict that the branching process will go extinct with probability p ext = min(r/s, 1), which we confirm with numerical simulations, (S4 Fig). Besides, we also find consistent results with theory as the estimated p ext decreases (resp. increases) with s (resp. r), in Fig 1. As intracellular stochastic processes only increase the risk for the infection to go extinct, we can use as the lower bound for the cumulative probability of extinction of our general process the cumulative probability of extinction for the sole inter-cellular process, the proof is detailed in S1 Text. Generally speaking, no matter how rapidly the virus replicates at the intracellular On each panel, the title indicates which parameter varies between the two sets. In each of this panel, we subtract the probability of extinction for the lowest of the two varying parameters to the probability of extinction of the highest of the two. We find consistent results as theory prediction: the probability of extinction decreases with λ, p, s, N 0 and increase with r. The viral capacity in cell (C), has little to no effect on the probability of infection. Line colors show the median (in red) and the first and third quartiles (in blue and green).

https://doi.org/10.1371/journal.pcbi.1009352.g001 level, its spread is bounded by the cellular dynamics of the host tissue. Therefore, we have

Globally, these results suggest that symmetric divisions can favor the persistence of the virus depending on the local context of the tissue. In particular, when s > r, e.g. when the epithelium switches from a maintenance behavior to a repair behavior [13, 14] , symmetric divisions appear to favor the spread of the virus, i.e. its persistence in the tissue. Additionally we can notice in Eq 2 that if the inter-cellular process is critical (i.e. r = s), then the general process is subcritical. Indeed, under such circumstances we have

In Fig 3, the proportion of the variance explained by symmetric divisions increases over time. Indeed, on average the first symmetric division occurs after (1/(r + s) cell generations. Hence, the early phase of the infection is mostly governed by intracellular aspect. But later in the infection, if the time between two symmetric events is sufficiently large, the intracellular dynamic has almost surely converged to one of the two states {0, C}. Under such circumstances, the inter-cellular aspect is the only one acting on the dynamic of the infection. We Effect of episome replication (λ) and asymmetry in episome partitioning (p) on the cumulative probability of extinction (p ext ). We here assume a Poisson scenario. p ext decreases with both λ and p, and reaches a threshold p ext = 1 − p for a given p when λ increases. This is consistent with the fact that if λ is sufficiently high, the main source of extinction is the first cell division of a stem cell containing 1 episome.

https://doi.org/10.1371/journal.pcbi.1009352.g002 observe a significant difference in the proportion of variance explained by the parameter r compared with s, the latter having a smaller slope. This difference is directly induced by the intracellular stochasticity that can terminate infections prematurely.

We then assess the contribution of the initial number of episomes and their distribution in one or more stem cells on the probability of extinction. Let us assume the branching process to start with a total of X 0 = ν × ω episomes, ðn; oÞ 2 N 2 ; n > o. We then look for the optimal distribution of these episomes from the virus' perspective. For this, we evaluate the probability of extinction under two extreme scenarios: i) ν cells each containing ω episomes, or ii) ω cells each containing ν viral copies. From branching theory, we know that for the inter-cellular level, the probability of extinction when starting with a 2 N identical and independent infected cells is equal to the probability of extinction of one cell raised to the power α.

Therefore, we only need to compute the probabilities of extinction for two sets of parameters that are identical, except for the number of episomes and number of infected stem cells at time t = 0, that are switched. We use numerical simulations to evaluate those probabilities, To determine which scenario favors virus persistence we compare the cumulative probabilities of extinction for each of pairs of parameters that only differ by the initial arrangement of the viral copies. Fig 4 shows the quartiles for the two extreme scenarios. This indicates a benefit for the virus to spread its copies to the maximum number of stem cells at the start of the infection, the bigger the difference between ν and ω the more significant this benefit is.

As immune mechanisms fail to completely explain the problem of spontaneous clearance in natural infections, Ryser [17] has suggested random events in the type of cell divisions could induce extinction and thus complement explanations based on immunity. Following this approach, we refine the model initially exposed by Ryser [17] and developed a model that also accounts for the risk of stochastic extinction due to virus partitioning. We hypothesized that an HPV infection could face extinction due to the sole effects of random processes occurring both at the intra-and inter-cellular level. Using a stochastic model comprised of two nested branching processes and numerical simulations, we find that random allocation of episomes during cell divisions drastically limit the persistence in the early phase of the epidemic. Later, if episomes persisted in some infected cells, randomness in the type of cell divisions can also lead to extinction, even though depending on the local context, symmetric divisions can also favor the spread of the virus in the tissue (when s > r). Notably, it was shown that the same cells can exhibit different proliferative patterns, switching from a balanced mode to an expanding mode if the tissue is wounded, and back to balanced mode once confluence is reached back [13, 14] . While the precise molecular mechanisms underlying those shifts are still unknown, studies suggest that actin signaling could play a key role in the fate of the keratinocytes [14] . The importance of each aspect varies throughout the infection, and more generally as HPV viral copies increase in number, the role of stochasticity in the explanation of extinction fades. Our results show the first stages of the infection matter a lot for the persistence of the virus, which is particularly important for viruses like HPV that have a low number of copies initially [18] .

Another outcome of our model is that many infections by HPV might clear naturally in a matter of weeks/months before even being noticed by the host immune system since episomes remain in small number in the basal layer. Adding to the pioneering work of Ryser [17] to highlight the role of stochastic processes in the clearance of HPV infections, we show that intracellular stochasticity greatly impacts extinction probability, especially in the early stages of the infection. Combining the two processes enables us to reach substantial clearance rates without involving immune response.

These results could explain why the antibody production against the virus is low after natural infections [10, 18] , as the infection might fade randomly within the first divisions of the initially infected cells. HPVs might have evolve to limit pressures exerted by randomness during the first cell divisions by amplifying the viral copies in the early phase of the epidemic. Notably in the case of HPV16, the E1 protein seems to play a decisive role in the establishment of a persistent infection until the number of episomes per cell reaches host capacity but the same protein was dispensable for the maintenance of the infection. Our model corroborates these observations: we show that the early phase of the infection can be detrimental to the infection, E1 protein could favor the maintenance of the infection either by acting on the probability to be distributed to the daughter stem cell (p parameter) or promote episome amplification (λ parameter). Upon reaching host cell capacity, even if the intracellular regime is critical, it is very unlikely for the virus to go extinct, hence E1 protein is dispensable during maintenance phase of the infection [38] .

Papillomaviruses have evolved over million years infecting various animals' epithelium [39] . From this long coevolution between the virus and its host, a large diversity of evolutionary strategies have emerged. For instance, HPV genotypes have developed specific host-and tissue-tropism that manifest through distinct pathologies [40] [41] [42] . HPV research mostly focuses on the Alpha genus since it contains all high-risk HPVs (HR-HPVs) that are susceptible to induce various forms of cancer in the genital areas and upper aerodigestive area [2, [43] [44] [45] . The Alpha genus also contains genotypes that cause benign infections and are classified as low-risk HPVs (LR-HPVs). Although benign, LR-HPV infections can cause hyperproliferative lesions such as genital warts. These HPVs have been less studied, compared to HR-HPV, despite accounting for a bigger number of types [46] . The distinct evolution of LR and HR-HPVs, which illustrates the persistence-virion production trade-off, has been investigated in the literature [47, 48] . The role of the oncogenic proteins E6 and E7 have been put forward to explain the capacity of the HR-HPV to induce cellular immortalization and cellular progression towards cancer development [49, 50] . However, it is unclear how these differences could impact the early phases of the infection and interfere with the stochastic pressures acting on the virus.

Globally our results suggest that remaining with a low copy number is a way to evade immunity [51] , but such a strategy is risky for the virus and is subject to strong stochastic forces that limit the persistence of the infection. Viruses committed to such strategy might alleviate these pressures exerted by randomness using proteins in the very early phase of the infection, as the E1 protein for HPV16 [38] . Unfortunately, these proteins might not be the most adequate targets for antiviral drugs as studies showed that the role of these are crucial in the very early steps of the infection but lose importance once cell carrying capacity is reached [38] . Therefore, we might not be able to detect HPV sufficiently early to inhibit E1 protein to prevent infection persistence.

Ryser [17] assumed a balance between the two types of symmetric divisions, thus restraining the inter-cellular dynamic to the critical regime. We allow the possibility for the tissue to deviate from that assumptions. Notably, we allow the probability of symmetric divisions in two stem cells to exceed the probability of symmetric divisions in two differentiated cells as keratinocytes can alternate between two modes of proliferation depending on local context. If local confluence is removed by scratch injury, cells switch to the expanding mode, characterized by an excess of cycling cells production, until confluence is reached again, after which they switch back to a balanced mode where similar proportions of proliferating and differentiating cells are produced, thus generating population asymmetry maintaining epidermal homeostasis [14] . The two modes of proliferation could differ by the rates of cell divisions or the probability of division in two stem cells (parameter s), or both. Our model predicts that HPVs greatly benefit from the expanding mode, which echoes previous empirical work on the link between HPV persistence and wound healing responses [16, [52] [53] [54] .

To further calibrate the models and delineate the true role of stochasticity, a comparison with longitudinal follow-ups of women recently infected by any type of HPV or presenting chronic latent infections would be usefull. Ideal time windows between two measurements should not exceed 1 month. Besides, we should perhaps reconsider how we label infections that cleared in less than 1 month. Many follow-ups treat them as "carriage", i.e. HPV from the partner and not a "true" infection but perhaps this underestimates the number of non-persistent infections. Overall to decipher this issue, we need accurate infection duration distributions over short time scales, ideally with virus loads data.

Comparing longitudinal follow-ups data and numerical simulations would require clarifying the estimation of the time between two cell divisions. We rely on old measures [25] to set the upper bound of the number of cell generations to compute in our three years simulations. This work estimates the rate of divisions around [0.03, 0.07] day -1 . These results are consistent with more recent work on the kinetics of epidermal maintenance stating an overall division rate of 1.1 per week [55] . To fit his model to follow-ups data, Ryser [17] estimated the rate of division using proxies (height of the stratified squamous epithelium expressed in number of cells, renewal time of the tissue, fraction of proliferative cells in the basal layer) and found a much higher division rate [0.29, 1] day -1 . Improving our knowledge of withinhost processes would be of primary interest to compare theoretical and experimental approaches. with finite intra-host capacity (plain dots, ribbons correspond to 95% confidence intervals) and according to theory (dashed lines). In most cases, the estimations of p ext under constrained host capacity do not diverge from theoretical predictions on the time window considered. Deviations emerge when the intra-host regime is slightly supercritical and/or the intra-host capacity is limited (<50). The colors indicate the row of the parameter matrix LHS_intra (see S1 Text) used to obtain the cumulative probability of extinction. (TIF)

We observe no significant differences between estimations and theoretical predictions. The colors indicate the row of the parameter matrix LHS_inter (see S1 Text) used to obtain the cumulative probability of extinction. (TIF) S5 Fig. Effect of episome replication (λ) and asymmetry in episome partitioning (p) on the cumulative probability of extinction (p ext ). We here assume a Dirac scenario. p ext decreases with both λ and p, and reaches a threshold p ext = 1 − p for a given p when λ increases. This is consistent with the fact that if λ is sufficiently high, the main source of extinction is the first cell division of a stem cell containing 1 episome. (TIF)

We display in red the estimations of p ext in the scenario where the initial few viral copies are spread in more cells and in blue the scenario where more viral copies are spread in less cells. The solid lines indicates the estimations of p ext and the ribbon the 95% confidence intervals. The cumulative probability of extinction is generally lower in the first scenario compared to the latter. Thus it is more beneficial for the virus to spread its copies in the maximum number of cells at the the start of the infection. On each facet, the title indicates the row of the parameter matrix LHS_A (see S1 Text). (TIF) S1 Text. Mathematical details and extra information on the simulation procedure. (ZIP) 

Global burden of cancers attributable to infections in 2012: a synthetic analysis. The Lancet Global Health

Worldwide burden of cancer attributable to HPV by site, country and HPV type

Cervical-Cancer Screening with Human Papillomavirus and Cytologic Cotesting

The Estimated Lifetime Probability of Acquiring Human Papillomavirus in the United States

Human Papillomavirus Testing in the Prevention of Cervical Cancer

Incidence and Duration of Cervical Human Papillomavirus 6, 11, 16, and 18 Infections in Young Women: An Evaluation from Multiple Analytic Perspectives

Origin and immunoescape of uterine cervical cancer. La Presse Mé dicale

Cell-Mediated Immune Response to Human Papillomavirus Infection

Use of Human Papillomavirus Type 6 Capsids to Detect Antibodies in People with Genital Warts

Immunobiology of HPV and HPV vaccines

Cell-Mediated Immune Responses to Human Papillomavirus 16 E6 and E7 Antigens as Measured by Interferon Gamma Enzyme-Linked Immunospot in Women With Cleared or Persistent Human Papillomavirus Infection

Mechanisms by which HPV Induces a Replication Competent Environment in Differentiating Keratinocytes

A Single Progenitor Population Switches Behavior to Maintain and Repair Esophageal Epithelium

Human keratinocytes have two interconvertible modes of proliferation

A single type of progenitor cell maintains normal epidermis

Molecular biology of human papillomavirus infection and cervical cancer

HPV clearance and the neglected role of stochasticity

The Biology and Life-Cycle of Human Papillomaviruses

Plasmid Partitioning by Human Tumor Viruses

Mechanistic mathematical models: An underused platform for HPV research

Revising ecological assumptions about Human papillomavirus interactions and type replacement

Could the human papillomavirus vaccines drive virulence evolution?

Biology of Human Papillomavirus Infections in Head and Neck Carcinogenesis

Epithelial Cell Responses to Infection with Human Papillomavirus

Autoradiographic analysis of cell proliferation kinetics in human genital tissues

Immune responses to human papillomavirus

Pathology and epidemiology of HPV infection in females

Branching Processes. The Annals of Mathematical Statistics

Branching Processes

Interaction of Bovine Papillomavirus E2 Protein with Brd4 Stabilizes Its Association with Chromatin

Approximate Is Better than "Exact" for Interval Estimation of Binomial Proportions. The American Statistician

Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates

The E1 Protein of Human Papillomavirus Type 16 Is Dispensable for Maintenance Replication of the Viral Genome

The clinical importance of understanding the evolution of papillomaviruses

Cutaneous and mucosal human papillomaviruses differ in net surface charge, potential impact on tropism

Epithelial Tropisms, and the Development of Neoplasia

Human papillomavirus molecular biology and disease association

Evidence for a Causal Association Between Human Papillomavirus and a Subset of Head and Neck Cancers

HPV: The global burden

Epidemiology of HPV-associated oropharyngeal cancer

The low-risk papillomaviruses

The carcinogenicity of human papillomavirus types reflects viral evolution

Evolution, Medicine, and Public Health

Functional roles of E6 and E7 oncoproteins in HPV-induced malignancies at diverse anatomical sites

HPV16 E7 Genetic Conservation Is Critical to Carcinogenesis

Latent papillomavirus infections and their regulation. Current Opinion in Virology

Molecular genetics of human cervical cancer: role of papillomavirus and the apoptotic cascade

Building Epithelial Tissues from Skin Stem Cells

Current understanding of the mechanism of HPV infection

Kinetics of cell division in epidermal maintenance

The authors acknowledge the IRD itrop HPC (South Green Platform) at IRD Montpellier for providing HPC resources that have contributed to the research results reported within this paper.