Abstract

To date, little attention has been paid to the potential role of stochastic cell dynamics in the time it takes to clear an HPV infection. In this study, we combine mechanistic mathematical models at the cellular level with epidemiological data at the population level to disentangle the respective roles of immune capacity and cell dynamics in the clearing mechanism. Our results suggest that chance—in form of the stochastic dynamics of basal stem cells—plays a critical role in the elimination of HPV-infected cell clones. In particular, we find that in immunocompetent adolescents with cervical HPV infections, the immune response may contribute less than 20% to virus clearance—the rest is taken care of by the stochastic proliferation dynamics in the basal layer. In HIV-negative individuals, the contribution of the immune response may be negligible.

Author Summary

The most common cancer types attributed to HPV are cervical and anal cancers, but HPV-related head and neck cancers are on the rise, too. Even though the lifetime risk of infection with HPV is as high as 80%, most infections clear spontaneously within 1—2 years, and only a small fraction progress to cancer. In order to identify who is at risk for HPV-related cancer, a better understanding of the underlying biology is of great importance. While it is generally accepted that the immune system plays a key role in HPV clearance, we investigate here a mechanism which could be equally important: the stochastic division dynamics of stem cells in the infected tissues. Combining mechanistic mathematical models at the cell-level with population-level data, we disentangle the contributions from immune system and cellular dynamics in the clearance process. We find that cellular stochasticity may play an even more important role than the immune system. Our findings shed new light onto open questions in HPV immunobiology, and may influence the way we vaccinate and screen individuals at risk of HPV-related cancers.

Introduction

Over 90% of cervical cancers are caused by HPV infections, and up to 60% of squamous cell carcinomas of the vulva, vagina, anus and penis are associated with high-risk types of HPV [1]. More recently, it has been shown that infection with HPV also plays a critical role in the genesis of certain head and neck cancers, particularly in cancers of oropharynX and base of tongue [2]. The incidence of these cancers in men has been increasing over the past decade, suggesting the emergence of a virus-related cancer epidemic [3]. Even though the lifetime risk of HPV infections is as high as 80% [4], most individuals clear the virus within 1—2 years [5]. However, if infection with a high-risk type of HPV persists, the viral genes can interfere with the cellular control mechanisms and trigger neoplastic changes, which can eventually develop into an invasive carcinoma [6].

In particular, the mechanisms of virus clearance are controversial [8]. Clearance of HPV infection is usually attributed to an effective immune response, and the observation of longer clearance times in immunocompromised individuals further corroborates this assumption [9]. On the other hand, the fact that development of antibodies preventing future reinfection after clearing of the virus (known as seroconversion) occurs only partially [10—14] suggests that mechanisms other than an effective immune response may contribute to viral clearance.

Across different organs (both anogenital and oropharyngeal), oncogenic types of HPV preferentially infect areas of stratified squamous epithelium (SSE), and these SSE are not just a static backdrop to the unfolding infection process [2, 15]. They have a relatively fast turnover rate and the entire thickness of the epithelium is renewed every few weeks. During the renewal process, stem cell-like progenitor cells (hereafter denoted as 8 cells) in the lowest layer of the tissue (the basal layer) produce commited daughter cells (denoted as D cells) that differentiate and move upwards into the intermediate and superficial layers, and eventually get sloughed off into the lumen [15], see Fig. 1. The critical role of the dynamic differentiation and maturation process in the viral life cycle is well established [16]. However, the hypothesis that stochastic dynamics in the basal layer could contribute significantly to the clearing of new infections has not been addressed elsewhere.

Several mouse studies have used fluorescent labeling to observe lineage dynamics over time, and have concluded that while 8 cell division is prominently asymmetric (yielding one 8 and one D cell), a small fraction of 8 cell divisions are symmetric, yielding either two stem cells or two differentiated daughter cells [18, 19]. Considering that HPV infections start with a small number of infected basal cells in the SSE [16] , it seems plausible that these stochastic division patterns in basal cells may have an impact on the persistence properties of the infection. To investigate the relevance of cellular proliferation patterns and tissue homeostasis on HPV infection dynamics, we develop in this study a stochastic model of HPV infection in the SSE. By explicitly accounting for the stochasticity in stem cell proliferation, as well as cytotoxic T-cell mediated elimination of infected basal cells, we investigate the potential role of chance in the viral clearing process. Combining the model with a longitudinal data set of cervical HPV infections, we provide evidence for the critical role of stochasticity in HPV clearance.

Methods

Model

