Viral and immune dynamics of genital human papillomavirus infections in young women with high temporal resolution

Abstract

Human papillomavirus (HPV) infections drive one in 20 new cancer cases, exerting a particularly high burden on women. Most anogenital HPV infections are cleared in less than two years, but the underlying mechanisms that favour persistence in around 10% of women remain largely unknown. Notwithstanding, it is precisely this information that is crucial for improving treatment, screening, and vaccination strategies. To understand viral and immune dynamics in non-persisting HPV infections, we set up an observational longitudinal cohort study with frequent on-site visits for biological sample collection. We enrolled 189 women aged from 18 to 25 and living in the area of Montpellier (France) between 2016 and 2020. We performed 974 on-site visits for a total of 1,619 months of follow-up. We collected data on virus load, local immune cell populations, local concentrations of cytokines, and circulating antibody titres. Using hierarchical Bayesian statistical modelling to simultaneously analyse the data from 164 HPV infections from 76 participants, we show that in two months after infection, HPV viral load in non-persisting infections reaches a plateau that lasts on average for 13 to 20 months (95% credibility interval) and is then followed by a rapid clearance phase. This first description of the dynamics of HPV infections comes with the identification of immune correlates associated with infection clearance, especially gamma-delta T cells and CXCL10 concentration. A limitation of this study on HPV kinetics is that many infection follow-ups are censored. Furthermore, some immune cell populations are difficult to label because cervical immunity is less well characterised than systemic immunity. These results open new perspectives for understanding the frontier between acute and chronic infections, and for controlling HPV-associated diseases, as well as for research on human cancers of infectious origin. Trial Registration: This trial was registered is registered at ClinicalTrials.gov under the ID NCT02946346. This study has been approved by the Comité de Protection des Personnes (CPP) Sud Méditerranée I (reference number 2016-A00712-49); by the Comité Consultatif sur le Traitement de l'Information en matière de Recherche dans le domaine de la Santé (reference number 16.504); by the Commission Nationale Informatique et Libertés (reference number MMS/ABD/ AR1612278, decision number DR-2016-488), by the Agence Nationale de Sécurité du Médicament et des Produits de Santé (reference 20160072000007).

Fig 1. Unsupervised clustering of flow cytometry data from cervical smears.

(A) Highlighting 11 homogeneous populations using UMAP clustering on data from samples without HPV or positive for a “focal” HPV infection. (B) Comparison of the clusters frequencies based on infection status, and (C) cluster annotation and FC values shown in panels A and B. FC were calculated by an abundance analysis with the diffcyt-DA-edgeR method adjusted with a Benjamini–Hochberg test. In panel B, each dot represents a single sample. Raw differential abundance results can be found in Table C in S1 Supplementary Materials. The code to generate this figure can be found in https://doi.org/10.57745/KJGOYZ. FC, fold change; HPV, human papillomavirus; UMAP, uniform manifold approximation and projection.

Fig 2. Local immune response associated with HPV infections.

(A) Association between HPV focal status and five cytokines. We show the posterior distributions of a Bayesian model in which the log10 of each cytokine or chemokine is a predictor of the HPV status (focal or uninfected), assuming a correlated random effect for each participant (see Tables F and S7 in S1 Supplementary Materials for details). Thick lines and shaded areas show the 80% and 95% credibility intervals respectively. (B) For the HPV focal infections only, correlation matrix between the local density of five cytokines and the proportion of the 11 cell clusters from Fig 2. β represents the regression coefficient of linear regression with each cell population as the response variable (see Methods). For example, a β of 0.14 can be interpreted as follows: “a 10 percentage point increase in the CD8 T cells frequency is associated with a 1.4 percentage point increase in the IFNγ concentration.” A similar correlation matrix for HPV–negative samples can be found in Fig D in S1 Supplementary Materials. The code to generate this figure can be found in https://doi.org/10.57745/KJGOYZ. HPV, human papillomavirus.

Fig 3. Modelling virus load kinetics for 158 HPV genital infections in 76 women.

The dots indicate the data points and, therefore, the window of observation for each participant. Each panel corresponds to one participant and shows the number of HPV genome copies per number of human genome copies resulting from a three-slopes hierarchical Bayesian model. The lines show the trajectory of the posterior median value and shaded area the 95% CrI. Open circles indicate values below the limit of detection, open triangles indicate samples where the virus was detected but the viral load could not be estimated (see the Methods). Letters before the anonymity number (above each panel) indicate whether the participant was vaccinated (V) or not (N). The code to generate this figure can be found in https://doi.org/10.57745/KJGOYZ. CrI, credibility interval; HPV, human papillomavirus.

Fig 4. Population predictions of viral load kinetics in HPV infections.

