Assessing the Impact of Prophylactic Vaccines on HPV Prevalence

ABSTRACT

A mathematical model considering female and male individuals is proposed to evaluate vaccination strategies applied to control of HPV transmission in human population. The basic reproductive number of the disease, R ₀ is given by the geometric mean of the basic reproductive number of female and male populations. The model has a globally asymptotically stable disease-free equilibrium whenever $R_0<1$. Furthermore, it has an unique endemic state when R ₀ exceeds unity which is globally asymptotically stable. Numerical simulations were done to compare several different vaccination schedules. The results showed that the vaccination strategies that do not include vaccination of men can only control the disease if more than 90% of women are vaccinated. The sensitivity analysis indicated that the relevant parameters to control HPV transmission, in order of importance, are vaccine efficacy times the fraction of population that is vaccinated, disease recovery-rate, and disease transmission rate. Therefore, health politics that promoting the increase of vaccine coverage, and screening for the disease in both population can improve disease control.

Keywords:
sensitivity analysis; basic reproductive number; ordinary differential equation

RESUMO

Um modelo matemático á proposto para avaliar as estratégias de controle da transmissão do HPV na população humana (homens e mulheres). O número reprodutivo básico da doença, R ₀, é dado pela média geométrica do número reprodutivo básico obtido nas populações femininas e masculinas. O modelo tem um equilíbrio livre da doença globalmente assintoticamente estável sempre que $R_0<1$. Além disso, tem um estado endêmico único quando $R_0>1$ que é globalmente assintoticamente estável. Simulações numéricas foram feitas para comparar vários esquemas de vacinação diferentes. Os resultados mostraram que as estratégias de vacinação que não incluem a vacinação de homens controlam a transmissão da doença se mais de 90 % das mulheres forem vacinadas. A análise de sensibilidade indicou que os parâmetros relevantes para controlar a transmissão do HPV, por ordem de importância, são a eficácia da vacina vezes a fração da população vacinada, a taxa de recuperação da doença e a taxa de transmissão da doença. Portanto, políticas de saúde que promovam o aumento da cobertura vacinal e a detecção da doença em ambas as populações podem melhorar o controle da doença.

Palavras-chave:
análise de sensibilidade; número reprodutivo básico; equações diferencias ordinárias

1 INTRODUCTION

Human Papillomavirus (HPV) is a sexually transmitted infection common in young sexually active individuals. It is projected that between 50 to 80% of these individuals will acquire the infection during their lives. Global prevalence in women is estimated to be around 10.4%, ranging from 2 to 44% ⁵. Most HPV infections are asymptomatic and self-limited. The natural recovery from infection could take from months to years, and the dominant immunological response involves neutralizing antibodies ¹¹. Around 10 to 20% of infections persist and they are strong associated with the development of cervical cancer in women as well as genital and anal cancer in men ^{4), (14}.

Bivalent, quadrivalent and nonavalent prophylactic vaccines against HPV types 16/18, 6/11/16/18, and 6/11/16/18/31/33/45/52/58, respectively, are already available ¹⁶, with vaccine efficacy ranging from 43 to 93% ¹³. The bivalent and quadrivalent vaccines were licenced in 2006 only for women, and from 2009 the quadrivalent vaccine was also licenced for men. The nonavalent vaccine was licenced for males and females in 2014 ¹.

Initially the vaccines were applied in a 3-dose schedule, but currently a 2-dose schedule with six months apart is enforced ⁷. The change in vaccination schedule was due to evidences that the antibody response generated after 2 dose is enough to prevent virus infection ^{21), (22}; besides vaccination cost is significantly reduced in this 2 dose-scheme. Vaccination target-individuals are between 9 to 14 years age ⁵ because exposure to infection is higher at younger ages, with a peak after the debut of sexual activity.

Data from cohort of vaccinated individuals suggest that the vaccine offers protection for at least 10-12 years ^{9), (15}. Moreover, unlike natural infection, the vaccine confers cross-protection against heterologous HPV types. It is interesting to mention that reduction of HPV 6/11/16/18 infection, genital warts, and cytological cervical abnormalities have been reported in the post-vaccination era ¹⁰.

Theoretical studies have been done to understand the epidemiological aspects of HPV infection. In particular, mathematical models have played an important role in estimating the impact of vaccination on HPV transmission. Some examples are ^{2), (3), (8), (17), (20}. In general, these works explore the impact of vaccination on HPV types that interacted in a synergistic or antagonist way. They concluded that vaccination can increase (competing types of virus) or decrease (synergistic interactions) the prevalence of the non-target type depending on the type of the interaction among virus ⁸. Both benefit of male vaccination and behaviour change were discussed by measuring its impact on R ₀; when waning immunity is included in the model, disease elimination is more difficult to be achieved ³. Heterogenity in disease tranmission rate associated with class age, and age-dependent impact of vaccination were also adressed. In this case, the model was parametrized using data from Toronto, Canadá. They concluded that vaccination of a single age cohort in one gender can eventually control the spread of the disease across age groups ². Combined impact of pap cytology screaning and mass vaccination was assessed to show that this approach can decrease disease prevalence faster that only vaccination ¹⁷. Susceptibility to HPV infection and risk of developing cancer were used to divide the population in two groups; and the vaccine was applied to the first group. The authors showed that in a scheme of three-doses vaccine, compliance is very important to achieve disease eradication, and sensitivity analysis showed that disease transmission probability is the most important parameter that impacts R ₀.

