Hindawi Abstract and Applied Analysis Volume 2019, Article ID 4567825, 13 pages https://doi.org/10.1155/2019/4567825 Research Article Optimal Control against the Human Papillomavirus: Protection versus Eradication of the Infection Fernando Saldaña ,1,2 Andrei Korobeinikov ,2,3 and Ignacio Barradas 1 1 Centro de Investigación en Matemáticas, 36240 Guanajuato, Mexico 2 Centre de Recerca Matemàtica, Campus de Bellaterra, 08193 Barcelona, Spain 3 Departament de Matemàtiques, Universitat Autònoma de Barcelona, Campus de Bellaterra, 08193 Barcelona, Spain Correspondence should be addressed to Fernando Saldaña; fernando.saldana@cimat.mx Received 4 March 2019; Accepted 9 April 2019; Published 9 May 2019 Academic Editor: Shengqiang Liu Copyright © 2019 Fernando Saldaña et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. We investigate the optimal vaccination and screening strategies to minimize human papillomavirus (HPV) associated morbidity and the interventions cost. We propose a two-sex compartmental model of HPV-infection with time-dependent controls (vaccination of adolescents, adults, and screening) which can act simultaneously. We formulate optimal control problems complementing our model with two different objective functionals. The first functional corresponds to the protection of the vulnerable group and the control problem consists of minimizing the cumulative level of infected females over a fixed time interval. The second functional aims to eliminate the infection, and, thus, the control problem consists of minimizing the total prevalence at the end of the time interval. We prove the existence of solutions for the control problems, characterize the optimal controls, and carry out numerical simulations using various initial conditions. The results and properties and drawbacks of the model are discussed. 1. Introduction majority of HPV-associated diseases [3]. In addition to the nonavalent and quadrivalent vaccines, there is also an FDA- Human papillomavirus (HPV) is the leading etiological fac- approved bivalent vaccine (Cervarix) that targets HPV types tor for the development of cervical cancer. According to 16 and 18. the World Health Organization (WHO) report, cervical The introduction of HPV vaccines is changing the epi- cancer is the fourth most frequent cancer in women with demiology of cervical cancer and other HPV-associated approximately 570,000 new cases reported in 2018 [1]. This diseases each year. Since HPV-associated morbidity and fact makes the design, implementation, and maintenance of mortality are usually considerably higher in females than effective preventing policies against cervical cancer a highly in males, in the majority of countries the primary target relevant public health problem. This consideration motivates group of HPV vaccination programs is adolescent school research on the design and development of cost-effective girls aged 9–14 [4]. However, considering that Gardasil-9 has vaccination and screening programs against HPV infection. been licensed for a broader age range in both males and On October 5th of 2018, the US Food and Drug Admin- females, it is essential to analyze the potential health benefits istration (FDA) approved the nonavalent vaccine Gardasil-9 of extending vaccination to males, as well as older individuals, against HPV infection for use in women and men aged 27 versus a higher economic cost of such programs. to 45 [2]. Earlier, the FDA approved the vaccine for persons To examine the dynamics of HPV spread and control, between the ages of 9 to 26. Gardasil-9 protects against numerous mathematical models have been proposed [5– HPV types 6, 11, 16, and 18 (that are already targeted by the 15]. While these and other recent studies analyzed different quadrivalent HPV vaccine Gardasil) and against the next five strategies for HPV control and brought important insight into most common oncogenic viral types, namely, HPV 31, 33, the problem, the majority of them did not take into account 45, 52, and 58. These nine HPV types are responsible for the the optimality of these interventions. The study of optimality, 2 Abstract and Applied Analysis however, is of particular relevance in the allocation of limited according to gender and infection status. We subdivide the public health resources. In this context, the application of total female population 𝑁𝑓 (𝑡) into four mutually exclusive optimal control theory has proven to be an important tool compartments: unvaccinated susceptibles 𝑆𝑓 (𝑡), vaccinated for evaluating and optimizing various detection, prevention, susceptibles 𝑉𝑓 (𝑡), and infectious females unaware 𝑈𝑓 (𝑡) and therapy, vaccination, and other intervention programs [16]. aware 𝐼𝑓 (𝑡) of their infection; thus So far, to the best of authors’ knowledge, for the case of HPV infection only a few studies have been done using 𝑁𝑓 (𝑡) = 𝑆𝑓 (𝑡) + 𝑉𝑓 (𝑡) + 𝑈𝑓 (𝑡) + 𝐼𝑓 (𝑡) . (1) the optimal control framework to assess the impact of vaccination programs [17, 18]. Brown and White [17] explored The total male population 𝑁𝑚 (𝑡) consists of three compart- how to target vaccination in the UK. This study highlights ments: unvaccinated susceptibles 𝑆𝑚 (𝑡), vaccinated suscepti- the importance of including a catch-up vaccination policy bles 𝑉𝑚 (𝑡), and infectious males 𝐼𝑚 (𝑡), so that in order to control the spread of the infection. In a recently published paper, Malik et al. [18] analyzed the optimal 𝑁𝑚 = 𝑆𝑚 (𝑡) + 𝑉𝑚 (𝑡) + 𝐼𝑚 (𝑡) . (2) control strategies for a vaccination program administering We assume that individuals enter the sexually active pop- the bivalent, quadrivalent, and nonavalent HPV vaccines to ulation as singles at a constant rate Λ with a sex ratio of 50:50. the female population. In particular, they explored scenarios Females and males leave the population by death or ceasing where the three vaccines are used simultaneously compared sexual activity at per capita rates 𝜇𝑓 and 𝜇𝑚 , respectively. with the case where the bivalent and quadrivalent vaccines In this study, we distinguish between vaccinations prior and were initially used and then, during the program, one or both after sexual initiation. In particular, a fraction 𝑤1 of females of them are replaced by the nonavalent vaccine. and a fraction 𝑤2 of males are vaccinated before they enter Since there is evidence that switching from the quadriva- the sexually active class and thus are recruited into their lent to the nonavalent is likely to be cost-saving [3], in this vaccinated compartment. Moreover, the susceptible sexually work we focus on the application of the nonavalent vaccine. active females and males are vaccinated at rates 𝑢1 and 𝑢2 The aim is to use optimal control theory to gain insight per unit of time, respectively. The vaccine reduces the force into the best combination of vaccination and screening to of infection by a factor 0 ≤ 𝜖 ≤ 1; thus, the vaccine is reduce the spread of HPV infection, as well as the cost of completely effective when 𝜖 = 0 and ineffective when 𝜖 = the intervention strategy. To address this issue, we construct a 1. Therefore 1 − 𝜖 is the vaccine effectiveness. The vaccine- mathematical model that allows vaccination prior to and after induced immunity wanes at a rate 𝜃. Both males and females sexual initiation in both males and females. Moreover, we can clear the infection naturally with per capita rates 𝛾𝑚 and also include into the model a screening program that allows 𝛾𝑓 , respectively. No permanent immunity is assumed after detection of HPV-associated diseases. that. In general, public health interventions have two basic Susceptible females become infected at a rate 𝛽𝑚 fol- objectives: protection of population groups that are at consid- lowing effective contact with an infectious male. After the erably higher risk of contracting or developing severe disease acquisition of HPV infection, a fraction 𝑝 of the infected and total elimination of the infection when it is potentially females may develop symptoms and become aware of their possible. Accordingly, in this paper, we consider both of these infection entering the 𝐼𝑓 (𝑡) class. The remaining infected objectives formulating two optimal control problems. Our females do not realize their infection and move to the ultimate goal is to explore where control policies aimed at the unaware class 𝑈𝑓 (𝑡). However, the screening program may protection of the vulnerable group and at the eradication of allow the unaware infected females to detect their infection the infection differ. at a rate 𝛼. The organization of this paper is as follows: The next Susceptible males can be infected by unaware infected section is devoted to formulation and analysis of an epidemic females 𝑈𝑓 (𝑡) at a rate 𝛽𝑓 and by aware infected females model for the transmission dynamics of HPV infection in a heterosexual population. In Section 3, we propose optimal 𝐼𝑓 (𝑡) at a rate 𝛽̃𝑓 . We assume that 0 < 𝛽̃𝑓 < 𝛽𝑓 since a control problems incorporating two objective functionals female conscious of her infection can take precautions like and simultaneous application of several controls into our condom use to reduce the possibility of transmission. Since model. In Section 4, we prove the existence of solutions and HPV infection is transmitted (in most cases) sexually, we characterization of the optimal controls via the Pontryagin model the transmission via the standard incidence [22]. maximum principle. Section 5 contains the results of numer- The model resulting from these assumptions is governed ical simulations of the model and the profiles of the optimal by the following system of differential equations: controls. The last section contains the conclusions and a Λ 𝐼 discussion of the obtained results. 𝑆𝑓̇ = (1 − 𝑤1 ) − 𝛽𝑚 𝑆𝑓 𝑚 − 𝑢1 𝑆𝑓 + 𝛾𝑓 𝑈𝑓 + 𝛾𝑓 𝐼𝑓 2 𝑁𝑚 2. Model Formulation and Analysis + 𝜃𝑉𝑓 − 𝜇𝑓 𝑆𝑓 , In this study, we propose and study a compartmental 𝐼 model for the transmission dynamics of the HPV types 𝑈̇ 𝑓 = (𝑆𝑓 + 𝜖𝑉𝑓 ) (1 − 𝑝) 𝛽𝑚 𝑚 − (𝛾𝑓 + 𝛼) 𝑈𝑓 𝑁𝑚 targeted by the nonavalent vaccine Gardasil-9. The model classifies the hosts’ population at time 𝑡 denoted by 𝑁(𝑡) − 𝜇𝑓 𝑈𝑓 , Abstract and Applied Analysis 3 𝐼𝑚 𝐼𝑓̇ = (𝑆𝑓 + 𝜖𝑉𝑓 ) 𝑝𝛽𝑚 + 𝛼𝑈𝑓 − 𝛾𝑓 𝐼𝑓 − 𝜇𝑓 𝐼𝑓 , Please note that system (5) always has a unique disease- 𝑁𝑚 free equilibrium 𝐸0 = (𝑆∗𝑓 , 𝑈𝑓∗ , 𝐼𝑓∗ , 𝑉𝑓∗ , 𝑆∗𝑚 , 𝐼𝑚 ∗ , 𝑉𝑚∗ )𝑇 given by Λ 𝐼 𝑉𝑓̇ = 𝑤1 + 𝑢1 𝑆𝑓 − 𝜖𝛽𝑚 𝑉𝑓 𝑚 − (𝜇𝑓 + 𝜃) 𝑉𝑓 , (1 − 𝑤1 ) 𝜇𝑓 + 𝜃 𝑤1 𝜇𝑓 + 𝑢1 2 𝑁𝑚 𝐸0 = ( , 0, 0, , 𝑢1 + 𝜇𝑓 + 𝜃 𝑢1 + 𝜇𝑓 + 𝜃 Λ 𝑈𝑓 𝐼𝑓 (6) ̇ = (1 − 𝑤2 ) 𝑆𝑚 − (𝛽𝑓 + 𝛽̃𝑓 ) 𝑆 − 𝑢2 𝑆𝑚 𝑇 2 𝑁𝑓 𝑁𝑓 𝑚 (1 − 𝑤2 ) 𝜇𝑚 + 𝜃 𝑤 𝜇 + 𝑢2 , 0, 2 𝑚 ) . 𝑢2 + 𝜇𝑚 + 𝜃 𝑢2 + 𝜇𝑚 + 𝜃 + 𝛾𝑚 𝐼𝑚 + 𝜃𝑉𝑚 − 𝜇𝑚 𝑆𝑚 , 𝑈𝑓 𝐼𝑓 In order to compute the basic reproduction number, we 𝐼𝑚̇ = (𝛽𝑓 + 𝛽̃𝑓 ) (𝑆𝑚 + 𝜖𝑉𝑚 ) − (𝛾𝑚 + 𝜇𝑚 ) 𝐼𝑚 , need the following components. During his infection period, 𝑁𝑓 𝑁𝑓 (𝛾𝑚 + 𝜇𝑚 )−1 , an infectious male produces on average 𝛽𝑚 [𝑆∗𝑓 + 𝑈 𝐼 𝜖𝑉𝑓∗ ] infections in females. Therefore, ̇ = 𝑤2 Λ − (𝛽𝑓 𝑓 + 𝛽̃𝑓 𝑓 ) 𝜖𝑉𝑚 + 𝑢2 𝑆𝑚 𝑉𝑚 2 𝑁𝑓 𝑁𝑓 𝛽𝑚 (𝑆∗𝑓 + 𝜖𝑉𝑓∗ ) 𝑇𝑚𝑓 = (7) − (𝜇𝑚 + 𝜃) 𝑉𝑚 , 𝛾𝑚 + 𝜇𝑚 (3) is the infection transfer from infectious males to females. An infectious female aware of her infection produces on average where all the parameters are assumed to be nonnegative. 𝛽̃𝑓 (𝑆∗𝑚 + 𝜖𝑉𝑚∗ ) infections during her infectious period (𝛾𝑓 + We are interested in the dynamics of model (3) over a finite time interval [0, 𝑇]. Therefore, we assume that the 𝜇𝑓 )−1 . As a result, population for both sexes is constant with values 𝑁𝑓 = Λ/2𝜇𝑓 and 𝑁𝑚 = Λ/2𝜇𝑚 . This assumption is reasonable since the 𝛽̃𝑓 (𝑆∗𝑚 + 𝜖𝑉𝑚∗ ) epidemic process is considerably faster than the demographic 𝑇𝐼𝑚𝑓 = (8) 𝛾𝑓 + 𝜇𝑓 one over a short period of time. Under the constant population assumption, we can is the infection transfer from the aware infectious females to express system (3) in terms of the proportions males. Analogously, 𝑆𝑓 𝑈𝑓 𝐼𝑓 𝛽𝑓 (𝑆∗𝑚 + 𝜖𝑉𝑚∗ ) 𝑆𝑓 = , 𝑈𝑓 = , 𝐼𝑓 = , 𝑇𝑈𝑚𝑓 = (9) 𝑁𝑓 𝑁𝑓 𝑁𝑓 𝛾𝑓 + 𝛼 + 𝜇𝑓 (4) 𝑉𝑓 𝑆 𝐼 is the infection transfer from the unaware infectious females 𝑉𝑓 = , 𝑆𝑚 = 𝑚 , 𝐼𝑚 = 𝑚 . to males. Moreover, a fraction 𝛼/(𝛾𝑓 + 𝛼 + 𝜇𝑓 ) of unaware 𝑁𝑓 𝑁𝑚 𝑁𝑚 infected females becomes aware of the infection via screen- ing; thus Omitting the bars above the variables for convenience, our normalized system becomes 𝛼 𝑇𝑓𝑚 = 𝑝𝑇𝐼𝑚𝑓 + (1 − 𝑝) (𝑇𝑈𝑚𝑓 + 𝑇𝑚 ) (10) 𝛾𝑓 + 𝛼 + 𝜇𝑓 𝐼𝑓 𝑆𝑓̇ = (1 − 𝑤1 ) 𝜇𝑓 − 𝛽𝑚 𝑆𝑓 𝐼𝑚 + 𝛾𝑓 𝑈𝑓 + 𝛾𝑓 𝐼𝑓 + 𝜃𝑉𝑓 is the infection transfer from females to males. − (𝑢1 + 𝜇𝑓 ) 𝑆𝑓 , Using the next-generation matrix 𝑈̇ 𝑓 = (𝑆𝑓 + 𝜖𝑉𝑓 ) (1 − 𝑝) 𝛽𝑚 𝐼𝑚 − (𝛾𝑓 + 𝛼) 𝑈𝑓 − 𝜇𝑓 𝑈𝑓 , 0 0 (1 − 𝑝) 𝑇𝑚𝑓 [ ] [ 𝑓 ] 𝐼𝑓̇ = (𝑆𝑓 + 𝜖𝑉𝑓 ) 𝑝𝛽𝑚 𝐼𝑚 + 𝛼𝑈𝑓 − 𝛾𝑓 𝐼𝑓 − 𝜇𝑓 𝐼𝑓 , [ K=[ 0 0 𝑝𝑇 𝑚 ], (11) ] [ 𝑚 𝛼 ] 𝑇𝑈𝑓 + 𝑇𝐼𝑚𝑓 𝑇𝐼𝑚𝑓 0 𝑉𝑓̇ = 𝑤1 𝜇𝑓 + 𝑢1 𝑆𝑓 − 𝜖𝛽𝑚 𝑉𝑓 𝐼𝑚 − (𝜇𝑓 + 𝜃) 𝑉𝑓 , 𝛾 + 𝛼 + 𝜇 (5) [ 𝑓 𝑓 ] ̇ = (1 − 𝑤2 ) 𝜇𝑚 − (𝛽𝑓 𝑈𝑓 + 𝛽̃𝑓 𝐼𝑓 ) 𝑆𝑚 + 𝛾𝑚 𝐼𝑚 + 𝜃𝑉𝑚 𝑆𝑚 it can be shown that the basic reproduction number R0 is the geometric mean of the terms 𝑇𝑚𝑓 and 𝑇𝑓𝑚 ; that is, − (𝑢2 + 𝜇𝑚 ) 𝑆𝑚 , 𝐼𝑚̇ = (𝛽𝑓 𝑈𝑓 + 𝛽̃𝑓 𝐼𝑓 ) (𝑆𝑚 + 𝜖𝑉𝑚 ) − (𝛾𝑚 + 𝜇𝑚 ) 𝐼𝑚 , 𝑓 R0 = √ 𝑇𝑚 ⋅ 𝑇𝑓𝑚 . (12) ̇ = 𝑤2 𝜇𝑚 − (𝛽𝑓 𝑈𝑓 + 𝛽̃𝑓 𝐼𝑓 ) 𝜖𝑉𝑚 + 𝑢2 𝑆𝑚 𝑉𝑚 As a consequence of Theorem 2 in [23], the disease-free equilibrium 𝐸0 is locally asymptotically stable if R0 < 1 and − (𝜇𝑚 + 𝜃) 𝑉𝑚 . unstable if R0 > 1. 4 Abstract and Applied Analysis 3. The Optimal Control Problems minimizing the number of cervical cancer cases and related diseases in females over the time interval [0, 𝑇]. For this case, Five controls are possible in this modeling framework; hence, the objective functional is we formulate optimal control problems incorporating these 𝑇 five controls, namely: 1 𝐽 = ∫ (𝐴 0 (𝑈𝑓 + 𝐼𝑓 ) + [𝐴 (𝑤2 + 𝑤22 ) 0 2 1 1 (i) 𝑤1 (⋅) and 𝑤2 (⋅) are the controls representing the (16) fraction of females and males, respectively, that are + 𝐴 2 (𝑢21 + 𝑢22 ) 2 + 𝐴 3 𝛼 ]) 𝑑𝑡. vaccinated prior to sexual initiation (ii) 𝑢1 (⋅) and 𝑢2 (⋅) are the controls representing the vacci- Secondly, having the aim of the ultimate elimination of the nation rates of susceptible sexually active females and infection in mind, we consider the problem of minimizing the males, respectively total prevalence 𝑃(𝑡) = 𝑈𝑓 (𝑡) + 𝐼𝑓 (𝑡) + 𝐼𝑚 (𝑡) at the terminal (iii) 𝛼(⋅) is the control representing the screening rate for time moment 𝑇. In this context, the objective functional takes females the following form: Because resources of a health system are limited, we assume ̃0 𝑃 (𝑇) + 1 𝐽̃ = 𝐴 that all the controls are subject to constraints 2 𝑇 (17) 0 ≤ 𝑤1 (𝑡) , 𝑤2 (𝑡) ≤ 1, ⋅ ∫ (𝐴 1 (𝑤12 + 𝑤22 ) + 𝐴 2 (𝑢21 + 𝑢22 ) + 𝐴 3 𝛼2 ) 𝑑𝑡. 0 0 ≤ 𝑢1 (𝑡) , 𝑢2 (𝑡) ≤ 𝑢𝑚𝑎𝑥 , (13) In these functionals, the parameter 𝐴 0 balances the proba- bility of developing cervical cancer and its treatment cost, 0 ≤ 𝛼 (𝑡) ≤ 𝛼𝑚𝑎𝑥 . ̃0 reflects the cost of the long-term presence of HPV whereas 𝐴 With these controls and taking into consideration the con- infection in the population. The meaning of the parameters stant population size assumption, we can write model (5) as 𝐴 1 , 𝐴 2 , and 𝐴 3 in both functionals is the same. Specifically, the following system of five equations: the weight parameter 𝐴 1 represents the cost associated with the vaccination of school girls and boys prior to their sexual 𝑈̇ 𝑓 = ((1 − 𝑈𝑓 − 𝐼𝑓 − 𝑉𝑓 ) + 𝜖𝑉𝑓 ) (1 − 𝑝) 𝛽𝑚 𝐼𝑚 initiation. The weight parameter 𝐴 2 represents the cost of administering the vaccine in sexually active individuals. − (𝛾𝑓 + 𝛼 (𝑡) + 𝜇𝑓 ) 𝑈𝑓 , Finally, 𝐴 3 represents the cost of screening and follow-up after an abnormal Pap test. 𝐼𝑓̇ = ((1 − 𝑈𝑓 − 𝐼𝑓 − 𝑉𝑓 ) + 𝜖𝑉𝑓 ) 𝑝𝛽𝑚 𝐼𝑚 + 𝛼 (𝑡) 𝑈𝑓 We are assuming that the cost of applying health poli- cies increases as a nonlinear function of the controls. The − (𝛾𝑓 + 𝜇𝑓 ) 𝐼𝑓 , quadratic terms penalize high levels of control administra- tion in comparison with the cost of low levels. Due to its 𝑉𝑓̇ = 𝑤1 (𝑡) 𝜇𝑓 + 𝑢1 (𝑡) (1 − 𝑈𝑓 − 𝐼𝑓 − 𝑉𝑓 ) mathematical convenience, the sum of the weighted squares (14) of the controls is probably the most common form of the − 𝜖𝛽𝑚 𝑉𝑓 𝐼𝑚 − (𝜇𝑓 + 𝜃) 𝑉𝑓 , objective functional in the literature: for functionals of this type, it is possible to reduce an optimal control problem 𝐼𝑚̇ = (𝛽𝑓 𝑈𝑓 + 𝛽̃𝑓 𝐼𝑓 ) ((1 − 𝐼𝑚 − 𝑉𝑚 ) + 𝜖𝑉𝑚 ) to a two-point boundary value problem. Furthermore, for such functionals, the optimal controls can be obtained as − (𝛾𝑚 + 𝜇𝑚 ) 𝐼𝑚 , explicit functions of the state and adjoint variables via the Pontryagin maximum principle. Although the financial cost ̇ = 𝑤2 (𝑡) 𝜇𝑚 − (𝛽𝑓 𝑈𝑓 + 𝛽̃𝑓 𝐼𝑓 ) 𝜖𝑉𝑚 𝑉𝑚 of public health interventions to prevent epidemics most likely grows linearly with the magnitude of the corresponding + 𝑢2 (𝑡) (1 − 𝐼𝑚 − 𝑉𝑚 ) − (𝜇𝑚 + 𝜃) 𝑉𝑚 , control [24], we choose objective functionals with a quadratic dependence on the controls following the current trend of the complemented by initial conditions literature [17, 18, 25–27]. 0 ≤ 𝑈𝑓 (0) , 𝐼𝑓 (0) , 𝑉𝑓 (0) , 𝐼𝑚 (0) , 𝑉𝑚 (0) ≤ 1. (15) 4. The Pontryagin Maximum Principle The set of admissible controls 𝐷(𝑇) is the set of Lebesgue measurable functions (𝑤1 (⋅), 𝑤2 (⋅), 𝑢1 (⋅), 𝑢2 (⋅), 𝛼(⋅))𝑇 which For control model (14) and the set of admissible controls for almost all 𝑡 ∈ [0, 𝑇] satisfy constraints (13). 𝐷(𝑇), we can formulate two different optimal control prob- We seek the optimal controls that reduce the incidence lems using the objective functionals (16) and (17). Here, we and prevalence of the infection and also the cost of using a prove the existence of solutions for these optimal control public health intervention over 𝑇 years. In the case where problems and we also obtain the optimality system in order the elimination of HPV infection is not possible due to social to find a numerical approximation of the optimal controls. or economic constraints, the public health system can focus, Let us introduce the states 𝑥 = (𝑈𝑓 , 𝐼𝑓 , 𝑉𝑓 , 𝐼𝑚 , 𝑉𝑚 )𝑇 instead, on the protection of the vulnerable group, that is, on and the controls 𝑐 = (𝑤1 , 𝑤2 , 𝑢1 , 𝑢2 , 𝛼)𝑇 and denote the Abstract and Applied Analysis 5 right-hand side of the control model (14) as the vector Our next aim is to obtain the so-called optimality system function 𝐹(𝑡, 𝑥, 𝑐). For the objective functional (16), the in order to find numerically optimal controls 𝑤1∗ (𝑡), 𝑤2∗ (𝑡), optimal control problem of the HPV epidemic model is 𝑢∗1 (𝑡), 𝑢∗2 (𝑡), 𝛼∗ (𝑡). We define the Hamiltonian 𝑇 𝐻 = 𝐿 (𝑡, 𝑥, 𝑐) + 𝜓1 𝑈̇ 𝑓 + 𝜓2 𝐼𝑓̇ + 𝜓3 𝑉𝑓̇ + 𝜓4 𝐼𝑚̇ min 𝐽 (𝑐) = min ∫ 𝐿 (𝑡, 𝑥 (𝑡) , 𝑐 (𝑡)) (18) (20) 𝑐∈𝐷(𝑇) 𝑐∈𝐷(𝑇) 0 ̇ . + 𝜓5 𝑉𝑚 where By the Pontryagin maximum principle [28], given the optimal control 𝑐∗ (⋅) that minimizes the Hamiltonian (20) for almost 𝐿 (𝑡, 𝑥, 𝑐) all 𝑡 ∈ [0, 𝑇] and the corresponding optimal solutions 𝑥(𝑡) = = 𝐴 0 (𝑈𝑓 + 𝐼𝑓 ) (𝑈𝑓 (𝑡), 𝐼𝑓 (𝑡), 𝑉𝑓 (𝑡), 𝐼𝑚 (𝑡), 𝑉𝑚 (𝑡))𝑇, there are adjoint variables (19) 𝜓𝑖 that satisfy the system 1 + (𝐴 1 (𝑤12 + 𝑤22 ) + 𝐴 2 (𝑢21 + 𝑢22 ) + 𝐴 3 𝛼2 ) 𝜓̇1 = −𝐴 0 + [(1 − 𝑝) 𝛽𝑚 𝐼𝑚 + 𝛾𝑓 + 𝛼∗ + 𝜇𝑓 ] 𝜓1 2 subject to the dynamics of the control model (14) with given + (𝑝𝛽𝑚 𝐼𝑚 − 𝛼∗ ) 𝜓2 + 𝑢∗1 𝜓3 initial conditions. We call a pair of states and controls (𝑥, 𝑐) satisfying both − 𝛽𝑓 [(1 − 𝐼𝑚 − 𝑉𝑚 ) + 𝜖𝑉𝑚 ] 𝜓4 + 𝛽𝑓 𝜖𝑉𝑚 𝜓5 , (14) and 𝑐 ∈ 𝐷(𝑇) a feasible pair. Theorem 4.1 in [28, Chapter III] ensures the existence of an optimal control and 𝜓̇2 = −𝐴 0 + (1 − 𝑝) 𝛽𝑚 𝐼𝑚 𝜓1 + (𝑝𝛽𝑚 𝐼𝑚 + 𝛾𝑓 + 𝜇𝑓 ) 𝜓2 the corresponding solution for this problem. In particular, this existence theorem states that the following conditions are + 𝑢∗1 𝜓3 − 𝛽̃𝑓 [(1 − 𝐼𝑚 − 𝑉𝑚 ) + 𝜖𝑉𝑚 ] 𝜓4 sufficient to guarantee the existence of an optimal control for (14): + 𝛽̃𝑓 𝜖𝑉𝑚 𝜓5 , (H1) 𝐹 is continuous, and there exist positive constants 𝐾1 𝜓3̇ = [(1 − 𝑝) 𝜓1 + 𝑝𝜓2 ] (1 − 𝜖) 𝛽𝑚 𝐼𝑚 (21) and 𝐾2 such that + (𝜖𝛽𝑚 𝐼𝑚 + 𝑢∗1 + 𝜇𝑓 + 𝜃) 𝜓3 , (a) |𝐹(𝑡, 𝑥, 𝑐)| ≤ 𝐾1 (1 + |𝑥| + |𝑐|) 𝜓̇4 = −𝛽𝑚 ((1 − 𝑈𝑓 − 𝐼𝑓 − 𝑉𝑓 ) + 𝜖𝑉𝑓 ) ̃ 𝑐) − 𝐹(𝑡, 𝑥, 𝑐)| ≤ 𝐾2 |𝑥̃ − 𝑥|(1 + |𝑐|) (b) |𝐹(𝑡, 𝑥, hold for all 𝑡 ∈ [0, 𝑇]. Moreover, 𝐹 can be written ⋅ [(1 − 𝑝) 𝜓1 + 𝑝𝜓2 ] + 𝜖𝛽𝑚 𝑉𝑓 𝜓3 as 𝐹(𝑡, 𝑥, 𝑐) = 𝜙(𝑡, 𝑥) + 𝜓(𝑡, 𝑥)𝑐 + (𝛽𝑓 𝑈𝑓 + 𝛽̃𝑓 𝐼𝑓 + 𝛾𝑚 + 𝜇𝑚 ) 𝜓4 + 𝑢∗2 𝜓5 , (H2) The set of admissible controls 𝐷(𝑇) is closed and convex. Moreover, there is at least one feasible pair 𝜓̇5 = (𝛽𝑓 𝑈𝑓 + 𝛽̃𝑓 𝐼𝑓 ) ((1 − 𝜖) 𝜓4 + 𝜖𝜓5 ) (𝑥(𝑡), 𝑐(𝑡)) satisfying both (14) and (13) + (𝑢∗2 + 𝜃 + 𝜇𝑚 ) 𝜓5 , (H3) 𝐿(𝑡, 𝑥, ⋅) is convex on 𝐷(𝑇), and 𝐿(𝑡, 𝑥, 𝑐) ≥ 𝐾3 |𝑐|𝛽 − 𝐾4 , 𝐾3 > 0, 𝐾4 ∈ R, 𝛽 > 1 with the corresponding conditions 𝜓𝑖 (𝑇) = 0 for 𝑖 = 1, . . . , 5. From the optimality conditions, we obtain the following Theorem 1. Consider the optimal control problem with control characterization for the optimal controls: model (14) and cost functional (16). Then there exist 𝑐∗ = 𝜇𝑓 (𝑤1∗ , 𝑤2∗ , 𝑢∗1 , 𝑢∗2 , 𝛼∗ ) ∈ 𝐷(𝑇) such that min 𝐽(𝑐)𝑐∈𝐷(𝑇) = 𝐽(𝑐∗ ). 𝑤1∗ (𝑡) = min {1, max {0, − 𝐴1 𝜓3 (𝑡)}} , Proof. Due to the a priori boundedness of the solutions 𝜇𝑚 𝑤2∗ (𝑡) = min {1, max {0, − 𝜓 (𝑡)}} , of the control model (14) and the differentiability of the 𝐴1 5 function 𝐹, it follows that (a) and (b) in (H1) are ensured 𝑢∗1 (𝑡) by suitable bounds on the partial derivatives of 𝐹 and on 1 𝐹(𝑡, 0, 0). Moreover, the state equations are linear with respect = min {𝑢𝑚𝑎𝑥 , max {0, − 𝐴2 (1 − 𝑈𝑓 (𝑡) − 𝐼𝑓 (𝑡) − 𝑉𝑓 (𝑡)) 𝜓3 (𝑡)}} , to the controls 𝑐, and thus 𝐹(𝑡, 𝑥, 𝑐) = 𝜙(𝑡, 𝑥) + 𝜓(𝑡, 𝑥)𝑐. The (22) existence of a feasible pair is guaranteed by the Caratheodory 𝑢∗2 (𝑡) theorem [29, pp. 182] for initial value problems. Moreover, 1 = min {𝑢𝑚𝑎𝑥 , max {0, − (1 − 𝐼𝑚 (𝑡) − 𝑉𝑚 (𝑡)) 𝜓5 (𝑡)}} , for bounded controls on a finite time interval, 𝐷(𝑇) is clearly 𝐴2 closed and convex, and, hence, (H2) is satisfied as well. The 𝛼∗ (𝑡) integrand (19) of the cost functional is positive and quadratic with respect to the controls; therefore, 𝐿(𝑡, 𝑥, ⋅) is convex on 1 = min {𝛼𝑚𝑎𝑥 , max {0, (𝜓 (𝑡) − 𝜓2 (𝑡)) 𝑈𝑓 (𝑡)}} . 𝐴3 1 𝐷(𝑇). Furthermore, defining 𝐾3 = min{𝐴 1 , 𝐴 2 , 𝐴 3 }, we have 𝐿(𝑡, 𝑥, 𝑐) ≥ (1/2)𝐾3 |𝑐|2 , thus verifying (H3) and completing The proof of the existence of solutions for the optimal the proof. control problem with the objective functional (17) follows the 6 Abstract and Applied Analysis Table 1: Parameters values of the control model (14) used in the simulations. Parameter Range Mean Value Units Source 1−𝜖 [0.8, 1] 0.9 adimensional [3] 1/𝜃 [5, 15] 10 year [19] 𝛽𝑚 [0.05, 5] 4.0 year−1 [3] 𝛽𝑓 [0.05, 5] 4.0 year−1 [3] 𝛽̃𝑓 [0.025, 2.5] 2.0 year−1 Estimated 1/𝛾𝑓 [0.83, 2] 1.3 year [20] 1/𝛾𝑚 [0.33, 1.2] 0.6 year [21] 𝑝 [0, 1] 0.2 adimensional Estimated 1/𝜇𝑓 [15, 50] 30 year Estimated 1/𝜇𝑚 [15, 50] 30 year Estimated same reasoning as in Theorem 1 and is, therefore, omitted. The that the aware infected females take measures, such as technical details for the computation of the optimality system condom use, that decrease the probability of infection and characterization of the optimal control associated with transfer. Thus 𝛽̃𝑓 ∈ [0.025, 2.5] 𝑦𝑒𝑎𝑟−1 . objective functional (17) are provided in the appendix. (iv) Females’ infectious period, 1/𝛾𝑓 . The mean duration of a new genital HPV infection in women is 14.8 months 5. Results for oncogenic types and 11.1 months for nononcogenic We start by exploring the control model (14) numerically. types [20]. Since the majority of the types targeted by The optimal control solutions can be obtained solving the Gardasil-9 are oncogenic, we assume that the female optimality system. In order to obtain approximations of infectious period is between 9 months and 2 years; our optimality systems, we use the forward-backward sweep thus 1/𝛾𝑓 ∈ [0.83, 2] 𝑦𝑒𝑎𝑟. method (FBSM) from [30]. This iterative scheme begins with (v) Male’s infectious period, 1/𝛾𝑚 . For men, most of the an initial guess on the controls; then, the state equations are HPV infections clear in less than 12 months. Studies solved forward in time. After that, since the transversality (see [21] and the references therein) report average conditions for our problems are posed at the final time 𝑇, the clearance time of 5.9 to 7.5 months. We set 1/𝛾𝑚 ∈ system for the adjoint variables is solved backward in time. [0.33, 1.2] year. The controls are updated using a convex combination of the (vi) Fraction of females that become aware of their infec- previous controls and the new control obtained substituting tion, 𝑝. Since 𝑝 is a fraction, 𝑝 ∈ [0, 1]. the states and adjoints into its characterization. This process is repeated until convergence criteria are satisfied. For solving (vii) Periods of sexual activity, 1/𝜇𝑓 , and 1/𝜇𝑚 . We assume the state and adjoint systems, we use a fourth order Runge- that the average duration of the sexual activity period Kutta scheme. for both, males and females, is between 15 and 50 years; therefore 𝜇𝑓−1 , 𝜇𝑚 −1 ∈ [15, 50] years. 5.1. Model Parameters. Parameters values, units, and source Estimates of weight parameters associated with the ob- of estimation are summarized in Table 1. We use these values jective functionals exhibit considerable variability in the unless otherwise is stated. Next, we outline their selection: literature. Thus, in the numerical simulations, several com- (i) Vaccine efficacy, 1 − 𝜖. Although clinical trials are binations of weight parameters are explored. still underway, there is evidence indicating efficacy of around 90% for the nonavalent HPV vaccine 5.2. Uncontrolled System and Constant Control. In order to Gardasil-9 [3]. We set 1 − 𝜖 ∈ [0.8, 1]. demonstrate the importance of the optimal controls, we firstly show the behavior of the model in the absence of control and (ii) Vaccine protection, 1/𝜃. The exact duration of vaccine under constant controls. The results are shown in Figure 1. It protection is unknown, but clinical trials [19] have is easy to see that, without control, the system converges to shown sustained efficacy for more than 5 years. We an endemic equilibrium state (see Figure 1(a)). With constant assume that the protection is from 5 to 15 years; thus controls (see Figure 1(b)) the infectious cases are reduced but 1/𝜃 ∈ [5, 15] years. the infection is not eradicated for this time interval. Besides, (iii) Transmission rates, 𝛽𝑚 , 𝛽𝑓 , and 𝛽̃𝑓 . Transmission rates the number of vaccinated individuals tends to a very high are the product of the average number of sexual con- proportion. tacts per unit of time and the probability of infection transfer per contact. We assume that both genders 5.3. Optimal Vaccination and Screening. For the numerical have between 1 and 5 sexual partners per year and per- computation of the optimal controls, we should complement act transmission probability range 5−100% [3]. We set the fixed parameters in Table 1 with values for the weights in 𝛽𝑖 ∈ [0.05, 5] 𝑦𝑒𝑎𝑟−1 with 𝑖 ∈ {𝑓, 𝑚}. We also assume the objective functionals. The cost of the nonavalent vaccine Abstract and Applied Analysis 7 0.7 0.7 0.6 0.6 0.5 0.5 Individuals Individuals 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0.0 0.0 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Time (years) Time (years) Uf + If Im Uf + If Im Vf Vm Vf Vm (a) (b) Figure 1: Behavior of model (14) in the absence of controls (a) and with constant controls (b). In case (b), the values of controls are 𝑤1 (𝑡) = 0.2, 𝑤2 (𝑡) = 0.1, 𝑢1 (𝑡) = 0.5, 𝑢2 (𝑡) = 0.2, and 𝛼(𝑡) = 0.2 for all 𝑡. The initial conditions are 𝑈𝑓 (0) = 0.05, 𝐼𝑓 (0) = 0.05, 𝑉𝑓 (0) = 0.1, 𝐼𝑚 (0) = 0.1, and 𝑉𝑚 (0) = 0.05. 0.0175 0.7 0.0150 0.6 0.0125 0.5 0.0100 Individuals Controls 0.4 0.0075 0.3 0.0050 0.2 0.0025 0.1 0.0000 0.0 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Time (years) Time (years) w1 u2 Uf + If Im w2  Vf Vm u1 (a) (b) Figure 2: Solution (a) and the optimal controls (b) of the control model (14) with objective functional (16) and weight parameters 𝐴 0 = 1, 𝐴 1 = 3, 𝐴 2 = 10, and 𝐴 3 = 3. Here, the initial conditions are 𝑈𝑓 (0) = 0.05, 𝐼𝑓 (0) = 0.05, 𝑉𝑓 (0) = 0.1, 𝐼𝑚 (0) = 0.1, and 𝑉𝑚 (0) = 0.05. Gardasil-9 is approximately $134.26 per dose [18]. However, individuals (including the cost of cancer induced by HPV) the associated weight parameters in functionals (16) and (17), is assumed less costly than mass vaccination, then, in terms namely, 𝐴 1 and 𝐴 2 , also include costs attributable to vaccine of minimizing cost, it can be logical to stop intervention distribution, storage, and other administration-associated programs. Therefore, in such a case, one can expect low costs. It can be expected that administering the vaccine in vaccination and screening rates as the optimal policies. This school preadolescents is less expensive than vaccination of consideration is confirmed by the results of the computations the sexually active adults; that is, 𝐴 1 < 𝐴 2 [17]. Note also shown in Figure 2 where we plotted the optimal control that the weights 𝐴 0 and 𝐴 ̃0 are related to the cost of having corresponding to the cost functional (16) for infected individuals for this; they are different from weights 𝐴1 𝐴 𝐴 𝐴 1 , 𝐴 2 , 𝐴 3 which reflect the cost of the control interventions. = 3, 2 = 10, 3 = 3. (23) ̃0 and their rations to 𝐴 1 , 𝐴 2 , and 𝐴0 𝐴0 𝐴0 Thus, the values of 𝐴 0 and 𝐴 𝐴 3 are particularly important for the optimal controls. For Please, observe that for this case the magnitude of the example, if 𝐴 0 , 𝐴 ̃0 < 𝐴 1 , 𝐴 2 , 𝐴 3 , that is, if treating infected controls is very small (Figure 2(b)) and, consequently, the 8 Abstract and Applied Analysis optimal controls are insufficient to reduce the spread of the controls for functional (16), this time the controls do not infection (Figure 2(a)). Therefore, in what follows, we assume start with high rates. This allows the number of infected that the costs of having infected individuals are greater individuals in the population to grow during the early phase than the costs of administering the vaccine and screening of the epidemic. In spite of this, the optimal controls manage programs. In particular, we assume that the values of weights to eliminate the infection at the end of the time interval. 𝐴 0 and 𝐴 ̃0 are two or three orders of magnitude bigger than Moreover, the infection levels in Figures 4(a), 4(c), and 4(e) the value of 𝐴 1 , 𝐴 2 , and 𝐴 3 . As a base case, we also set 𝑢𝑚𝑎𝑥 = are clearly higher than those in Figures 3(a), 3(c), and 3(e), 3 and 𝛼𝑚𝑎𝑥 = 3. Before going any further to computations, respectively. Therefore, even though the controls for both we stress that, for problems of infectious disease control, the functionals achieve the elimination of the infection, the initial conditions can play an important role in the search of controls for protection of the vulnerable group do it more effi- optimal controls. However, most studies overlook this issue. ciently. Nevertheless, there is considerably more vaccination In this paper, we consider different initial conditions and in this case than for the case of eradication of the infection. explore how the optimal control varies depending on them. Thus, the controls for functional (16) are presumably more We choose three different sets of initial conditions to reflect expensive than the controls for functional (17). precise infection levels, namely, Figure 5 presents a more direct comparison of relevant model outcomes under objective functionals (16) and (17). (i) Low prevalence levels for both females and males Figure 5(a) shows the total infected population as function (Figures 3(a), 3(b), 4(a), and 4(b)) of time under the optimal controls. Figure 5(b) shows the (ii) High prevalence for females and low prevalence for function 𝐴 1 (𝑤1∗ (𝑡)+𝑤2∗ (𝑡))+𝐴 2 (𝑢∗1 (𝑡)+𝑢∗2 (𝑡))+𝐴 3 𝛼∗ (𝑡) which males (Figures 3(c), 3(d), 4(c), and 4(d)) represents the cost of applying the optimal controls (without (iii) Low prevalence for women and high prevalence for the cost of having infected individuals in the population). males (Figures 3(e), 3(f), 4(e), and 4(f)) Figure 5(c) shows the reproduction number in the presence of vaccination and screening R0 as function of time using the We start the numerical simulations considering the control time series of the optimal controls. The blue dashed lines are model (14) with the objective functional (16) that aims to the outcomes associated with functional (16) whose aim is to protect the vulnerable group. In other words, we seek to protect the vulnerable group. Meanwhile, the red solid lines minimize the number of infections and cervical cancer cases are associated with functional (17) whose aim is to eradicate in females. In this scenario, the optimal controls are charac- the infection. terized by system (22). In Figure 3, we show the numerical From Figure 5 it is clear that the total prevalence is lower, solution obtained via the FBSM. Figures 3(b), 3(d), and 3(f) for almost all times, when the objective is the protection illustrate the profiles of the optimal vaccination and screening of the vulnerable group (Figure 5(a)). Probably, this is a rates, whereas Figures 3(a), 3(c), and 3(e) show the dynamics consequence of the huge difference in the magnitude of the of the state variables. For the three initial conditions, we vaccination rates for both functionals at the beginning of the observe that 𝑤2 (𝑡) ≤ 𝑤1 (𝑡) and 𝑢2 (𝑡) ≤ 𝑢1 (𝑡). This implies that time interval. Nevertheless, the cost of applying the control the vaccination of females is more significant than increasing policies is less, on average, for the objective of eradication male’s immunization coverage. Nevertheless, we remark that, of the infection (Figure 5(b)). Yet, in the second half of the regardless of the independence of the objective functional time interval, there is more vaccination occurring for the (16) from male prevalence, the magnitude of the vaccination objective of eradication, resulting in greater costs for the rates for males and females are nearly identical. Thus, to controls associated with functional (17). This also explains a certain extent, males’ vaccination can be as important as why the value of the reproduction number R0 is less for the females’ vaccination. Furthermore, comparing Figures 3(b), objective of protection at the beginning of the time interval, 3(d), and 3(f), one can see that the optimal control profiles but after a while, its value increases exceeding the value of the for the screening and vaccination rates are qualitatively and R0 associated with the objective of eradication (Figure 5(c)). quantitatively similar. In other words, the optimal controls do Finally, we stress that the cost of having infected individuals not depend on the initial conditions. Please note that for all (which is of paramount importance) is not considered in the scenarios explored, the optimal controls start with very Figure 5(b), therefore it is still uncertain which policy has the high vaccination rates which quickly reduce the number of highest cost-effectiveness under more general hypotheses. infected individuals to very low levels. These controls are then gradually reduced in value but are maintained enough time to 6. Discussion and Conclusion achieve the ultimate eradication of the infection. Our next task is to investigate the profiles of the optimal In this paper, we investigate the optimal control policies vaccination and screening rates with the cost functional (17) against HPV infection. To address this problem, we intro- that aims to minimize the total prevalence at the final time duced a two-sex compartmental epidemic model for the 𝑇. We compute the numerical solution of the two-point transmission dynamics of HPV in a heterosexual population. boundary value problem for the maximum principle via the We incorporated five controls into this model which are able FBSM (Figure 4). It is notable that also for the functional (17), to act simultaneously. These five controls comprise all prac- the optimal vaccination rates for both sexes are remarkably tical intervention policies used against HPV. Using the sex- similar for each of the initial conditions, respectively (see specific infection transfers that we found, we computed the Figures 4(b), 4(d), and 4(f)). Even so, in contrast with the basic reproduction number R0 in the presence of vaccination Abstract and Applied Analysis 9 3.0 0.8 2.5 0.6 2.0 Individuals Controls 0.4 1.5 1.0 0.2 0.5 0.0 0.0 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Time (years) Time (years) Uf + If Im w1 u2 Vf Vm w2  u1 (a) (b) 3.0 0.8 2.5 0.6 2.0 Individuals Controls 0.4 1.5 1.0 0.2 0.5 0.0 0.0 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Time (years) Time (years) Uf + If Im w1 u2 Vf Vm w2  u1 (c) (d) 3.0 0.8 2.5 0.6 2.0 Individuals Controls 0.4 1.5 1.0 0.2 0.5 0.0 0.0 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Time (years) Time (years) Uf + If Im w1 u2 Vf Vm w2  u1 (e) (f) Figure 3: (b, d, f) Optimal vaccination and screening rates as functions of time for the optimal control problem with cost functional (16) for 30 years. Here, 𝐴 0 = 103 , 𝐴 1 = 3, 𝐴 2 = 10, and 𝐴 3 = 3. (a, c, e) Behavior of the state variables for the optimal controls in the same row. (a, b) Initial condition 𝑥1 (0) = (0.05, 0.05, 0.01, 0.05, 0.01)𝑇 . (c, d) Initial condition 𝑥2 (0) = (0.2, 0.1, 0.01, 0.05, 0.01)𝑇 . (e, f) Initial condition 𝑥3 (0) = (0.05, 0.05, 0.01, 0.3, 0.01)𝑇 . and screening. Using the van den Driessche theorem [23], we usually a fragile strategy in a sense that it does not consider established the local stability of the disease-free equilibrium the changing dynamics of the disease’s spread and can be for values of R0 less than one. This stability implies that, nonoptimal in terms of cost-efficacy. if initial infection levels are sufficiently low, maintaining A key issue in the implementation of HPV-vaccination R0 < 1 ensures disease’s elimination and ultimate eradication is whether vaccinating both sexes is more effective than of the virus. However, such a constant control scheme is vaccinating the females only (as it is currently done in 10 Abstract and Applied Analysis 0.8 1.2 0.7 0.6 1.0 0.5 0.8 Individuals Controls 0.4 0.6 0.3 0.4 0.2 0.2 0.1 0.0 0.0 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Time (years) Time (years) w1 u2 Uf + If Im w2  Vf Vm u1 (a) (b) 0.8 0.6 0.7 0.5 0.6 0.4 Individuals 0.5 0.4 Controls 0.3 0.3 0.2 0.2 0.1 0.1 0.0 0.0 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Time (years) Time (years) w1 u2 Uf + If Im w2  Vf Vm u1 (c) (d) 0.8 0.6 0.7 0.5 0.6 0.4 Individuals 0.5 Controls 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0.0 0.0 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Time (years) Time (years) w1 u2 Uf + If Im w2  Vf Vm u1 (e) (f) Figure 4: (b, d, f) Optimal vaccination and screening rates as functions of time for the optimal control problem with cost functional (17) for 30 years. Here 𝐴 ̃0 = 103 , 𝐴 1 = 3, 𝐴 2 = 10, and 𝐴 3 = 3. (a, c, e) Behavior of the state variables for the optimal controls in the same row. (a, b) Initial condition 𝑥1 (0) = (0.05, 0.05, 0.01, 0.05, 0.01)𝑇 . (c, d) Initial condition 𝑥2 (0) = (0.2, 0.1, 0.01, 0.05, 0.01)𝑇 . (e, f) Initial condition 𝑥3 (0) = (0.05, 0.05, 0.01, 0.3, 0.01)𝑇 . most countries). Moreover, although vaccinating individuals theory to identify an effective distribution strategy for HPV- before their potential exposure to HPV infection can be vaccines. We considered two different optimal control prob- a strategy with a positive outcome, vaccination of older lems. Firstly, we studied a rather usual problem in the individuals may be necessary if vaccine-induced protec- literature: minimizing the cumulative number of infected tion is not lifelong. Here, we applied the optimal control females and the cost of intervention policies over a given Abstract and Applied Analysis 11 0.5 70 3.0 60 2.5 Total prevalence 0.4 50 2.0 0.3 40 Cost ℛ0 30 1.5 0.2 20 1.0 0.1 10 0.5 0.0 0 0 5 10 15 20 25 30 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Time (years) Time (years) Time (years) Cost (Protection) ℛ0 (Protection) Prevalence (Protection) ℛ0 (Eradication) Cost (Eradication) Prevalence (Eradication) (a) (b) (c) Figure 5: Comparison of the total prevalence (a), the cost of applying the optimal control interventions (b), and the value of R0 as function of the optimal controls (c). The blue dashed lines are associated with functional (16) and the red solid lines are associated with functional (17). Here, 𝐴̃0 = 𝐴 0 = 103 , 𝐴 1 = 3, 𝐴 2 = 10, 𝐴 3 = 3, and the initial condition is 𝑥3 (0) = (0.05, 0.05, 0.01, 0.3, 0.01)𝑇 . time interval. Such a problem corresponds to a policy that Nevertheless, due to a possible nonlinear growth of the costs aimed at protection of the most vulnerable group of the with respect to the fraction of the vaccinated individuals, even population, that is, females at risk of developing invasive if the vaccination coverage of females is already high, it may cancer. Secondly, we considered the problem of minimizing still be cost-effective to include males into the vaccination the total infection level in both sexes at the end of a given program. Besides, there can be a group who is negligent or time interval at the lowest possible costs. This problem has a reluctant towards the vaccination; thus, a 100% vaccination more ambitious objective, namely, the ultimate eradication of coverage can be difficult to achieve even in one sex. Finally, HPV. Such a goal is possible if the basic reproduction number we like to note that the impacts of vaccine coverage on the is comparatively low. The possibility of total eradication profiles of the optimal vaccination rates and the dependency of cervical cancer worldwide is currently discussed in the of principal properties of the controls on the forms of scientific community [31]. objective functionals remain open questions and deserve Results of our numerical simulations show that vaccina- further study. tion should be administered both prior and after the sexual debut and in both, females and males. For both functionals, Appendix the simulations support the intuitively expected outcome that females’ vaccination should be the priority; however, at Here, we compute the adjoint system and the characterization the same time, we show that males’ vaccination is almost of the optimal controls considering the objective functional as beneficial and required. If the objective is to protect (17): the vulnerable group, the numerical experiments imply that initially, higher vaccination rates should be applied which can ̃0 𝑃 (𝑇) + 1 𝐽̃ (𝑐) = 𝐴 gradually be reduced. This policy leads to a fast reduction 2 in the number of infectious individuals. For the objective of (A.1) 𝑇 eradication of the infection, the optimal controls begin at low ⋅ ∫ (𝐴 1 (𝑤12 + 𝑤22 ) + 𝐴 2 (𝑢21 + 𝑢22 ) + 𝐴 3 𝛼2 ) 𝑑𝑡. 0 rates allowing the prevalence to increase. Yet, the optimal vaccination and screening rates increase over time achieving In this case, the optimal control problem of the HPV the elimination of the infection at the end of the time interval. epidemic model is We have to stress that our results differ from the majority of previous studies where the suggested preferred strategy min 𝐽̃ (𝑐) (A.2) 𝑐∈𝐷(𝑇) was the vaccination of the preadolescent females prior to their sexual initiation. We believe that such a preference for subject to the dynamics of the control model (14) with given the young females’ vaccination that prevails in the literature initial conditions. arose because the majority of authors did not consider the The Pontryagin maximum principle converts the optimal problem of HPV control as an optimal control problem. control problem (A.2) into a problem of minimizing point- This also implies that the “young females only” vaccination wise the Hamiltonian is neither the optimal nor the most cost-effective policy. Furthermore, the efficacy of an HPV reduction program may ̃ = 1 [𝐴 1 (𝑤2 + 𝑤2 + 𝐴 2 (𝑢2 + 𝑢2 ) + 𝐴 3 𝛼2 ] 𝐻 1 2 1 2 also depend on the actual vaccine coverage. Our results were 2 (A.3) obtained for initially low vaccine coverage. This corresponds to a setting where HPV vaccination is just being introduced. + 𝜓1 𝑈̇ 𝑓 + 𝜓2 𝐼𝑓̇ + 𝜓3 𝑉𝑓̇ + 𝜓4 𝐼𝑚̇ + 𝜓5 𝑉𝑚 ̇ 12 Abstract and Applied Analysis with respect to the controls. Here, 𝜓𝑖 for 𝑖 = 1, . . . , 5 are Acknowledgments adjoint variables that satisfy the system of differential equa- tions Fernando Saldaña acknowledges the support by CONACyT, Mexico, under grant number 331294. Andrei Korobeinikov 𝜓1̇ = [(1 − 𝑝) 𝛽𝑚 𝐼𝑚 + 𝛾𝑓 + 𝛼∗ + 𝜇𝑓 ] 𝜓1 is supported by Spain’s Ministerio de Ciencia, Innovación y Universidades grant MTM2015-71509-C2-1-R. + (𝑝𝛽𝑚 𝐼𝑚 − 𝛼∗ ) 𝜓2 + 𝑢∗1 𝜓3 − 𝛽𝑓 [(1 − 𝐼𝑚 − 𝑉𝑚 ) + 𝜖𝑉𝑚 ] 𝜓4 + 𝛽𝑓 𝜖𝑉𝑚 𝜓5 , References [1] WHO, “Cervical cancer,” 2018, https://www.who.int/cancer/ 𝜓2̇ = (1 − 𝑝) 𝛽𝑚 𝐼𝑚 𝜓1 + (𝑝𝛽𝑚 𝐼𝑚 + 𝛾𝑓 + 𝜇𝑓 ) 𝜓2 prevention/diagnosis-screening/cervical-cancer/en/. [2] FDA, “Fda approves expanded use of gardasil-9 to include indi- + 𝑢∗1 𝜓3 − 𝛽̃𝑓 [(1 − 𝐼𝑚 − 𝑉𝑚 ) + 𝜖𝑉𝑚 ] 𝜓4 viduals 27 through 45 years old,” 2018, https://bit.ly/2zVRxa9. + 𝛽̃𝑓 𝜖𝑉𝑚 𝜓5 , [3] M. Brisson, J.-F. Laprise, H. W. Chesson et al., “Health and economic impact of switching from a 4-valent to a 9-valent hpv vaccination program in the united states,” Journal of the 𝜓3̇ = [(1 − 𝑝) 𝜓1 + 𝑝𝜓2 ] (1 − 𝜖) 𝛽𝑚 𝐼𝑚 (A.4) National Cancer Institute, vol. 108, no. 1, 2016. + (𝜖𝛽𝑚 𝐼𝑚 + 𝑢∗1 + 𝜇𝑓 + 𝜃) 𝜓3 , [4] K. Seto, F. Marra, A. Raymakers, and C. A. Marra, “The cost effectiveness of human papillomavirus vaccines: a systematic review,” Drugs, vol. 72, no. 5, pp. 715–743, 2012. 𝜓4̇ = −𝛽𝑚 ((1 − 𝑈𝑓 − 𝐼𝑓 − 𝑉𝑓 ) + 𝜖𝑉𝑓 ) [5] G. P. Garnett, J. J. Kim, K. French, and S. J. Goldie, “Chapter 21: ⋅ [(1 − 𝑝) 𝜓1 + 𝑝𝜓2 ] + 𝜖𝛽𝑚 𝑉𝑓 𝜓3 Modelling the impact of HPV vaccines on cervical cancer and screening programmes,” Vaccine, vol. 24, no. 3, pp. S178–S186, + (𝛽𝑓 𝑈𝑓 + 𝛽̃𝑓 𝐼𝑓 + 𝛾𝑚 + 𝜇𝑚 ) 𝜓4 + 𝑢∗2 𝜓5 , 2006. [6] A. T. Newall, P. Beutels, J. G. Wood, W. J. Edmunds, and C. 𝜓5̇ = (𝛽𝑓 𝑈𝑓 + 𝛽̃𝑓 𝐼𝑓 ) ((1 − 𝜖) 𝜓4 + 𝜖𝜓5 ) R. MacIntyre, “Cost-effectiveness analyses of human papillo- mavirus vaccination,” The Lancet Infectious Diseases, vol. 7, no. 4, pp. 289–296, 2007. + (𝑢∗2 + 𝜃 + 𝜇𝑚 ) 𝜓5 , [7] E. H. Elbasha, E. J. Dasbach, and R. P. Insinga, “Model for ̃0 for 𝑖 = 1, 2, 4, and assessing human papillomavirus vaccination strategies,” Emerg- with transversality conditions 𝜓𝑖 (𝑇) = 𝐴 ing Infectious Diseases, vol. 13, no. 1, pp. 28–41, 2007. 𝜓𝑘 (𝑇) = 0 for 𝑘 = 3, 5. For this control problem, applying [8] J. J. Kim and S. J. Goldie, “Cost effectiveness analysis of includ- the maximum principle from [32], we obtain the following ing boys in a human papillomavirus vaccination programme characterization for the optimal controls: in the United States,” British Medical Journal, vol. 339, article 𝜇𝑓 b3884, 2009. 𝑤1∗ (𝑡) = min {1, max {0, − 𝜓3 (𝑡)}} , [9] E. H. Elbasha and E. J. Dasbach, “Impact of vaccinating boys 𝐴1 and men against HPV in the United States,” Vaccine, vol. 28, no. 𝜇𝑚 𝑤2∗ (𝑡) = min {1, max {0, − 𝜓 (𝑡)}} , 42, pp. 6858–6867, 2010. 𝐴1 5 [10] M. A. Smith, J.-B. Lew, R. J. Walker, J. M. L. Brotherton, 𝑢∗1 (𝑡) C. Nickson, and K. Canfell, “The predicted impact of HPV 1 vaccination on male infections and male HPV-related cancers = min {𝑢𝑚𝑎𝑥 , max {0, − (1 − 𝑈𝑓 (𝑡) − 𝐼𝑓 (𝑡) − 𝑉𝑓 (𝑡)) 𝜓3 (𝑡)}} , in Australia,” Vaccine, vol. 29, no. 48, pp. 9112–9122, 2011. 𝐴2 (A.5) 𝑢∗2 (𝑡) [11] J. A. Bogaards, M. Kretzschmar, M. Xiridou, C. J. L. M. Meijer, J. Berkhof, and J. Wallinga, “Sex-specific immunization for 1 sexually transmitted infections such as human papillomavirus: = min {𝑢𝑚𝑎𝑥 , max {0, − (1 − 𝐼𝑚 (𝑡) − 𝑉𝑚 (𝑡)) 𝜓5 (𝑡)}} , 𝐴2 insights from mathematical models,” PLoS Medicine, vol. 8, no. 𝛼∗ (𝑡) 12, article e1001147, 2011. 1 [12] T. Malik, J. Reimer, A. Gumel, E. Elbasha, and S. Mahmud, = min {𝛼𝑚𝑎𝑥 , max {0, (𝜓 (𝑡) − 𝜓2 (𝑡)) 𝑈𝑓 (𝑡)}} . “The impact of an imperfect vaccine and Pap cytology screening 𝐴3 1 on the transmission of human papillomavirus and occurrence of associated cervical dysplasia and cancer,” Mathematical Data Availability Biosciences and Engineering, vol. 10, no. 4, pp. 1173–1205, 2013. [13] A. A. Alsaleh and A. B. Gumel, “Analysis of risk-structured The software and data of computation used to support the vaccination model for the dynamics of oncogenic and warts- findings of this study are available from the corresponding causing HPV types,” Bulletin of Mathematical Biology, vol. 76, author upon request. no. 7, pp. 1670–1726, 2014. [14] A. M. Guerrero, A. J. Genuino, M. Santillan et al., “A cost-utility Conflicts of Interest analysis of cervical cancer screening and human papillomavirus vaccination in the Philippines,” BMC Public Health, vol. 15, no. The authors declare that they have no conflicts of interest. 1, article 730, 2015. Abstract and Applied Analysis 13 [15] F. Saldaña and I. Barradas, “The role of behavioral changes and prompt treatment in the control of STIs,” Infectious Disease Modelling, vol. 4, pp. 1–10, 2019. [16] O. Sharomi and T. Malik, “Optimal control in epidemiology,” Annals of Operations Research, vol. 251, no. 1-2, pp. 55–71, 2017. [17] V. L. Brown and K. A. Jane White, “The role of optimal control in assessing the most cost-effective implementation of a vaccination programme: HPV as a case study,” Mathematical Biosciences, vol. 231, no. 2, pp. 126–134, 2011. [18] T. Malik, M. Imran, and R. Jayaraman, “Optimal control with multiple human papillomavirus vaccines,” Journal of Theoretical Biology, vol. 393, pp. 179–193, 2016. [19] W. K. Huh, E. A. Joura, A. R. Giuliano et al., “Final efficacy, immunogenicity, and safety analyses of a nine-valent human papillomavirus vaccine in women aged 16–26 years: a ran- domised, double-blind trial,” The Lancet, vol. 390, no. 10108, pp. 2143–2159, 2017. [20] N. Muñoz, F. Méndez, H. Posso et al., “Incidence, duration, and determinants of cervical human papillomavirus infection in a cohort of Colombian women with normal cytological results,” The Journal of Infectious Diseases, vol. 190, no. 12, pp. 2077–2087, 2004. [21] G. M. Anic and A. R. Giuliano, “Genital HPV infection and related lesions in men,” Preventive Medicine, vol. 53, no. 1, pp. S36–S41, 2011. [22] M. Martcheva, An Introduction to Mathematical Epidemiology, vol. 61 of Texts in Applied Mathematics, Springer, New York, NY, USA, 2015. [23] P. van den Driessche and J. Watmough, “Reproduction numbers and sub-threshold endemic equilibria for compartmental mod- els of disease transmission,” Mathematical Biosciences, vol. 180, pp. 29–48, 2002. [24] E. Grigorieva, E. Khailov, and A. Korobeinikov, “Optimal control for an SEIR epidemic model with nonlinear incidence rate,” Studies in Applied Mathematics, vol. 141, no. 3, pp. 353–398, 2018. [25] B. Buonomo, D. Lacitignola, and C. Vargas-De-Leon, “Quali- tative analysis and optimal control of an epidemic model with vaccination and treatment,” Mathematics and Computers in Simulation, vol. 100, pp. 88–102, 2014. [26] K. O. Okosun, O. Rachid, and N. Marcus, “Optimal control strategies and cost-effectiveness analysis of a malaria model,” BioSystems, vol. 111, no. 2, pp. 83–101, 2013. [27] A. Camacho and S. Jerez, “Bone metastasis treatment modeling via optimal control,” Journal of Mathematical Biology, pp. 1–30, 2018. [28] W. H. Fleming and R. W. Rishel, Deterministic and Stochastic Optimal Control, vol. 1, Springer, New York, NY, USA, 1975. [29] D. L. Lukes, Differential Equations: Classical to Controlled, Elsevier, 1982. [30] S. M. Lenhart and J. T. Workman, Optimal Control Applied to Biological Models, CRC Press, 2007. [31] I. Baussano and F. Bray, “Modelling cervical cancer elimination,” The Lancet Public Health, vol. 4, no. 1, pp. e2–e3, 2019. [32] L. S. Pontryagin, Mathematical Theory of Optimal Processes, Routledge, 2018. Advances in Advances in Journal of The Scientific Journal of Operations Research Hindawi Decision Sciences Hindawi Applied Mathematics Hindawi World Journal Hindawi Publishing Corporation Probability and Statistics Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 http://www.hindawi.com www.hindawi.com Volume 2018 2013 www.hindawi.com Volume 2018 International Journal of Mathematics and Mathematical Sciences Journal of Hindawi Optimization Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 Submit your manuscripts at www.hindawi.com International Journal of Engineering International Journal of Mathematics Hindawi Analysis Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 Journal of Advances in Mathematical Problems International Journal of Discrete Dynamics in Complex Analysis Hindawi Numerical Analysis Hindawi in Engineering Hindawi Differential Equations Hindawi Nature and Society Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 International Journal of Journal of Journal of Abstract and Advances in Stochastic Analysis Hindawi Mathematics Hindawi Function Spaces Hindawi Applied Analysis Hindawi Mathematical Physics Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018