(A) Representation of the six parameters governing the descriptive models: plateau virus load ψ_vl, time of the infection ψ_t0, growth phase duration ψ_gr, plateau duration ψ_p, clearance phase duration ψ_cl, and initial viral load v₀. (B) The lines and points show the observed kinetics. The shaded areas show the population quantiles of the viral load through time simulated from the posterior distributions, with the median in yellow and the 10% and 90% deciles in grey (for all three, shaded areas indicate the 90% CrI). The time alignment was done using the posterior median time of infection. (C) Standard deviation of the random effects for the virus load and the infection duration related to the heterogeneity among virus genotypes (in black) or among hosts (in orange). (D) Correlation between these two random effects. The line represents the 95% quantiles, the thick box the interquartile range, and the dot the median of the posterior distribution. The code to generate this figure can be found in https://doi.org/10.57745/KJGOYZ. CrI, credibility interval; HPV, human papillomavirus.

Fig 5. Modelling local immune dynamics during HPV infections.

(A) Illustration of the parameters governing the descriptive models for a hypothetical case with two genotypes (A in dark blue and B in light blue) infecting an individual. For each infection, we have virus kinetics (top) and the immune dynamics model (bottom) is governed by four parameters: the peak of the response ψ_tm, the time to the peak ψ_delay, the beginning of the infection t_start, and its end t_clear. The first 2 are estimated, whereas the last 2 originate from the viral kinetics model. At the individual level, the effects of each infection are summed (in purple), assuming a basal level ψ₀. For the cytokines, this trajectory reflects the log of the normalised concentration. The trajectory can also decrease throughout the infection (i.e., ψ_tm∈ℝ). (B) Illustration of the additional steps for the FCM model. We consider the same two virus genotypes A and B and two hypothetical cell populations: 1, whose trajectory computation is illustrated in panel A, and 2, generated in the same way, with the same t_start, and t_clear, but different ψs. For the FCM, the frequencies of each cell population are considered on the centred log-ratio (clr) scale. To compare our model to the count data, the clr trajectory is transformed back into frequencies using the softmax function and then multiplied by the total number of observed cells in the sample (see the Methods). The code to generate this figure can be found in https://doi.org/10.57745/KJGOYZ. FCM, flow cytometry; HPV, human papillomavirus.

Fig 6. Local immune dynamics during HPV infections for ten randomly selected participants.

Each row represents a local immune variable, and each column represents one participant. The dots indicate the observed data, the lines indicate the trajectory of the posterior median value, and the shaded area the 95% CrI from the bootstrapped posterior distribution. We also show the total protein concentration and total number of leukocytes as indicators of the sample quality for the cytokines assay and the FCM, respectively. The code to generate this figure can be found in https://doi.org/10.57745/KJGOYZ. CrI, credibility interval; FCM, flow cytometry; HPV, human papillomavirus.

Fig 7. Estimates of the key metric for local immune dynamics during the HPV infections.

(A) Aggregated posterior distribution of the median peak time of each immune variable, compared to the prior distribution. The x-axis was scaled to the estimated average infection duration, i.e., 16 months (with a 95% CrI of 13 to 20). The colour indicates the median peak fold change compared to the baseline value. (B) Distribution of the median leukocyte frequency or cytokine concentration fold change compared to the baseline, at the time of the peak. The colours indicate the value of the baseline frequency. The vertical line represents the 95% quantiles, the thick box the interquartile range, and the dot the median of the bootstrapped posterior distribution. The code to generate this figure can be found in https://doi.org/10.57745/KJGOYZ. CrI, credibility interval; HPV, human papillomavirus.

Fig 8. Associations between HPV infection intensity and the immune response.

(A) Geometric mean IgG and IgM titers during the infection for vaccinated participants (yellow) or not (red). Antibodies were only measured in 56 infections (see the main text). (B) Regression coefficients of a linear model between IgG serostatus and viral load AUC. For each infection, “0” indicates seronegativity, “1” seropositivity in at least one sample, and “NA” missing information (because of censored follow-ups). The first value indicates the status before the infection, and the second, the status after the infection, whence the notations “0 0,” “0 1,” or “NA 1.” All vaccinated individuals were seropositive and treated as a separate class (“vax”). The reference is “0 0,” i.e., individuals who never seroconverted during the infection. (C) Regression coefficients of a linear model between the geometric mean antibody titer measured during infection and viral load AUC. (D) Distribution of the bootstrapped correlation between the individual effect on viral load AUC and the individual effect on local immune variable mean value during the infection. The panel colours indicate the median peak time (also shown in Fig 7A) for an average infection. In A, the three components of the box and whisker represent the median, interquartile range, and 95% range. The other panels represent the median aggregated posterior distribution from 500 bootstrapped regressions from the viral kinetics fit (the dot), the interquartile range (the thick line), and the 95% CrI (the thin line). Immune response variables and viral load AUC were centred and reduced to compare the effects. The code to generate this figure can be found in https://doi.org/10.57745/KJGOYZ. CrI, credibility interval; HPV, human papillomavirus.