In this work we formulate a simple model that mimics HPV transmission in a sex-structured population where 2-dose vaccination schemes are implemented. Our purpose is to compare the impact of different vaccination programs on disease prevalence. Following the epidemiology of the infection, the female and male populations is divided in susceptible, infected, vaccinated and recovered individuals. The HPV basic reproductive number associated to the each sex is obtained, numerical simulations exemplify disease transmission dynamics and captures disease prevalence, and a sensitivity analysis highlights the most important parameters associated to R ₀. The full analysis of the steady states is done. This simple model captures disease dynamics and can be used to discuss strategies of HPV control, highlighting the important parameters related to the disease transmission, geographic diferences in disease prevalence, and target groups.

The paper is structured as follows. In section 2 a model is formulated to describe the HPV transmission given a vaccination schedule. It is assumed that the vaccine or the natural infection induces long-lasting immunity; besides, two doses are required to protect individuals from getting the disease. In section 3 the existence and stability of the model equilibria is investigated. In section 4 numerical simulations and a sensitivity analysis are carried out to evaluate the disease control strategies. Finally, the discussion is given in section 5.

2 FORMULATION OF THE MODEL

We denote by ${\overline{N}}_h$ the human population sexually active that comprise individuals between 9 to 69 years at time t. This population is further sub-divided into two groups: females, ${\overline{N}}_f$ and males, ${\overline{N}}_m$ which are recruited at constant rates Λ_f , and Λ_m , respectively, and are subject to a constant mortality rate µ. Then, ${\overline{N}}_f$, and ${\overline{N}}_m$ grow according to

$\frac{d{\overline{N}}_f}{dt}={\Lambda }_f-\mu {\overline{N}}_f,\frac{d{\overline{N}}_m}{dt}={\Lambda }_m-\mu {\overline{N}}_m.$ (2.1)

The model subdivide the female (f) and male (m) populations according to the epidemiology of the infection into susceptible individuals ${\overline{S}}_f$ and ${\overline{S}}_m$; individuals that receive the first and second dose of the vaccine, ${\overline{V}}_f$ and ${\overline{V}}_m$; infected individuals with HPV, ${\overline{I}}_f$, ${\overline{I}}_m$; and recovered individuals, ${\overline{Z}}_f$, ${\overline{Z}}_m$. We want to emphasize that the reason why the model only consider a single class of vaccinated for each sex is because both doses are applied before the individuals become sexually active, and the probability of being infected during the period between the two vaccinations is practically null.