There is a large overlap among HPV types found in lesions of different sites, and HPV-16 is the most common type found in all HPV-related invasive cancers [20]. In addition, the viral replication strategy is essentially the same across affected sites [21, 22]. On the other hand, there are some organ-specific differences with respect to the biology of the affected stratified squamous epithelia. In fact, cervical, anal and oropharyngeal infections are usually restricted to a confined metaplastic transformation zone that separates columnar and squamous regions of the epithelium, whereas infections of e.g. the vulva, vagina and penis do not take place in such a transformation zone [23, 24]. Nevertheless, the bottom-up renewal dynamics (as explained below) of the affected epithelia are very similar, and the parametric model developed here can be applied to different tissue types by virtue of adjusting the relevant parameters, such as density of stem cells in the basal layer and regeneration time of the epithelium.

In physiological equilibrium, SSE (Fig. 1A) are dynamic tissues, and their entire thickness is renewed every few weeks. The renewal time AT is site-specific and varies between 3—6 weeks in the cervix [22, 25] and 2—3 weeks in the oral mucosa [26]. The regeneration is a bottom-up process as depicted in Fig. 1B, and the following dynamic model has been proposed based on lineage-tracing experiments [18, 19] (Fig. 1B). Stem cell-like progenitor cells (8) in the basal layer divide at a rate A, and during each division, there are three possible outcomes: with probability 1 — 2r (where 1’ << 1), 8 cells divide asymmetrically into one 8 and one differentiated daughter cell (D); with probability 2r, the division is symmetric, resulting in either two 8 or two D cells. Since space in the basal layer is limited, D cells leave the basal layer at rate F, and move upwards into the intermediate layers. Once fully matured and differentiated, D cells have reached the superficial layers, where they are eventually shed off to make space for new cells. Mathematically, these dynamics are summarized as a two-type branching process, see also Fig. 1B, where (Z) signifies cell death [18, 27]. In homeostasis, the basal layer consists of a conserved fraction p of S cells, which means that proliferation and migration rates satisfy the relationship pi = F(1 — p). Since the 8 cells in (1) undergo a critical branching process, their progeny will eventually go extinct (see also the discussion of 8*< cells below). However, since the basal layer is not compartmentalized and there is a large pool of 8 cells, this is very unlikely to occur within a human life time. In addition, it has been shown that there are regulatory mechanisms for stem cell fate [19, 28], and it is conceivable that similar mechanisms prevent the total number of 8 cells to fluctuate significantly. Finally, even though this model was originally developed based on mouse experiments, it has since been corroborated in human SSE [29].

New HPV infections arise in the basal layer of the SSE: physical ruptures in the tissue allow virions to reach the bottom layer of the SSE, where they infect residing 8 and D cells [16]. Due to the high turnover rate of SSE, non-dividing and upward moving D cells are lost from the epithelium within a few weeks, and hence persistent infections require the infection of 8 cells. Hereafter, we denote infected 8 and D cells by 8* and D*, respectively. Furthermore, we denote by nX(t) the number of cells of type X E {8, D, 8*, D*} present at time t. Since the viral count is kept at very low copy numbers (10—100) in the basal layer, and there is only minimal viral gene expression [30], we assume that the host cell dynamics are not affected by the presence of the virus. In particular, it has been shown that HPV-infected cells only acquire a selective growth advantage once the viral DNA has been integrated into the host DNA [31], which occurs at later, symptomatic stages of the infection. Consequently, the dynamics of infected cells in the early stages are still governed by (1), with 8 and D replaced by 8* and D*, respectively. From this it is easy to see that the dynamics of the 8* cell population in the basal layer are governed by the continuous-time critical branching process It follows from the theory of branching processes [32] that the probability an infected clone (starting from one infected cell) survives until time t is In particular, the clone will die out With probability 1, though the time until this occurs has infinite expected value.

