i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 305 — #1 i i i i i i Tema Tendências em Matemática Aplicada e Computacional, 20, N. 2 (2019), 305-321 © 2019 Sociedade Brasileira de Matemática Aplicada e Computacional www.scielo.br/tema doi: 10.5540/tema.2019.020.02.0305 Assessing the Impact of Prophylactic Vaccines on HPV Prevalence F. AZEVEDO1, L. ESTEVA2 and C.P. FERREIRA3* Received on October 19, 2018 / Accepted on February 18, 2019 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, R0, 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 R0 < 1. Furthermore, it has an unique endemic state when R0 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. 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 in- fection during their lives. Global prevalence in women is estimated to be around 10.4%, ranging from 2 to 44% [5]. Most HPV infections are asymptomatic and self-limited. The natural recov- ery from infection could take from months to years, and the dominant immunological response involves neutralizing antibodies [11]. 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]. *Corresponding author: Claudia Ferreira – E-mail: claudia.pio@unesp.br – https://orcid.org/ 0000-0002-9404-6098 1 Faculdade de Computação e Engenharia Elétrica, UNIFESSPA, 68507-590, Marabá, Pará, Brasil E-mail: franazevedo@unifesspa.edu.br 2 Facultad de Ciencias, UNAM, 04510 Ciudad Universitaria, DF, México. E-mail: lesteva@ciencias.unam.mx 3São Paulo State University (UNESP), Institute of Biosciences, Department of Biostatistics, 18618-689 Botucatu, SP, Brazil. E-mail: claudia.pio@unesp.br https://orcid.org/0000-0002-9404-6098 https://orcid.org/0000-0002-9404-6098 i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 306 — #2 i i i i i i 306 ASSESSING THE IMPACT OF PROPHYLACTIC VACCINES ON HPV PREVALENCE 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 [16], with vac- cine efficacy ranging from 43 to 93% [13]. 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 [1]. Initially the vaccines were applied in a 3-dose schedule, but currently a 2-dose schedule with six months apart is enforced [7]. 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 [5] 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 [10]. 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 ex- plore 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 (synergis- tic interactions) the prevalence of the non-target type depending on the type of the interaction among virus [8]. Both benefit of male vaccination and behaviour change were discussed by mea- suring its impact on R0; when waning immunity is included in the model, disease elimination is more difficult to be achieved [3]. 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 [2]. Combined impact of pap cytology screaning and mass vaccination was assessed to show that this approach can decrease disease prevalence faster that only vaccination [17]. 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 R0. 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 Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 307 — #3 i i i i i i AZEVEDO, ESTEVA and FERREIRA 307 prevalence, and a sensitivity analysis highlights the most important parameters associated to R0. 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 trans- mission 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 get- ting 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 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, N̄ f and males, N̄m which are recruited at constant rates Λ f , and Λm, respectively, and are subject to a constant mortality rate µ . Then, N̄ f , and N̄m grow according to dN̄ f dt = Λ f −µN̄ f , dN̄m dt = Λm−µN̄m. (2.1) The model subdivide the female (f) and male (m) populations according to the epidemiology of the infection into susceptible individuals S̄ f and S̄m; individuals that receive the first and second dose of the vaccine, V̄f and V̄m; infected individuals with HPV, Ī f , Īm; and recovered individuals, Z̄ f , 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 ε1, and ε2 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−µη with η = 6 denotes the probability of an individual to survive up to six months. Therefore, Λiε1σi1qε2σi2 , 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 Λiε1σi1(1−qε2σi2), 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 /Nm, and pmcmβm/N f , respectively, where ci 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 Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 308 — #4 i i i i i i 308 ASSESSING THE IMPACT OF PROPHYLACTIC VACCINES ON HPV PREVALENCE women to men, respectively, and pi = e−µη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 indi- viduals 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 compart- mental model. According to the above assumptions we obtain the following system of ordinary differential equations that mimics the dynamics of HPV transmission and vaccination: dS̄ f dt = Λ f (1− ε1σ f1qε2σ f2)− p f c f β f N̄m S̄ f Īm−µ S̄ f dV̄f dt = Λ f ε1σ f1qε2σ f2 −µV̄f dĪ f dt = p f c f β f N̄m S̄ f Īm− (γ f +µ)Ī f dZ̄ f dt = γ f Ī f −µZ̄ f (2.2) dS̄m dt = Λm(1− ε1σm1qε2σm2)− pmcmβm N̄ f S̄m Ī f −µ S̄m dV̄m dt = Λmε1σm1qε2σm2 −µV̄m dĪm dt = pmcmβm N̄ f S̄m Ī f − (γm +µ)Īm dZ̄m dt = γm Īm−µZ̄m. From equation (2.1) we can see that N̄ f (t)→ Λ f µ = N f , and N̄m(t)→ Λm µ = Nm when t→ ∞. Therefore, N f and Nm are upper bounds of N̄ f (t) and N̄m(t), respectively, provided that N̄ f (t) ≤ N f and N̄m(t) ≤ Nm. Furthermore, if N̄i(t) > Ni then N̄i(t) will decrease to Ni, i = f ,m. Then, the feasible region is given by Ω = {(S̄ f ,V̄f , Ī f , Z̄ f , S̄m,V̄m, Īm, Z̄m) ∈ R8 +|N̄ f ≤ N f , N̄m ≤ Nm}. To facilitate the analysis and interpretation of results, we take proportions Si = S̄i/Ni, Vi = V̄i/Ni, Ii = Īi/Ni, Zi = Z̄i/Ni, Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 309 — #5 i i i i i i AZEVEDO, ESTEVA and FERREIRA 309 S f Z fVf Λ(1 − a f ) Λ a f f λ γ f I f I m a m )Λ(1 − µ µ µ µ µ µ V S fΛ a λ γ Z m m m m m µ µ 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 S̄i, vaccinated V̄i, infected Īi, and recovered Z̄i individuals, with i = f ,m; ai = ε1σi1qε2σi2 , and λi = piciβi Īi/N̄i. with Si,Vi, Ii, i = f ,m satisfying the following ODE system dS f dt = µ(1− ε1σ f1qε2σ f2)− p f c f β f S f Im−µS f dVf dt = µε1σ f1qε2σ f2 −µVf dI f dt = p f c f β f S f Im− (γ f +µ)I f (2.3) dSm dt = µ(1− ε1σm1qε2σm2)− pmcmβmSmI f −µSm dVm dt = µε1σm1qε2σm2 −µVm dIm dt = pmcmβmSmI f − (γm +µ)Im, and Zi = 1−Si−Vi− Ii. Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 310 — #6 i i i i i i 310 ASSESSING THE IMPACT OF PROPHYLACTIC VACCINES ON HPV PREVALENCE 3 ANALYSIS OF THE MODEL 3.1 Disease-free equilibrium and R0 The model (2.3) has the disease free equilibrium, P0 = (1−a f ,a f ,0,1−am,am,0) where a f = ε1σ f1qε2σ f2 and am = ε1σm1qε2σ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≤ Si ≤ 1−ai, i = f ,m. The Basic Reproductive Number, denoted by R0, 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 R0 < 1, the disease dies out; otherwise if R0 > 1 the disease persists. In [6], 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 matri- ces: 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 = ( 0 b f (1−a f ) bm(1−am) 0 ) and T =  γ f +µ 0 0 γm +µ  , with b f = p f c f β f and bm = pmcmβm, respectively. Then, it follows that R0 = ρ(KT−1) = √ R f Rm (3.1) where R f = b f (1−a f ) γ f +µ , Rm = bm(1−am) γm +µ (3.2) and ρ denotes the spectral radius of KT−1. Using Theorem 2 of [23], the following result is established. Lemma 1. The disease-free equilibrium, P0 of model (2.3) is locally-asymptotically stable (LAS) whenever R0 < 1, and unstable if R0 > 1. Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 311 — #7 i i i i i i AZEVEDO, ESTEVA and FERREIRA 311 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 elimina- tion 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 R0 < 1, then the disease-free equilibrium, P0, of model (2.3) is globally-asymptotically stable (GAS) in Ω. Proof To prove this statement consider the following function defined in Ω. V = AI f + Im, (3.3) where A = bm(1−am) γ f +µ . The orbital derivative of V is given by V̇ = A [ b f S f Im− (γ f +µ)I f ] + [ bmSmI f − (γm +µ)Im ] = (γm +µ) [ Ab f S f γm +µ −1 ] Im +(γ f +µ) [ bmSm A(γ f +µ) −1 ] I f (3.4) ≤ (γm +µ) [ Ab f (1−a f ) γm +µ −1 ] Im +(γ f +µ) [ bm(1−am) A(γ f +µ) −1 ] I f = (γm +µ)(R2 0−1)Im which is less or equal to zero in Ω. We see that V̇ = 0 if and only if Im = 0. From system (2.3) it can be seen that Im = 0 implies I f (t)→ 0, and Si(t)→ 1− ai, i = f ,m when t → ∞. Therefore, the largest compact invariant subset contained in {(S f ,Vf , I f ,Sm,Vm, Im) ∈ Ω |L̇ = 0} is the disease-free equilibrium P0. Therefore, from the LaSalle-Lyapunov Theorem [12] it follows that all trajectories starting in Ω approach P0. The epidemiological implication of the above result is that the disease will be eliminated from the community if R0 < 1, regardless of the initial number of infected individuals. 3.2 Existence and stability of endemic equilibrium point The endemic equilibrium is given by P1 = (S∗f ,a f , I∗f ,S ∗ m,am, I∗m) where S∗f = µ(1−a f )− (γ f +µ)I∗f µ , S∗m = µ(1−am)− (γm +µ)I∗m µ , I∗m = µbm(1−a f )I∗f (γm +µ)(µ +bmI∗f ) , I∗f = µ(γm +µ)(R f Rm−1) bm(γm +µ +b f (1−am)) . This equilibrium corresponds to the state where the disease is endemic in both populations of men and women. Then we have the following result. Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 312 — #8 i i i i i i 312 ASSESSING THE IMPACT OF PROPHYLACTIC VACCINES ON HPV PREVALENCE Theorem 2. For RmR f > 1, or equivalently R0 > 1, there exists a unique endemic equilibrium P1. Global stability of the endemic equilibrium P1 is proved next using a suitable Lyapunov function. Theorem 3. Let R0 > 1, then the endemic equilibrium P1 of system (2.3) is globally asymptotically stable. Proof. Assume R0 > 1. To prove the global stability of the equilibrium P1 = (S∗f ,a f , I∗f ,S ∗ m,am, I∗m), we will show that the function W = d1(S f −S∗f −S∗f ln S f S∗f )+d2(Vf −V ∗f −V ∗f ln Vf V ∗f ) + d3(I f − I∗f − I∗f ln I f I∗f )+d4(Sm−S∗m−S∗m ln Sm S∗m ) + d5(Vm−V ∗m−V ∗m ln Vm V ∗m )+d6(Im− I∗m− I∗m ln Im I∗m ) with di > 0, i = 1, ...,6, is a Lyapunov function for the equilibrium P1, that is: i) W(X)≥ 0, where X = (S f ,Vf , I f ,Sm,Vm, Im) ∈Ω, and W(P1) = 0, ii) The orbital derivative Ẇ(X)< 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 Ẇ(X) = d1 ( 1− S∗f S f )( µ(1−a f )−b f S f Im−µS f ) +d2 ( 1− V ∗f V f )( µa f −µV f ) + d3 ( 1− I∗f I f )( b f S f Im− (γ f +µ)I f ) +d4 ( 1− S∗m Sm )( µ(1−am)−bmSmI f −µSm ) + d5 ( 1− V ∗m Vm ) (µam−µVm)+d6 ( 1− I∗m Im )( bmSmI f − (γm +µ)Im ) . From system (2.3) at equilibrium we obtain the following relations (1−a f )µ = b f S∗f S∗m +µS∗f a f = V ∗f γ f +µ = b f S∗f I∗m I∗f (1−am)µ = bmS∗mI∗f +µS∗m am = V ∗m γm +µ = bmS∗mI∗f I∗m . Substituting the parameters above in the derivative (3.5), and taking d1 = d3, d4 = d6, d4 = b f S∗f I∗m bmS∗mI∗f , Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 313 — #9 i i i i i i AZEVEDO, ESTEVA and FERREIRA 313 after several calculations and simplifications, we obtain that Ẇ can be written as Ẇ(X) = −d1µ (S f −S∗f ) 2 S f −d2µ (Vf −V ∗f ) 2 Vf −d4µ (Sm−S∗m) 2 Sm −d5µ (Vm−V ∗m) 2 Vm Vm −4d1b f S∗f I∗m [ 1 4 ( S∗f S f + S f ImI∗f S∗f I∗mI f + S∗m Sm + SmI f I∗m S∗mI∗f Im ) −1 ] . (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 ,Vf , I f ,Sm,Vm, Im) = P1. It follows that Ẇ(X)≤ 0, and Ẇ(X) = 0 if and only if X = P1. Therefore, trajectories in the interior of Ω approach P1. 4 NUMERICAL SIMULATIONS AND SENSITIVITY ANALYSIS In this section we present the numerical simulations of system (2.3) and a sensitivity analysis of R0 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 t=0) S f proportion of susceptible females at time t S f = 0.73 Vf proportion of vaccinated females at time t Vf = 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 Sm proportion of susceptible males at time t Sm = 0.64 Vm proportion of vaccinated males at time t Vm = 0.0 Im proportion of infected males at time t Im = 0.36 Zm proportion of recover males at time t Zm = 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, Sm = 0.64, Im = 0.36, and Vf = Vm = 0. In Figures 2 and 3 the disease prevails in female and male populations because R0 > 1, while in Figures 4 and 5, R0 < 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 Rm = 0.98 < 1; on the other hand, the disease can dies out from the population, even when R f > 1 or Rm > 1 as it is shown in Figure 4, where R f = 1.26. Figure 6 shows the partial rank correlation coefficient (PRCC) obtained from a global sensitiv- ity analysis [18]. The input parameters were sampled using the latin hypercube sampling with Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 314 — #10 i i i i i i 314 ASSESSING THE IMPACT OF PROPHYLACTIC VACCINES ON HPV PREVALENCE 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 100 200 300 400 500 600 P ro p o rt io n o f th e p o p u la ti o n Time (days) Sf Sm Vf Vm 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0 100 200 300 400 500 600 P ro p o rt io n o f th e p o p u la ti o n Time (days) If Im Figure 2: Temporal dynamics of the populations. The parameters are σ f 1×σ f 2 = 0.4, σm1× σm2 = 0.28, ε1× ε2 = 0.9, c f ×β f = 0.012 days−1, cm×βm = 0.006 days−1, µ = 3.65× 10−5 days−1, η f = ηm = 16.2 years, γ f = γm = 0.0028 days−1, and η = 6 months. The corresponding thresholds are R f = 2.19,Rm = 1.28, and R0 = 1.67 0 0.2 0.4 0.6 0.8 1 0 100 200 300 400 500 600 P ro p o rt io n o f th e p o p u la ti o n Time (days) Sf Sm Vf Vm 0 0.005 0.01 0.015 0.02 0.025 0.03 0 100 200 300 400 500 600 P ro p o rt io n o f th e p o p u la ti o n Time (days) If Im Figure 3: Temporal dynamics of the populations. The parameters are σ f 1×σ f 2 = 0.4, σm1× σm2 = 0.4, ε1× ε2 = 0.9, c f ×β f = 0.012 days−1, cm×βm = 0.0053 days−1, µ = 3.65× 10−5 days−1, η f = ηm = 16.2 years, γ f = γm = 0.0028 days−1, and η = 6 months. The corresponding thresholds are R f = 2.19,Rm = 0.98, and R0 = 1.46. Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 315 — #11 i i i i i i AZEVEDO, ESTEVA and FERREIRA 315 0 0.1 0.2 0.3 0.4 0.5 0.6 0 100 200 300 400 500 P ro p o rt io n o f th e p o p u la ti o n Time (days) Sf Sm Vf Vm 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 2 4 6 8 10 12 14 P ro p o rt io n o f th e p o p u la ti o n Time (days) If Im Figure 4: Temporal dynamics of the populations. The parameters are σ f 1×σ f 2 = 0.5, σm1× σm2 = 0.6, ε1× ε2 = 0.9, c f ×β f = 0.008 days−1, cm×βm = 0.0053 days−1, µ = 3.65× 10−5 days−1, η f = ηm = 16.2 years, γ f = γm = 0.0028 days−1, and η = 6 months. The corresponding thresholds are R f = 1.26,Rm = 0.7, and R0 = 0.94. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 100 200 300 400 500 P ro p o rt io n o f th e p o p u la ti o n Time (days) Sf Sm Vf Vm 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 2 4 6 8 10 P ro p o rt io n o f th e p o p u la ti o n Time (days) If Im Figure 5: Temporal dynamics of the populations. The parameters are σ f 1×σ f 2 = 0.7, σm1× σm2 = 0.6, ε1× ε2 = 0.9, c f ×β f = 0.008 days−1, cm×βm = 0.0053 days−1, µ = 3.65× 10−5 days−1, η f = ηm = 16.2 years, γ f = γm = 0.0028 days−1, and η = 6 months. The corresponding thresholds are R f = 0.85,Rm = 0.7, and R0 = 0.77. Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 316 — #12 i i i i i i 316 ASSESSING THE IMPACT OF PROPHYLACTIC VACCINES ON HPV PREVALENCE 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×10−5 - days−1 c f β 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 σ f 1×σ f 2 proportion of vaccinated females with 2-dose 0.40 0 - 1.0 - σm1×σ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 ε1× ε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 ε1ε2σ f1σ f2 = ε1ε2σm1σm2 = Unif(0,0.9) [1,16,19,22], c f β f =Unif(0.0073,0.0123) days−1 [4,14], cmβm=Unif(0.0053,0.0097) days−1 [4, 14], η f = ηm =Unif(14,20) years [19], and γ f = γm =Unif(0.0013,0.0055) days−1 [11, 21] where Unif means the uniform probability distribution. The other parameters are µ = 3.65× 10−5 days−1 and η = 6 months [9]. We run 5000 simula- tions that correspond to different input parameter sets, and the output chosed for the sensitivity analysis is R0. Figure 7 shows the disease prevalence versus the proportion of immunized females, σ f1σ f2 (as- suming that ε1ε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, σm1σm2 . Figures 8a-b show the disease prevalence versus the transmission rates c f β f and cmβm. The others parameters in the simulations are µ = 3.65×10−5, ε1× ε2 = 0.9, η f = ηm = 16.2 years, and η = 6 months. Each curve corresponds to a chosen recover rate, γ f , and it is assumed that γm = γ f . 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 Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 317 — #13 i i i i i i AZEVEDO, ESTEVA and FERREIRA 317 σf1σf2ε1ε2 c fβf τf γf cmβm dummy P R C C − 1 .0 − 0 .5 0 .0 0 .5 1 .0 Figure 6: Sensitivity analysis using R0 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. 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 I* f + I * m σf1 σf2 Rm=1.7, σm1σm2=0.0 Rm=1.4, σm1σm2=0.2 Rm=1.1, σm1σm2=0.4 Rm=0.8, σm1σm2=0.6 Rm=0.5, σm1σm2=0.8 Figure 7: Disease prevalence versus proportion of immunized girls. The parameters are σm1× σm2 = 0.6, ε1× ε2 = 0.9 c f × β f = 0.011 days−1, cm× βm = 0.006 days−1, µ = 3.65× 10−5 days−1, η f = ηm = 16.2 years, γ f = γm = 0.0028 days−1, and η = 6 months. Each curve corresponds to a given proportion of immunized boys. Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 318 — #14 i i i i i i 318 ASSESSING THE IMPACT OF PROPHYLACTIC VACCINES ON HPV PREVALENCE 0 0.005 0.01 0.015 0.02 0.025 0.03 0.008 0.01 0.012 I* f + I * m cfβf (a) 0 0.005 0.01 0.015 0.02 0.025 0.03 0.006 0.008 cmβm (b) 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) cmβm = 0.007 days−1 and in (b) c f β f = 0.0096 days−1. The other parameters are µ = 3.65× 10−5 days−1, ε1× ε2 = 0.92, σ f 1× σ f 2 = 0.4, σm1×σm2 = 0.28 η f = ηm = 16.2 years, and η = 6 months. 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 sum- marize the important parameters related to disease transmission. We found an expression for the basic reproductive number of the disease, R0, corresponding to the average of the basic repro- ductive numbers of female and male population R f , and Rm, respectively, which in turn depend Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 319 — #15 i i i i i i AZEVEDO, ESTEVA and FERREIRA 319 on the parameters that characterize virus-host dynamics and vaccination efficacy. Thus, if R0 < 1 the disease dies out; otherwise if R0 > 1 it persists. Persistence of the disease occurs when one of the populations (or both) is an infection source, i.e. Rm > 1 or R f > 1 (see Figures 1-4). Our findings show that control strategies that include even a small percentage of boys increase signif- icantly 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, ε1ε2, or vaccination compliance, σ fiσ fi i = f ,m, is more efficient than decreasing transmission rate c f β f or cmβm. In particular, the contribution of male and female transmission rates in R0, 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 R0; 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. 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 reprodu- tivo básico da doença, R0, é 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 R0 < 1. Além disso, tem um estado endêmico único quando R0 > 1 que é globalmente assintoticamente estável. Simulações numéricas foram feitas para comparar vários esquemas de vacinação diferentes. Os resul- tados 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. Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 320 — #16 i i i i i i 320 ASSESSING THE IMPACT OF PROPHYLACTIC VACCINES ON HPV PREVALENCE REFERENCES [1] Center for disease control and prevention. https://www.cdc.gov/nchs/products/databriefs/db280.htm. [2] M. Al-arydah & R.J. Smith. An age-structured model of human papilloma virus vaccination. Mathematics and Computers in Simulation, 82 (2011), 629–652. [3] 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. [4] A.N. Burchell, F. Coutleé, P. Tellier, J. Hanley & E.L. Franco. Genital Transmission of Human Papil- lomavirus in Recently Formed Heterosexual Couples. The Journal of Infectious Diseases, 204 (2011), 1723–1729. [5] 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. [6] O. Diekmann & J.A.P. Heesterbeek. “Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation”. Wiley (2000). [7] 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. [8] E.H. Elbasha & A.P.E. Galvani. Vaccination against multiple HPV types. Math. Biosci, 197 (2005), 88–117. [9] 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. [10] 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 Quadri- valent Human Papillomavirus Vaccine: A Systematic Review of 10 Years of Real-world Experience. Clinical Infectious Diseases, 63 (2016), 519–527. [11] 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. [12] J. Hale & H. Ko¸cak. “Dynamics and Bifurcations”. Springer-Verlag (1991). [13] D.M. Harpera & L.R. De Mars. HPV vaccines - A review of the first decade. Gynecologic Oncology, 146 (2017), 196–204. [14] B.Y. Hernandez, L.R. Wilkens, X. Zhu, P. Thompson, M. McDuffi, Y.B. Shvetsov, L.E. Kamem- oto, J. Killeen, L. Ning & M.T. Goodman. Transmission of Human Papillomavirus in Heterosexual Couples. Emerging Infectious Diseases, 14 (2008), 888–894. Tend. Mat. Apl. Comput., 20, N. 2 (2019) i i “A7-1286-6755-1-LE” — 2019/9/11 — 10:19 — page 321 — #17 i i i i i i AZEVEDO, ESTEVA and FERREIRA 321 [15] S.K. Kjaer, M. Nygard, J. Dillner, J.B. Marshall, D. Radley, M. Li, C. Munk, B.T. Hansen, L.G. Sig- urdardottir, 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. [16] 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. [17] 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. [18] 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. [19] 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. [20] 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. [21] M. Stanley. Immune responses to human papillomavirus. Vaccine, 24 (2006), DOI:10.1016/j.vaccine.2005.09.002. [22] M. Stanley. HPV - immune response to infection and vaccination. Infectious Agents and Cancer, 5 (2010), 19–25. [23] 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. Tend. Mat. Apl. Comput., 20, N. 2 (2019) Introduction Formulation of the model Analysis of the model Disease-free equilibrium and R0 Existence and stability of endemic equilibrium point Numerical Simulations and Sensitivity Analysis Discussion