We assume that a proportion σ_f1 and σ_m1 , respectively, of women and men in the age range of 9 to 14 are vaccinated for the first time, and six months later, a proportion σ_f2 , σ_m2 of these women and men receive the second dose. We denote by (₁ and (₂ the efficiency of the first and second dose of vaccine, that is, the fraction of individuals that are protected against the infection. According to the demographic assumptions given above, $q=e^{-\mu \eta }$ with $\eta =6$ denotes the probability of an individual to survive up to six months. Therefore, ${\Lambda }_i{\epsilon }_1{\sigma }_{i_1}q{\epsilon }_2{\sigma }_{i_2}$, with $i=f,m$, are the number of individuals per unit of time of each sex that are successfully vaccinated for the first time, survived up to six months (interval between vaccine doses), and completed the vaccination schedule taking an efficient second dose. These individuals are fully protected and have lifelong immunity. Observe that ${\Lambda }_i{\epsilon }_1{\sigma }_{i_1}\left(1-q{\epsilon }_2{\sigma }_{i_2}\right),i=f,m$ are the number of individuals of each sex that only get the first doses.

Susceptible women and men can be infected at rates p _f c _f β _f /N _m , and p _m c _m β _m /N _f , respectively, where c _i is the number of sexual encounters per year, β _f , β _m are the infection transmissions (probability that an infectious contact gives rise to an infection) from men to women and from women to men, respectively, and $p_i=e^{-\mu {\eta }_i}$, with $i=f,m$, are the proportion of individuals that survives until age η _i and become sexually active at this age. This susceptible group comprises individuals that did not receive the vaccine, individuals that received only one dose, and individuals that receive two doses but do not developed protective immunity. The infected women and men recover at rates γ _f , and γ _m , respectively. Figure 1 shows the diagram of the compartmental model. According to the above assumptions we obtain the following system of ordinary differential equations that mimics the dynamics of HPV transmission and vaccination:

$$\begin{gathered} \frac{d{\overline{S}}_f}{dt}={\Lambda }_f\left(1-{\epsilon }_1{\sigma }_{f1}q{\epsilon }_2{\sigma }_{f2}\right)-\frac{p_f{c}_f{\beta }_f}{{\overline{N}}_m}{\overline{S}}_f{\overline{I}}_m-\mu {\overline{S}}_f \\ \frac{{\displaystyle d{\overline{V}}_f}}{{\displaystyle dt}}={\Lambda }_f{\epsilon }_1{\sigma }_{f1}q{\epsilon }_2{\sigma }_{f2}-\mu {\overline{V}}_f \\ \frac{{\displaystyle d{\overline{I}}_f}}{{\displaystyle dt}}=\frac{{\displaystyle p_f{c}_f{\beta }_f}}{{\displaystyle {\overline{N}}_m}}{\overline{S}}_f{\overline{I}}_m-\left({\gamma }_f+\mu \right){\overline{I}}_f \\ \frac{{\displaystyle d{\overline{Z}}_f}}{{\displaystyle dt}}={\gamma }_f{\overline{I}}_f-\mu {\overline{Z}}_f \\ \frac{{\displaystyle d{\overline{S}}_m}}{{\displaystyle dt}}={\Lambda }_m\left(1-{\epsilon }_1{\sigma }_{m1}q{\epsilon }_2{\sigma }_{m2}\right)-\frac{{\displaystyle p_m{c}_m{\beta }_m}}{{\displaystyle {\overline{N}}_f}}{\overline{S}}_m{\overline{I}}_f-\mu {\overline{S}}_m \\ \frac{{\displaystyle d{\overline{V}}_m}}{{\displaystyle dt}}={\Lambda }_m{\epsilon }_1{\sigma }_{m1}q{\epsilon }_2{\sigma }_{m2}-\mu {\overline{V}}_m \\ \frac{{\displaystyle d{\overline{I}}_m}}{{\displaystyle dt}}=\frac{{\displaystyle p_m{c}_m{\beta }_m}}{{\displaystyle {\overline{N}}_f}}{\overline{S}}_m{\overline{I}}_f-\left({\gamma }_m+\mu \right){\overline{I}}_m \\ \frac{{\displaystyle d{\overline{Z}}_m}}{{\displaystyle dt}}={\gamma }_m{\overline{I}}_m-\mu {\overline{Z}}_m. \end{gathered}$$ (2.2)

Figure 1:
Diagram of the compartmental model. The dashed lines highlights disease transmission between female (f) and male (m) populations. Recruitment occurs among individuals from 9 to 14 years, and disease transmission comprises only individuals sexually active (from 15 to 69 years old). The epidemiological classes are susceptible ${\overline{S}}_i$, vaccinated ${\overline{V}}_i$, infected ${\overline{I}}_i$, and recovered ${\overline{Z}}_i$ individuals, with $i=f,m;a_i={\epsilon }_1{\sigma }_{i_1}q{\epsilon }_2{\sigma }_{i_2}$, and ${\lambda }_i=p_i{c}_i{\beta }_i{\overline{I}}_i/{\overline{N}}_i$.

From equation (2.1) we can see that

${\overline{N}}_f\left(t\right)\to \frac{{\Lambda }_f}{\mu }=N_f,\text{and}{\overline{N}}_m\left(t\right)\to \frac{{\Lambda }_m}{\mu }=N_m$

when $t\to \infty$. Therefore, N _f and N _m are upper bounds of ${\overline{N}}_f\left(t\right)$ and ${\overline{N}}_m\left(t\right)$, respectively, provided that ${\overline{N}}_f\left(t\right)\le N_f$ and ${\overline{N}}_m\left(t\right)\le N_m$. Furthermore, if ${\overline{N}}_i\left(t\right)>N_i$ then ${\overline{N}}_i\left(t\right)$ will decrease to $N_i,i=f,m$. Then, the feasible region is given by

$\Omega =\left\{\left({\overline{S}}_f,{\overline{V}}_f,{\overline{I}}_f,{\overline{Z}}_f,{\overline{S}}_m,{\overline{V}}_m,{\overline{I}}_m,{\overline{Z}}_m\right)\in R_{+}^8|{\overline{N}}_f\le N_f,{\overline{N}}_m\le N_m\right\}.$

To facilitate the analysis and interpretation of results, we take proportions

$S_i={\overline{S}}_i/N_i,V_i={\overline{V}}_i/N_i,I_i={\overline{I}}_i/N_i,Z_i={\overline{Z}}_i/N_i,$

with $S_i,V_i,I_i,i=f,m$ satisfying the following ODE system

$$\begin{gathered} \frac{d{S}_f}{dt}=\mu \left(1-{\epsilon }_1{\sigma }_{f1}q{\epsilon }_2{\sigma }_{f2}\right)-p_f{c}_f{\beta }_f{S}_f{I}_m-\mu S_f \\ \frac{{\displaystyle d{V}_f}}{{\displaystyle dt}}=\mu {\epsilon }_1{\sigma }_{f1}q{\epsilon }_2{\sigma }_{f2}-\mu V_f \\ \frac{{\displaystyle d{I}_f}}{{\displaystyle dt}}=p_f{c}_f{\beta }_f{S}_f{I}_m-\left({\gamma }_f+\mu \right)I_f \\ \frac{{\displaystyle d{S}_m}}{{\displaystyle dt}}=\mu \left(1-{\epsilon }_1{\sigma }_{m1}q{\epsilon }_2{\sigma }_{m2}\right)-p_m{c}_m{\beta }_m{S}_m{I}_f-\mu S_m \\ \frac{{\displaystyle d{V}_m}}{{\displaystyle dt}}=\mu {\epsilon }_1{\sigma }_{m1}q{\epsilon }_2{\sigma }_{m2}-\mu V_m \\ \frac{{\displaystyle d{I}_m}}{{\displaystyle dt}}=p_m{c}_m{\beta }_m{S}_m{I}_f-\left({\gamma }_m+\mu \right)I_m, \end{gathered}$$ (2.3)

and $Z_i=1-S_i-V_i-I_i$.

3 ANALYSIS OF THE MODEL

3.1 Disease-free equilibrium and R ₀

The model (2.3) has the disease free equilibrium, $P_0=(1–a_f,a_f,0,1–a_m,a_m,0)$ where

$a_f={\epsilon }_1{\sigma }_{f1}q{\epsilon }_2{\sigma }_{f2}\text{and}{a}_m={\epsilon }_1{\sigma }_{m1}q{\epsilon }_2{\sigma }_{m2}.$

This equilibrium corresponds to the state where both populations of men and women are free of HPV infection. From the first and fourth equations of system (2.3) it can be seen readily that $0\le S_i\le 1–a_i,i=f,m$.

The Basic Reproductive Number, denoted by R ₀, represents the average number of secondary cases that one primary infective individual generates over the course of its infectious period in a whole susceptible population. It is a threshold condition that determines when an epidemy can occur, or when a disease can remain endemic. If $R_0<1$, the disease dies out; otherwise if $R_0>1$ the disease persists. In ⁶, the basic reproductive number of a disease is mathematically defined as the spectral ratio of the next-generation operator associated to the disease-free equilibrium. The next-generation operator corresponding to model (2.3) is given by the product of two matrices: the non-negative matrix of the new infection terms in the susceptible population denoted by K, and the inverse of the matrix of the transfer terms denoted by T, which are given by

$K=\left(\begin{array}{cc}0& b_f\left(1-a_f\right)\\ b_m\left(1-a_m\right)& 0\end{array}\right)$

and

$T=\left(\begin{array}{cc}{\gamma }_f+\mu & 0\\ 0& {\gamma }_m+\mu \end{array}\right),$

with

$b_f=p_f{c}_f{\beta }_{f}and{b}_m=p_m{c}_m{\beta }_m,$

respectively. Then, it follows that

$R_0=\rho \left(K{T}^{-1}\right)=\sqrt{R_f{R}_m}$ (3.1)

where

$R_f=\frac{b_f\left(1-a_f\right)}{{\gamma }_f+\mu },R_m=\frac{b_m\left(1-a_m\right)}{{\gamma }_m+\mu }$ (3.2)

and ρ denotes the spectral radius of KT ^-1.

Using Theorem 2 of ²³, the following result is established.

Lemma 1. The disease-free equilibrium, P ₀ of model (2.3) is locally-asymptotically stable (LAS) whenever $R_0<1$, and unstable if $R_0>1$.

The above lemma shows that VPH infection can be effectively controlled if the initial populations of the model are in the basin of attraction of the disease-free equilibrium. To ensure the elimination of disease regardless of initial population sizes, a global stability proof for the disease-free equilibrium is needed. This is done below, using a suitable Lyapunov function.

Theorem 1. Let $R_0<1$, then the disease-free equilibrium, P ₀ of model (2.3) is globally-asymptotically stable (GAS) in Ω.

Proof To prove this statement consider the following function defined in Ω.

$V=A{I}_f+I_m,$ (3.3)

where $A=\frac{b_m\left(1-a_m\right)}{{\gamma }_f+\mu }$. The orbital derivative of V is given by

$$\begin{gathered} \dot{V}=A\left[b_f{S}_f{I}_m-\left({\gamma }_f+\mu \right)I_f\right]+\left[b_m{S}_m{I}_f-\left({\gamma }_m+\mu \right)I_m\right] \\ =\left({\gamma }_m+\mu \right)\left[\frac{A{b}_f{S}_f}{{\gamma }_m+\mu }-1\right]I_m+\left({\gamma }_f+\mu \right)\left[\frac{b_m{S}_m}{A\left({\gamma }_f+\mu \right)}-1\right]I_f \\ \le \left({\gamma }_m+\mu \right)\left[\frac{A{b}_f\left(1-a_f\right)}{{\gamma }_m+\mu }-1\right]I_m+\left({\gamma }_f+\mu \right)\left[\frac{b_m\left(1-a_m\right)}{A\left({\gamma }_f+\mu \right)}-1\right]I_f \\ =\left({\gamma }_m+\mu \right)\left(R_0^2-1\right)I_m \end{gathered}$$ (3.4)

which is less or equal to zero in Ω. We see that $\dot{V}=0$ if and only if $I_m=0$. From system (2.3) it can be seen that $I_m=0$ implies $I_f\left(t\right)\to 0$ and $S_i\left(t\right)\to 1-a_i,i=f,m$ when $t\to \infty$. Therefore, the largest compact invariant subset contained in $\left\{(S_f,V_f,I_f,S_m,V_m,I_m)\in \Omega |\dot{L}=0\right\}$ is the disease-free equilibrium P ₀. Therefore, from the LaSalle-Lyapunov Theorem ¹² it follows that all trajectories starting in Ω approach P ₀.

The epidemiological implication of the above result is that the disease will be eliminated from the community if $R_0<1$, regardless of the initial number of infected individuals.

3.2 Existence and stability of endemic equilibrium point

The endemic equilibrium is given by $P_1=\left(S_f^{*},a_f,I_f^{*},S_m^{*},a_m,I_m^{*}\right)$ where

$$\begin{gathered} S_f^{*}=\frac{\mu \left(1-a_f\right)-\left({\gamma }_f+\mu \right)I_f^{*}}{\mu }, \\ {S}_m^{*}=\frac{\mu \left(1-a_m\right)-\left({\gamma }_m+\mu \right)I_m^{*}}{\mu }, \\ {I}_m^{*}=\frac{\mu b_m\left(1-a_f\right)I_f^{*}}{\left({\gamma }_m+\mu \right)\left(\mu +b_m{I}_f^{*}\right)}, \\ {I}_f^{*}=\frac{\mu \left({\gamma }_m+\mu \right)\left(R_f{R}_m-1\right)}{b_m\left({\gamma }_m+\mu +b_f\left(1-a_m\right)\right)}. \end{gathered}$$

This equilibrium corresponds to the state where the disease is endemic in both populations of men and women. Then we have the following result.

Theorem 2. For $R_m{R}_f>1$, or equivalently $R_0>1$, there exists a unique endemic equilibrium P ₁.

Global stability of the endemic equilibrium P ₁ is proved next using a suitable Lyapunov function.

Theorem 3. Let $R_0>1$, then the endemic equilibrium P ₁ of system (2.3) is globally asymptotically stable.

Proof. Assume $R_0>1$. To prove the global stability of the equilibrium $P_1=\left(S_f^{*},a_f,I_f^{*},S_m^{*},a_m,I_m^{*}\right)$ we will show that the function

$$\begin{gathered} W=d_1\left(S_f-S_f^{*}-S_f^{*}\ln \frac{S_f}{S_f^{*}}\right)+d_2\left(V_f-V_f^{*}-V_f^{*}\ln \frac{V_f}{V_f^{*}}\right) \\ +d_3\left(I_f-I_f^{*}-I_f^{*}\ln \frac{I_f}{I_f^{*}}\right)+d_4\left(S_m-S_m^{*}-S_m^{*}\ln \frac{S_m}{S_m^{*}}\right) \\ +d_5\left(V_m-V_m^{*}-V_m^{*}\ln \frac{V_m}{V_m^{*}}\right)+d_6\left(I_m-I_m^{*}-I_m^{*}\ln \frac{I_m}{I_m^{*}}\right) \end{gathered}$$

with $d_i>0,i=1,...,6$, is a Lyapunov function for the equilibrium P ₁, that is:

1. $W\left(X\right)\ge 0$, where $X=(S_f,V_f,I_f,S_m,V_m,I_m)\in \Omega$, and $W\left(P_1\right)=0$,

2. The orbital derivative $\dot{W}\left(X\right)<0$.

It is clear that the function W satisfies condition i). To prove ii) we calculate the orbital derivative of W which is given by