Even though HPV is equipped with molecular mechanisms that facilitate immune evasion after infection, it is generally assumed that clearance of the Virus is the result of a successful immune response [25, 33]. Initially, detection of the infection triggers an innate immune response which targets the Virions that are released at the surface, as well as infected cells in the superficial layers. However, neither of these responses can permanently clear the Virus: elimination of basal cells require the presence of HPV-specific cytotoxic T-cells which are recruited during an adaptive immune response [24]. To model the specific targeting of infected basal cells, we need to introduce T cell mediated elimination of infected 8* and D* cells. For this purpose, we assume that the number of T cells recruited is proportional to the number of infected cells in the basal layer, 115* + 1119*. This implies that each infected cell is targeted at a constant rate (,4, referred to hereafter as the immune capacity,

Finally, since the number of basal cells is assumed to be conserved on average, the eliminated cell needs to be replaced. Therefore, elimination of an infected cell triggers proliferation of an 8 cell with probability 193, or of a 8*< cell with probability 193* = 1 — 195, see Fig. 2. The nature of these probabilities will be discussed in Results.

Next, we describe the replacement dynamics that ensue after elimination of infected basal cells. The underlying premise for the following replacement rules is conservation (on average) of basal stem cells. While there is, to our knowledge, no direct eXperimental evidence for conservation of basal stem cells during viral clearance, it has been shown that stem cell fate can change temporarily and reversibly after perturbation [19] , see also [28]. With this in mind, and as illustrated in Fig. 2A, we assume that D* cell elimination triggers division of a stem cell (infected or uninfected) in asymmetric or symmetric fashion, identically to the spontaneous division events depicted in Fig. 1B. On the other hand, 8* cell elimination is assumed to trigger a symmetric division, resulting in two identical 8 copies if an 8 cell is triggered, and into two 8* copies if an 8* cell is prompted to divide (Fig. 2B). In fact, if 8* division were to trigger asymmetric division of either an 8 or an 8* cell, this would violate our assumption of epithelial cell homeostasis because the total number of stem cells in the tissue would not be conserved over time.

In summary, the infection dynamics in the basal layer are described by a 4-dimensional continuous-time Markov process with variables 8, 8*, D and D*. The process is governed by the spontaneous division and migration events as in Fig. 1B, in conjunction with the immune-mediated events depicted in Fig. 2. Furthermore, we shall assume that each infection starts with an inoculum of no infected stem cells 8* at time t: 0.

Data

The longitudinal REACH data set [34] used in this study is publicly available and distributed by Sociometrics Corporation (http://www.socio.com/reachdata.php). Access to the nonsensitive part of the data (standard data) used in our study was granted by Sciometrics Corporation without IRB or other ethics approval.

To calibrate the models we used longitudinal infection data from The Reaching for Excellence in Adolescent Care and Health (REACH) project of the Adolescent Medicine HIV/AIDS Research Network [34]. Between 1996—2000, the REACH study followed 578 HIV-infected and HIV-uninfected adolescents (ages 13—18) in 13 US cites. A detailed description of the study is found in [35], see also [36]. In particular, the 411 female participants in the cohort were tested for HPV every 6 months, and we extracted the corresponding longitudinal data set from the master file. Not all data points were suitable for our purposes, so we made exclusions according to the following criteria: 31 participants did not test positive for HPV during the study; 31 participants had less than 2 valid HPV tests; 15 participants had less than one valid HPV test after the first HPV positive test; 6 participants had missing visit dates; 15 participants had incomplete information about the HPV subtypes. Among the remaining 313 females included in our analysis, 212 were HIV-positive, and 121 were HIV-negative. None of the participants changed HIV-status during the study period.

Results

Estimates from the parametric data analysis

For this purpose, we made the assumption of a well-mixed basal layer: upon removal of an infected cell, division of an 8 cell occurs with probability p3 = nS/(ns + 113*), and division of an 8* cell with probability 195* = nS*/(n3 + 113*). In other words, we assumed that the spatial arrangement of cells in the 2D basal layer can be ignored (the opposite end of the spectrum—a spatially clustered population of infected cells—is discussed below). Since the relative size of the infected population compared to the entire basal layer is small throughout the infection, 113 >> 113*, we can approximate 195 m 1 and 195* m 0. As outlined in section 2 in 81 Text, it follows that the 8* cell dynamics reduce to a subcritical branching process, The probability of survival to time tfor this process is, according to results in [32], In particular, addition of the immune capacity transforms the ~ 1/ t decay in (3) into an exponential decay.

Thereby, we faced the issue of non-identifia-bility of the model, a common problem in statistical inference. To understand where these issues arise, we first consider the probability density function f for the persistence time of the infection (see section 3 in 81 Text for its derivation)

From (7) we see that the values of l and 1’ cannot by inferred individually, and the best we can do is infer their product, a E Ar. Even though there are no further apparent identifiability issues, we found that for no large enough, the density (7) only depends on the ratio 05/ no (see section 4 in SI Text). As a consequence, a and no cannot be inferred individually, and we perform the inference over [,4 and no for fixed values of 05, across a prior range of biologically meaningful values a E [0.01, 0.25] d_1(see section 5 in SI Text for a justification of this range). In addition to the identif1ability issues, the MLE required the derivation of a nonstandard likelihood function that takes into account the different combinations of data types: infections were either present at the time of the first visit (prevalent infections), or they were initiated after the first visit (incident infections); some individuals were lost to followup before clearing the virus (right-censoring), and both the time of initiation and the time of clearance were only determined up to the between-visit intervals (interval-censoring). The derivation of the corresponding likelihood function is found in section 3 in SI Text.

In fact, it has been shown that longitudinal HPV studies bear a significant amount of misclassif1cations due to short-term variation [37], and that apparently cleared infections can reappear after variable amounts of time [38, 39]. The time before reappearance of seemingly cleared infections could be interpreted as a latency period during which the infection temporarily regresses to subdetection levels. However, molecular evidence for this latency mechanism has so far only been established in animal models [40]. Therefore, we decided to interpret the first of two consecutive negative test results as the time of clearance of the infection.

In what follows, the maximum likelihood estimates are denoted by a hat (A) on the parameter name, and subscripts (—) and (+) are used to refer to the HIV-negative and HIV-positive cohorts, respectively. As explained above, the number of initially infected cells no is a linear function of a, which varies over the prior range [0.01, 0.25]. The inferred ranges for the initial number of infected cells are n07_(a) E [5, Across the prior range of a, the inferred number of initially infected cells is slightly higher (but of the same order of magnitude) in HIV-positive compared to HIV-negative individuals:

p _.__<Ioom=_<m .

I_<I:momz<m

p _.__<I_uom_=<m

Parameter inference and clearance time distributions. A, B MLE inference for the initial number of infected cells no and immune capacity ,uforthe HIV-negative and HIV-positive cohorts, over the prior range of a. Estimators are denoted by a (A) symbol. Due to identifiability issues (see text), the estimators [so are linear functions of (1, whereas the estimators 11 are constant overthe prior range of a. The [so estimates are similar between the HIV-negative and HIV-positive cohorts (A), but there is a 100-fold difference forthe 11 estimates

C The parametric clearance time distributions for both cohorts are derived using equation (7) and the estimates from panels A and B. The distributions are insensitive to the choice of a, see section 6 in Text SI (plots shown for a = 0.14). The band inside the box is the median, the bottom and top of the box are 1st and 3rd quartile, respectively, and the whiskers correspond to the 5th and 95th percentiles, respectively.

+ (a) > fi07_(a), for all a. To our knowledge, there is no experimental data that would allow us to assess the validity of these model predictions. Regarding the immune capacity (,4, we find a stark difference between the cohorts: the estimated capacity [L in the HIV-negative cohort (fi_ = 1.4 - 10—3d‘1) is two orders of magnitude larger than the estimated capacity [1+ in the

Since the clearance time distributions were found to be insensitive to a over the prior range (see section 6 in Text SI), the distribution in Fig. 3C is only shown for an intermediate value of a = 0.14. Due to the reduced immune capacity in the HIV-positive cohort, its median clearance time is considerably larger (689 days) than the median time in the HIV-negative cohort (340 days).

Stochasticity vs immune response

Therefore, we compared the model-based persistence distributions for varying immune capacities (,4. As shown in Fig. 4A, the median time to clearance is a decreasing function of (,4, and the distributions become more localized with increasing ],4. However, comparing the distributions for [,4 = 0 and ,u /,11_ = 1 (where ,u = [L is the estimated immune capacity of HIV-negative individuals), the contribution of the immune response appears to be small in comparison to the contribution of the stochastic cell dynamics (compare the box plots for [,4 = 0 and ,u /,11_ = 1 in Fig. 4A). This is particularly clear when plotting the clearance probability as a function of time as shown in Fig. 4B. In particular, comparing the (,u /,11_ = 0)curve with the (,u /,11_ = 1)curve after 2 years, the clearance probability without immune response (0.66) is only ~ 17% smaller than the clearance probability with normal immune capacity (0.79). In other words, the stochastic dynamics contribute to as much as ~ 83% of the viral clearing mechanism in healthy individuals, and the contribution from the immune system is comparatively small.

Space and the impact of clustering

In this situation, elimination of an infected cell prompts division of an 8 cell with high probability, j ustify-ing the approximation p3 = 1 — 195* = 1. As a consequence, the persistence distribution could be derived analytically (6), rendering the model amenable to MLE inference. To assess whether the ensuing results were an artifact of the well-mixing assumption, we developed the following alternative model that takes into account the spatial clustering of infected cells.

That is, all the infected cells are inside the disk, whereas the outside is populated only by uninfected cells. Since the number of D* cells is roughly proportional to the number of 8* cells (see section 2.2 in Text

Immune capacity in the branching process model. A The model-derived distribution of the clearance time is shown as a function of varying immune capacity ,u (note that [1, corresponds to the immune capacity of HIV-negative individuals). The box-and-whisker plots correspond to the 5th, 25th, 50th, 75th and 95th percentiles, see caption of Fig. 3 for details. B The clearance probability as a function of time is shown for different levels of immune capacity.

Accordingly, Whenever an infected cell in the interior of the disk is eliminated by a T-cell, the probability to trigger an 8 cell division is given by the ratio of disk circumference to disk area: p8 = min { 1/ , #13,, 1} and 195* = 1 — 195. Under these assumptions, the 8* cell dynamics are now decoupled from the 8 cell dynamics, but they still depend on the D* cell dynamics, see also section 2 in 81 Text for details.

V

Immune capacity in the spatial model. A The model-derived clearance time distribution is shown as a function of varying immune capacity ,u for the spatial model. The box-and-whisker plots correspond to the 5th, 25th, 50th, 75th and 95th percentiles, see capture of Fig. 3. The dotted lines are the three quartiles from u = ft, in Fig. 4A. B The clearance probability as a function of time is shown for different levels of immune capacity. All results in this figure are estimates based on 1000 simulations per ,u-value.

As in Fig. 4A for the well-mixed model, we investigated the impact of increasing immune capacity [,4 on the clearance time distribution in Fig. 5A. We make the following observations. First, time to clearance is generally longer in the branching process model: the three dotted horizontal lines correspond to the three quartiles for the (M/fi_ = 1)distribution in Fig. 4. Only for the (,u/fi_ = 8) -distri-bution, which corresponds to an 8-fold increase in immune capacity, are all three quartiles of the spatial model below the corresponding quartiles of the branching process model. Second, the impact of the immune capacity on the clearance time for the clustered model is even weaker than in the well-mixed model. Whereas the well-mixed model yields a decrease in median time to clearance for increasing (,4, small [,4 values yield a slight increase in median clearance time for the spatial version. This is due to the fact that, in contrast to the branching process model, elimination of an infected cell can trigger division of an 8* cell (with probability 195* > 0), therefore compensating for the loss of the infected cell and delaying clearance. The relative insensitivity of the persistence time distribution to [,4 is further illustrated in Fig. 5B, where we observe that the clearance probability is only slightly increased for small [,4 values. Finally, since the prior estimates of several model parameters have a relatively large interval of uncertainty (see section 5 in 81 Text), we performed a combined sensitivity analysis. By means of a Monte-Carlo simulation (with the parameters r, a, p and [,4 drawn from their prior ranges), we computed the corresponding persistence time distribution, and found that it did not substantially differ from the fixed parameter distribution (see section 7 in 81 Text for details).

Discussion

To date, little attention has been paid to the potential role of homeostatic cell dynamics in clearing HPV infections. In this study, we combined mechanistic mathematical models at the cellular level with epidemiological data at the population level to disentangle the respective roles of immune capacity and cell dynamics in the clearing mechanism. Our results suggest that chance—in form of the stochastic dynamics of basal stem cells—plays a critical role in the elimination of HPV-infected cell clones. In particular, we found that in individuals with normal immune capacity (HIV-negative cohort), the immune response may contribute to less than 20% of the clearing task overall—the rest is taken care of by the random succession of symmetric and asymmetric stem cell divisions. Furthermore, in immunocompromised individuals (HIV-positive cohort) the contribution of the immune response is likely to be negligible.

First, in view of the high prevalence of HPV infections and the relatively small risk of persistent infections that eventually lead to malignant disease, the identification of predictive markers for persistence would be valuable [8]. However, if stochasticity does indeed play a key role in viral clearance, and if the major difference between individuals who clear effectively and individuals who develop persistent infections is largely a matter of chance, then there may not be any predictive markers to discover. Hence, we may want to rephrase the question, and ask if there is a way of modulating the cellular dynamics to achieve an increase in the clearance probability. Our results suggests that by increasing either the probability of a symmetric division (r) or the proliferation frequency (A) through a locally administered drug, time to clearance and risk of progression could be substantially reduced.

Since estrogen stimulates [42] and progesterone inhibits [43] epithelial proliferation, it seems plausible that a decrease in cervical proliferation could be caused directly via increased progesterone levels, and indirectly via loss of the estrogenic mid-cycle peak. The resulting decrease in proliferation (smaller 1) would imply an increase in time to clearance and a higher risk of progression to cancer. While the same reasoning would imply an increased risk of cervical cancer in progestin-only users, the effect of progestin-only contraceptives on HPV persistence and cervical cancer development is less consistent in the literature [44, 45]. This highlights the need for future research into the influence of sex steroids on the natural history of oncogenic HPV infection.

To date, the existence of latent infections has been demonstrated in animal models, and it is assumed to occur in HPV infections as well. The current theory of latency is based on the assumption that the virus stays present inside long-lived basal stem cells [48]. But while the notion of such long-lived, asymmetrically dividing and slow-cy-cling stem cells is consistent with a theory of epithelial homeostasis developed in the 1970’s [49] , it is not aligned with the new paradigm that is based on fast-cycling stem cells that divide both asymmetrically and symmetrically [18, 19, 29]. According to our model, which is based on this more recent theory of homeostasis, viral latency is again a stochastic phenomenon and occurs if the number of infected cells becomes very small (latent period) before growing back to a detectable size. A more thorough discussion of the latency issue will be the subject of future work.

In two recent studies [50, 51] , deterministic (partial) differential equation models were used to study evolutionary and ecological aspects of HPV infections and competition between coexisting HPV types. To our knowledge, we are the first to develop a stochastic model of HPV infection that couples stem cell dynamics with viral infection and immune response. In addition, the methods introduced here provide a useful tool in the parametric analysis of longitudinal data sets that contain both prevalent (present at study begin) and incident (initiation after study begin) infections, as well as right-censoring (study exit before viral clearance) and interval-censoring (duration of infection only known up to an interval). In fact, 70% of the individuals in the analyzed REACH data set had an unknown time of initiation, rendering conventional nonparametric approaches problematic (see section 1 in S1 Text for details). Thanks to the mechanistic models introduced and analyzed in this study, we were able to account for the unknown time lapse between infection initiation and study entry. Finally, the approach employed in this study may prove useful in other situations. In fact, mathematical models at the tissue-level are often difficult to parametrize because sample sizes in pathology studies are generally small and exhibit large between-patient variation. By combining longitudinal population-level data with cell-level mechanistic models as done in this study, insights can be gained across the scales.

First, it is known that there can be time-lags between inoculation and productive infection [22]. Since these lag times vary widely among individuals, and since we wanted to avoid adding to the complexity of the model, we set the incubation period to zero. Second, since infected cells acquire a selective growth advantage only at later, symptomatic stages of the infection [31], we assumed that the presence of viral DNA did not alter the proliferation rates of infected stem cells. In addition, there is, to our knowledge, no experimental evidence regarding HPV-mediated modulation of symmetric and asymmetric division patterns in infected tissues. Third, we assumed that the interactions between virus and immune system are independent of the specific HPV strains, and that there are no synergistic or competitive effects among co-infecting types, see also [51]. Since we believe that adding these more subtle aspects would not change the main conclusion of the importance of stochasticity, we did not incorporate them into the current model. However, we plan to address these issues in future work. Fourth, a more realistic alternative to the clustered model version is provided by explicitly spatial models with lattice-based voter dynamics [52, 53]. Such a spatial model extension is subject of ongoing work. Fifth, even though the stratified squamous epithelia at different anogenital and oropharyngeal sites affected by HPV are qualitatively similar, we parametrized our model for cervical infections, and our insights regarding the role of stochastic stem cell proliferation in viral clearance may not apply to other organs. Finally, our model predicts extinction of infection with probability 1 due to the subcritical nature of the process. This is not in contradiction with the observation that a small fraction of infections persist and progress. In fact, progression from HPV infection to sustained neoplastic growth is associated with cellular changes triggered by the viral genome. These transformations are themselves stochastic processes, and hence progression only takes place in the small group of individuals where the oncogenic transformation takes place before extinction of the infected population.

Supporting Information

A priori parameter estimates for cervical squamous epithelium; 6. Sensitivity of clearance time to choice of a; 7. Sensitivity analysis for spatial model (PDF)

Acknowledgments

MDR would like to thank Profs. Stanley Robboy and Mark Stoler for stimulating discussions and correspondences on the pathology of the female reproductive tract.

Author Contributions

Performed the experiments: MDR. Analyzed the data: MDR. Wrote the paper: MDR ERM RD.