The results we obtained show that models assuming that seropositive individuals are fully and permanently protected from reinfection with HPV-16 are clearly inferior to the other models making no such strong assumptions. This conclusion is based on DIC scores. It is important to realise that DIC does not detect a ‘correct’ model in terms of HPV-16 transmission mechanism. Instead, it provides a quantitative model ranking which discourages complexity and is based on the ability of models under consideration (among which the ‘correct’ model may not even be present) to be fitted to the data. Hence, if a simpler model can be calibrated to the data at least as well as a more complex model, it will get a better DIC score. To receive a better DIC ranking, a more complex model would have to justify its complexity by producing a notably better fit than its simpler competitors. To further clarify the context in which our results should be viewed, we mention that our results can be meaningfully interpreted only if we completely rely on the available data – should these be extended or replaced, our results would inevitably change too. Another important aspect is that the DIC ranking factors in how well the models can be fitted to all data at once, for both males and females. If we, for instance, restricted ourselves to only calibrating the models to HPV seroprevalence, the resulting model ranking would likely be different.
As is evident from Table 1, our ‘best’ model is SIS2, closely followed by SIRS4. The difference in DIC scores between the two models is not substantial and hence does not imply that SIS2 is clearly preferable. We should note that the reason why SIS2 outscored SIS1 is inclusion of seroreversion. Indeed, it is the only difference between the models. Seroreversion in SIS2 is implemented with the help of two additional parameters (rsr,m and rsr,f), as compared with SIS1, and nonetheless, it improved the fit substantially enough to overcome penalisation for extra parameters and get ahead of SIS1 by 8.5 points. The benefits of seroreversion in SIS1 are predictable since without it, SIS1 can not capture declining seroprevalence in older females. For the same reason, SIRS4 provided a significant improvement over SIRS3. We see that seroreversion in SIS and SIRS models is crucial in terms of improving the fit to data, even though the rate of seroreversion is low.
Although the highest ranking models SIS2 and SIRS4 have different structures, as we mentioned in Results, the fitted durations of full natural immunity in SIRS4 are very short. Hence, this model is approaching a limit case when it almost becomes SIS2 (see Figure 1 and Figure 3).
It is our view, given what is currently known about immunity (in particular, the reported association between seropositivity and reductions in the number of incident infections in seropositive individuals [34, 35]), that the protection mechanism assumed in SIS2 may be a more realistic representation of naturally acquired protective immunity than a short but full immunity as in SIRS4.
It is important to note that nearly all information available regarding the possible association of seropositivity with protective immunity has come from studies of females. The only study of males in this context that we are aware of  suggests that for males seropositivity is possible without any immunity. No substantiated inferences in regard to the existence of protective immunity in males resulted from our study: SIS2 was not sensitive to variations in the degree of immunity in males. To increase sensitivity, the amount of data used for model specification and calibration and/or their accuracy should be increased, which we expect to happen in future, when, for example, HPV DNA seroprevalence data for males become available.
Our models have a number of limitations. In particular, we assumed the duration of immunity to be the same for all ages, which is unlikely to be true in reality, and the probability of seroconversion to be independent of an individual’s age though there is some evidence to the contrary . Additionally, compartmental models are inherently biased in certain respects. Because they assume a sexual contact is effectively instantaneous, to achieve better fit to real data, compartmental models need to compensate for a somewhat lowered level of sexual activity by maintaining higher probabilities of transmission and longer durations of infection (see  for detailed discussion). It is also important to remember that considerable uncertainty remains in our understanding of HPV natural history which influence our specification of priors for model parameters. Also, reliability of data obtained from sexual behavior surveys may be arguable. Finally, the results of this study rely on the data we calibrated our models to, which had their own limitations (see [17, 39] for discussion). Perhaps, the most evident limitation is that HPV-16 prevalence data only covered women aged 15–39.