$$\begin{gathered} \dot{W}\left(X\right)=d_1\left(1-\frac{S_f^{*}}{S_f}\right)\left(\mu \left(1-a_f\right)-b_f{S}_f{I}_m-\mu S_f\right)+d_2\left(1-\frac{V_f^{*}}{V_f}\right)\left(\mu a_f-\mu V_f\right) \\ +d_3\left(1-\frac{I_f^{*}}{I_f}\right)\left(b_f{S}_f{I}_m-\left({\gamma }_f+\mu \right)I_f\right)+d_4\left(1-\frac{S_m^{*}}{S_m}\right)\left(\mu \left(1-a_m\right)-b_m{S}_m{I}_f-\mu S_m\right) \\ +d_5\left(1-\frac{V_m^{*}}{V_m}\right)\left(\mu a_m-\mu V_m\right)+d_6\left(1-\frac{I_m^{*}}{I_m}\right)\left(b_m{S}_m{I}_f-\left({\gamma }_m+\mu \right)I_m\right). \end{gathered}$$

From system (2.3) at equilibrium we obtain the following relations

$$\begin{gathered} \left(1-a_f\right)\mu =b_f{S}_f^{*}{S}_m^{*}+\mu S_f^{*} \\ {a}_f=V_f^{*} \\ {\gamma }_f+\mu =\frac{b_f{S}_f^{*}{I}_m^{*}}{I_f^{*}} \\ \left(1-a_m\right)\mu =b_m{S}_m^{*}{I}_f^{*}+\mu S_m^{*} \\ {a}_m=V_m^{*} \\ {\gamma }_m+\mu =\frac{b_m{S}_m^{*}{I}_f^{*}}{I_m^{*}}. \end{gathered}$$

Substituting the parameters above in the derivative (3.5), and taking

$d_1=d_3,d_4=d_6,d_4=\frac{b_f{S}_f^{*}{I}_m^{*}}{b_m{S}_m^{*}{I}_f^{*}},$

after several calculations and simplifications, we obtain that $\dot{W}$ can be written as

$$\begin{gathered} \dot{W}\left(X\right)=-d_1\mu \frac{{\left(S_f-S_f^{*}\right)}^2}{S_f}-d_2\mu \frac{{\left(V_f-V_f^{*}\right)}^2}{V_f}-d_4\mu \frac{{\left(S_m-S_m^{*}\right)}^2}{S_m}-d_5\mu \frac{{\left(V_m-V_m^{*}\right)}^2}{V_m}{V}_m \\ -4d_1{b}_f{S}_f^{*}{I}_m^{*}\left[\frac{1}{4}\left(\frac{S_f^{*}}{S_f}+\frac{S_f{I}_m{I}_f^{*}}{S_f^{*}{I}_m^{*}{I}_f}+\frac{S_m^{*}}{S_m}+\frac{S_m{I}_f{I}_m^{*}}{S_m^{*}{I}_f^{*}{I}_m}\right)-1\right]. \end{gathered}$$ (3.5)

Using the fact that the arithmetic mean is greater or equal to the geometric mean, it can be readily seen that the last term is less than zero, and it is zero if and only if $(S_f,V_f,I_f,S_m,V_m,I_m)=P_1$. It follows that $\dot{W}\left(X\right)\le 0$, and $\dot{W}\left(X\right)=0$ if and only if $X=P_1$. Therefore, trajectories in the interior of Ω approach P ₁.

4 NUMERICAL SIMULATIONS AND SENSITIVITY ANALYSIS

In this section we present the numerical simulations of system (2.3) and a sensitivity analysis of R ₀ using the parameter values given in Table 1 (unless otherwise stated) to evaluate the prevalence of the disease as well as the effectiveness of several vaccination strategies.

Table 1:
Summary of the model variables, their description and initial conditions of the simulations ^{1), (5}.

Variable	Description	Initial condition (time $\mathrm{t}=0$)
S f	proportion of susceptible females at time t	$S_f=0.73$
V f	proportion of vaccinated females at time t	$V_f=0.0$
I f	proportion of infected females at time t	$I_f=0.27$
Z f	proportion of recover females at time t	$Z_f=0$
S m	proportion of susceptible males at time t	$S_m=0.64$
V m	proportion of vaccinated males at time t	$V_m=0.0$
I m	proportion of infected males at time t	$I_m=0.36$
Z m	proportion of recover males at time t	$Z_m=0$

Figures 2 to 5 show the temporal evolution of each epidemiological class for distinct set of parameters that represent different scenarious of disease transmission. In all figures the initial conditions are $S_f=0.73$, $I_f=0.27$, $S_m=0.64$,$I_m=0.36$, and $V_f=V_m=0$. In Figures 2 and 3 the disease prevails in female and male populations because $R_0>1$, while in Figures 4 and 5, $R_0<1$, and the disease dies out. It is interesting to notice in Figure 3 that disease is endemic for both male and female population although $R_m=0.98<1$; on the other hand, the disease can dies out from the population, even when $R_f>1$ or $R_m>1$ as it is shown in Figure 4, where $R_f=1.26$.

Figure 2:
Temporal dynamics of the populations. The parameters are ${\sigma }_{f1}\times {\sigma }_{f2}=0.4,{\sigma }_{m1}\times {\sigma }_{m2}=0.28,{\epsilon }_1\times {\epsilon }_2=0.9,c_f\times {\beta }_f=0.012$ days^-1, $c_m\times {\beta }_m=0.006$ days^-1, $\mu =3.65\times {10}^{-5}$ days^-1, ${\eta }_f={\eta }_m=16.2$ years, ${\gamma }_f={\gamma }_m=0.0028$ days^-1, and $\eta =6$ months. The corresponding thresholds are $R_f=2.19,R_m=1.28$, and $R_0=1.67$

Figure 3:
Temporal dynamics of the populations. The parameters are ${\sigma }_{f1}\times {\sigma }_{f2}=0.4,{\sigma }_{m1}\times {\sigma }_{m2}=0.4,{\epsilon }_1\times {\epsilon }_2=0.9,c_f\times {\beta }_f=0.012$ days^-1, $c_m\times {\beta }_m=0.0053$ days^-1, $\mu =3.65\times {10}^{-5}$ days^-1, ${\eta }_f={\eta }_m=16.2$ years, ${\gamma }_f={\gamma }_m=0.0028$ days^-1, and $\eta =6$ months. The corresponding thresholds are $R_f=2.19,R_m=0.98$, and $R_0=1.46$

Figure 4:
Temporal dynamics of the populations. The parameters are ${\sigma }_{f1}\times {\sigma }_{f2}=0.5,{\sigma }_{m1}\times {\sigma }_{m2}=0.6,{\epsilon }_1\times {\epsilon }_2=0.9,c_f\times {\beta }_f=0.008$ days^-1, $c_m\times {\beta }_m=0.0053$ days^-1, $\mu =3.65\times {10}^{-5}$ days^-1, ${\eta }_f={\eta }_m=16.2$ years, ${\gamma }_f={\gamma }_m=0.0028$ days^-1, and $\eta =6$ months. The corresponding thresholds are $R_f=1.26,R_m=0.7$, and $R_0=0.94$

Figure 5:
Temporal dynamics of the populations. The parameters are ${\sigma }_{f1}\times {\sigma }_{f2}=0.7,{\sigma }_{m1}\times {\sigma }_{m2}=0.6,{\epsilon }_1\times {\epsilon }_2=0.9,c_f\times {\beta }_f=0.008$ days^-1, $c_m\times {\beta }_m=0.0053$ days^-1, $\mu =3.65\times {10}^{-5}$ days^-1, ${\eta }_f={\eta }_m=16.2$ years, ${\gamma }_f={\gamma }_m=0.0028$ days^-1, and $\eta =6$ months. The corresponding thresholds are $R_f=0.7$, and $R_0=0.77$.

Figure 6 shows the partial rank correlation coefficient (PRCC) obtained from a global sensitivity analysis ¹⁸. The input parameters were sampled using the latin hypercube sampling with ${\epsilon }_1{\epsilon }_2{\sigma }_{f_1}{\sigma }_{f_2}={\epsilon }_1{\epsilon }_2{\sigma }_{m_1}{\sigma }_{m_2}=\mathrm{Unif}(0,0.9)$ ^{1),(16),(19),(22}, $c_f{\beta }_f=\mathrm{Unif}(0.0073,0.0123)$ days^{-1(4), (14}, $c_m{\beta }_m=\mathrm{Unif}(0.0053,0.0097)$ days^{-1 (4), (14} ${\eta }_f={\eta }_m=\mathrm{Unif}(14,20)$ years ¹⁹, and ${\gamma }_f={\gamma }_m=\mathrm{Unif}(0.0013,0.0055)$ days^{-1 (11), (21} where Unif means the uniform probability distribution. The other parameters are $\mu =3.65\times {10}^{-5}$ days^-1 and $\eta =6$ months ⁹. We run 5000 simulations that correspond to different input parameter sets, and the output chosed for the sensitivity analysis is R ₀.

Figure 6:
Sensitivity analysis using R ₀ as the output. The input parameters were chosen from an uniform distribution using the latin hypercube sampling. A negative-control, dummy parameter, was used to assign a zero value for a sensitivity index.

Table 2:
Summary of model parameters, their description and units ^{1), (4), (9), (11), (14), (16), (19), (21), (22}.

Parameter	Description	Value	Interval	Unit
μ	natural mortality rate	$3.65\times {10}^{-5}$	-	days-1
cf βf	transmission rate from male to female per person, per month	0.29	0.22 - 0.37	months-1
cm βm	transmission rate from female to male per person, per month	0.21	0.16 - 0.29	month-1
${\sigma }_{f1}\times {\sigma }_{f2}$	proportion of vaccinated females with 2-dose	0.40	0 - 1.0	-
${\sigma }_{m1}\times {\sigma }_{m2}$	proportion of vaccinated males with 2-dose	0.28	0 - 1.0	-
γf , γm	recovered rates for female and male	1.0	0.5 - 2.0	years-1
${\epsilon }_1\times {\epsilon }_2$	vaccine efficacy	0.9	0.45 - 1.0	-
ηf , ηm	age of sexual initiation for female and male	16.2	14 - 20	years
η	interval between the first and second vaccine dose	6.0	6 - 12	months

Figure 7 shows the disease prevalence versus the proportion of immunized females, ${{\sigma }_f}_1{{\sigma }_f}_2$ (assuming that ${\epsilon }_1{\epsilon }_2=0.9$), where the disease prevalence is defined as the proportion of infected females and males in the steady state, $I_f^{*}+I_m^{*}$. Each curve in the figure corresponds to a chosen proportion of immunized males, ${\sigma }_{m_1}{{\sigma }_m}_{{}_2}$.

Figure 7:
Disease prevalence versus proportion of immunized girls. The parameters are σm ₁ × σ _m2 = 0.6, (₁ × (₂ = 0:9, c _f × β _f = 0:011 days^-1, c _m × β _m = 0.006 days^-1, $\mu =3.65\times {10}^{-5}$ days^-1, η _f = η _m = 16.2 years, γ _f = γ _m = 0.0028 days^-1, and $\eta =6$ months. Each curve corresponds to a given proportion of immunized boys.

Figures 8a-b show the disease prevalence versus the transmission rates c _f β _f and c _m β _m . The others parameters in the simulations are $\mu =3.65\times {10}^{-5}$, ${\epsilon }_1\times {\epsilon }_2=0.9,{\eta }_f={\eta }_m=16.2$ years, and $\eta =6$ months. Each curve corresponds to a chosen recover rate, γ _f , and it is assumed that ${\gamma }_m={\gamma }_f$.

Figure 8:
Disease prevalence versus transmission rate. From top to bottom: γ _f = γ _m =0.0013, 0.0018, 0.0023, 0.0028, 0.0033 in days^-1. In (a) c _m β _m = 0.007 days^-1 and in (b) c _f β _f = 0.0096 days^-1. The other parameters are $\mu =3.65\times {10}^{-5}$ days^-1, (₁ × (₂ = 0.92, σ f ₁ × σ f ₂ = 0.4, σm ₁ × σ _m2 = 0.28

5 DISCUSSION

Routine vaccination against HPV started in 2006 as a 3-dose schedule targeting women from 11 to 25 years. But, although all effort made by the health authorities to convince people about the importance of vaccination, the number of individuals taking the third dose was very small. Due to this, and to the discovery that with a 2-dose schedule the antibody titters are enough to promote protection for at least 10 to 12 years, a 2-dose scheme was implemented for females from 9 to 14 years with an interval of 6 months between doses. Besides, vaccination was also implemented in men with the objective of increasing herd immunity.

In this work we model HPV infection dynamics considering female and male population from 9 to 69 years old. The recruitment of susceptible individuals comprises boys and girls between 9 to 14 years old. This class of age is suppose do not yet had it sexual debut (this is the premise of the vaccination strategy currently in course). The proportion of the susceptible population that participates in the transmission chain consist of: (i) individuals that did not receive the vaccine, (ii) individuals that received only 1-dose, (iii) and individuals that receive 2-dose, but did not reach levels of protective immunity. Individuals that are successfully immunized are assumed to be life-long protected (the HPV vaccination was started in 2006 and until now the individuals immunized still showing high level of antibodies). Infected individuals recover (spontaneous or with treatment) from infection after a period of time; natural infection also induces long-lasting immunity.

Our main objective was to explore the effectiveness of a vaccine schedule implemented in young population. For this, a deeply search in the literature was done to construct Table 1 which summarize the important parameters related to disease transmission. We found an expression for the basic reproductive number of the disease, R ₀, corresponding to the average of the basic reproductive numbers of female and male population R _f , and R _m , respectively, which in turn depend on the parameters that characterize virus-host dynamics and vaccination efficacy. Thus, if $R_0<1$ the disease dies out; otherwise if $R_0>1$ it persists. Persistence of the disease occurs when one of the populations (or both) is an infection source, i.e. $R_m>1$ or $R_f>1$ or (see Figures 1-4). Our findings show that control strategies that include even a small percentage of boys increase significantly the herd immunity, and therefore reduce the disease prevalence (see Figure 6). Therefore, the best strategy to control HPV is to vaccinate girls and boys.

We performed a sensitivity analysis to identify the key parameters affecting disease prevalence (see Figure 5). Our study indicates that increasing both vaccine efficacy, (₁(₂, or vaccination compliance, ${\sigma }_{f_i}{\sigma }_{f_i}i=f,m$, is more efficient than decreasing transmission rate c _f β _f or c _m β _m . In particular, the contribution of male and female transmission rates in R ₀, are almost equal, which emphasize the importance of including the boys as a target population for vaccination. The second more important parameter in disease control is early clearance of the virus either by the immune system, or by treatment (by reducing γ _i , see Figures 5 and 7). Finally, a delayed of age of sexual debut negatively impacts R ₀; and it can explain the variability of disease prevalence among regions as well as social/economic groups. This information could help to design the minimum age of vaccination, and to emphasize vaccination campaign among target groups.

ACKNOWLEDGEMENTS

CPF thanks grant 16/23738-3, São Paulo Research Foundation (FAPESP). LE acknowledgments support from project IN-112713, PAPIIT-UNAM. Part of this work was developed during the visit of CPF and FA to Facultad de Ciencias at UNAM. Both thanks project IN-112713, PAPIIT-UNAM.

REFERENCES

• ¹
  Center for disease control and prevention. https://www.cdc.gov/nchs/products/databriefs/db280.htm
  » https://www.cdc.gov/nchs/products/databriefs/db280.htm
• ²
  M. Al-arydah & R.J. Smith. An age-structured model of human papilloma virus vaccination. Mathematics and Computers in Simulation, 82 (2011), 629-652.
• ³
  V. Brown & K.A.J. White. The HPV vaccination strategy: could male vaccination have a significant impact? Comput. Math. Methods. Med, 11 (2010), 223-237.
• ⁴
  A.N. Burchell, F. Coutleé, P. Tellier, J. Hanley & E.L. Franco. Genital Transmission of Human Papillomavirus in Recently Formed Heterosexual Couples. The Journal of Infectious Diseases, 204 (2011), 1723-1729.
• ⁵
  A.N. Burchell, R.L. Winer, S. de Sanjosé & E.L. Franco. Epidemiology and transmission dynamics of genital HPV infection. Vaccine, 24 (2006), DOI:10.1016/j.vaccine.2006.05.031.
  » https://doi.org/10.1016/j.vaccine.2006.05.031
• ⁶
  O. Diekmann & J.A.P. Heesterbeek. “Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation”. Wiley (2000).
• ⁷
  R. Donken, T.M.S. Klooster, R.M. Schepp, F.R.M. van der Klis, M.J. Knol, C.J.L.M. Meijer & H.E. de Melker. Immune Responses After 2 Versus 3 Doses of HPV Vaccination up to 4 1/2 Years After Vaccination: An Observational Study Among Dutch Routinely Vaccinated Girls. JID, 215 (2017), 359-367.
• ⁸
  E.H. Elbasha & A.P.E. Galvani. Vaccination against multiple HPV types. Math. Biosci, 197 (2005), 88-117.
• ⁹
  C. Fraser, J.E. Tomassini, L. Xi, G. Golm, M. Watson, A.R. Giuliano, E. Barr & K.A. Ault. Modeling the long-term antibody response of a human papillomavirus (HPV) virus-like particle (VLP) type 16 prophylactic vaccine. Vaccine, 25 (2007), 4324-4333.
• ¹⁰
  S.M. Garland, S.K. Kjaer, N. Muñoz, S.L. Block, D.R. Brown, M.J. DiNubile, B.R. Lindsay, B.J. Kuter, G. Perez, G.D.A.J. Saah, R. Drury, D. Das & C. Velicer. Impact and Effectiveness of the Quadrivalent Human Papillomavirus Vaccine: A Systematic Review of 10 Years of Real-world Experience. Clinical Infectious Diseases, 63 (2016), 519-527.
• ¹¹
  P.E. Gravitt & R.L. Winer. Natural History of HPV Infection across the Lifespan: Role of Viral Latency. Viruses, 9 (2017), DOI:10.3390/v9100267.
  » https://doi.org/10.3390/v9100267
• ¹²
  J. Hale & H. Ko¸cak. “Dynamics and Bifurcations”. Springer-Verlag (1991).
• ¹³
  D.M. Harpera & L.R. De Mars. HPV vaccines - A review of the first decade. Gynecologic Oncology, 146 (2017), 196-204.
• ¹⁴
  B.Y. Hernandez, L.R. Wilkens, X. Zhu, P. Thompson, M. McDuffi, Y.B. Shvetsov, L.E. Kamemoto, J. Killeen, L. Ning & M.T. Goodman. Transmission of Human Papillomavirus in Heterosexual Couples. Emerging Infectious Diseases, 14 (2008), 888-894.
• ¹⁵
  S.K. Kjaer, M. Nygard, J. Dillner, J.B. Marshall, D. Radley, M. Li, C. Munk, B.T. Hansen, L.G. Sigurdardottir, M. Hortlund, L. Tryggvadottir, A. Joshi, R. Das & A.J. Saah. A 12-Year Follow-up on the Long-Term Effectiveness of the Quadrivalent Human Papillomavirus Vaccine in 4 Nordic Countries. Clinical Infectious Diseases , 66 (2018), 339-345.
• ¹⁶
  M. Lehtinen, C. Lagheden, T. Luostarinen, T. Eriksson, D. Apter, K. Harjula, M. Kuortti, K. Natunen, J. Palmroth, T. Petaja, E. Pukkala, M. Siitari-Mattila, F. Struyf, P. Nieminen, J. Paavonen, G. Dubin & J. Dil. Ten-year follow-up of human papillomavirus vaccine efficacy against the most stringent cervical neoplasia end-point-registry-based follow-up of three cohorts from randomized trials. BMJ Open, 7 (2017), DOI:10.1136/bmjopen-2017-015867.
  » https://doi.org/10.1136/bmjopen-2017-015867
• ¹⁷
  T. Malik, J. Reimer, A. Gumel, E.H. Elbasha & S. Mahmud. The impact of an imperfect vaccine and pap cytology screening on the transmission of human papillomavirus and occurrence of associated cervical dysplasia and cancer. Math Biosci Eng., 13 (2013), 1173-1205.
• ¹⁸
  S. Marino, I.B. Hogue, C.J. Ray & D.E. Kirschner. A methodology for performing global uncertainty and sensitivity analysis in systems biology. J Theor Biol., 254 (2008), 178-196.
• ¹⁹
  E.Y. Petrosky, G. Liu, S. Hariri & L.E. Markowitz. Human Papillomavirus Vaccination and Age at First Sexual Activity, National Health and Nutrition Examination Survey. Clin Pediatr (Phila), 56 (2017), 363-370.
• ²⁰
  O. Sharomi & T. Malik. A model to assess the effect of vaccine compliance on Human Papillomavirus infection and cervical cancer. Applied Mathematical Modelling, 47 (2017), 528-550.
• ²¹
  M. Stanley. Immune responses to human papillomavirus. Vaccine, 24 (2006), DOI:10.1016/j.vaccine.2005.09.002.
  » https://doi.org/10.1016/j.vaccine.2005.09.002
• ²²
  M. Stanley. HPV - immune response to infection and vaccination. Infectious Agents and Cancer, 5 (2010), 19-25.
• ²³
  P. van den Driessche & J. Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci., 180 (2002), 29-48.

Publication Dates

• Publication in this collection
  16 Sept 2019
• Date of issue
  May-Aug 2019

History

• Received
  19 Oct 2018
• Accepted
  18 Feb 2019