key: cord-0255492-tbzvzvf4 authors: Omame, A.; Okuonghae, D.; Nwafor, U. E.; Odionyenma, B. U. title: A Co-infection model for HPV and Syphilis with Optimal Control and Cost-Effectiveness Analysis date: 2020-09-11 journal: nan DOI: 10.1101/2020.09.09.20191635 sha: 760cdb97c1feaef89513cfde733b61d9cfa5dc3c doc_id: 255492 cord_uid: tbzvzvf4 In this work, we develop and present a co-infection model for human papillomavirus (HPV) and syphilis with cost-effectiveness optimal control analysis. The full co-infection model is shown to undergo the phenomenon of backward bifurcation when a certain condition is satisfied. The global asymptotic stability of the disease-free equilibrium of the full model is shown textbf{not to exist}, when the associated reproduction number is less than unity. The existence of endemic equilibrium of the syphilis-only sub-model is shown to exist and the global asymptotic stability of the disease-free and endemic equilibria of both the syphilis-only sub-model and HPV-only sub-model were established. The global asymptotic stability of disease-free equilibrium of the HPV-only sub-model is also proven. Numerical simulations of the optimal control model showed that the optimal control strategy which implements syphilis treatment controls for singly infected individuals is the most cost-effective of all the control strategies in reducing the burden of HPV and syphilis co-infections. (TB), in the presence of HPV vaccination and condom use as well as TB treatment. Their results showed that TB-only treatment strategy can significantly bring down the burden of HPV and the co-infection of the two diseases in a population. Nwankwo and Okuonghae [24] investigated the effect of syphilis treatment in a population where HIV Syphilis are co-circulating, and showed that targeted syphilis treatment could significantly curb the burden of the co-infection of the two diseases. Also, an optimal control model for two-strain tuberculosis and HIV/AIDS co-infection with cost-effectiveness analysis was studied by Agusto and Adekunle [5] . They showed that the control strategy that combines prevention of treatment failure in drug-sensitive TB infectious individuals and the treatment of individuals with drug-resistant TB was the most cost-effective in reducing the burden of the co-infection of HIV/AIDS and two strains of tuberculosis. To the best of the authors' knowledge, no robust optimal control mathematical model has been developed to capture the combined effect HPV vaccination and syphilis treatment on the control of HPV and syphilis co-infection, despite the availability of epidemiological evidences supporting the coinfections of both diseases. This study assesses the combined effect of these aforementioned strategies, using optimal control analysis, in order to determine ways of bringing down the bnurden of co-infection of the two diseases. The model is based on the main assumptions stated below: i To avoid model complexity, the primary and secondary stages of syphilis infection are joined together and are referred to as "early stage of syphilis infection". In a similar manner, the latent and tertiary stages of infection are referred to combined together and known as "late stage of syphilis infection". ii. individuals infected with syphilis infection is susceptible to infection with HPV and vice versa. [34] . iii. individuals dually infected with HPV and syphilis can transmit either HPV or syphilis but not mixed infections. The paper is organized as follows. The model is formulated in Section 2, together with the presentation of its basic properties. The sub-models are analysed in Section 3. Qualitative and quatitative analyses of the full co-infection model with constant and time dependent controls are presented in Section 5 while Section 6 gives the concluding remarks. depicted in Figure 1 , the associated state variables and parameters are well described in Tables 1 and 2 . = σ 4 H I + γ 2 λ h L − (τ 4 + r 5 + µ + δ s2 + δ h4 )H L dP E dt = ε 2 λ s P − (r 6 + µ + σ 5 )P E + ρ 2 r 3 H E dP I dt = σ 5 P E − (τ 5 + r 7 + µ + σ 6 )P I + ρ 3 r 4 H I dP L dt = σ 6 P I − (τ 6 + r 8 + µ + δ s3 )P L + ρ 4 r 5 H L dC E dt = λ s C + χ 1 r 6 P E − (µ + δ C + σ 7 + α 2 )C E dC I dt = χ 2 r 7 P I − (µ + δ C + σ 8 + α 3 + τ 7 )C I + σ 7 C E dC L dt = χ 3 r 8 P L − (µ + δ C + δ s4 + α 4 + τ 8 )C L + σ 8 C I where, λ s = β s [I + C I + η l (L + C L ) + θ h {H I + ω P P I + η l (H L + ω P P L )}] N λ h = β h [H + H E + ω P (P + P E ) + θ s {H I + ω P P I + η l (H L + ω P P L )}] N The parameter θ h (θ h ≥ 1) is a modification parameter accounting for the increased infectiousness of dually infected individuals due to active HPV, η l (η l ≥ 1) is a modification parameter accounting for increased infectiousness of infected individuals in the latent stage of syphilis, in comparison to those in the early stage of syphilis infection. ω p (ω p ≤ 1) is a modification parameter accounting for the reduced infectiousness of infected individuals due to persistent HPV infection while θ s (θ s ≥ 1) is a modification parameter accounting for the increased infectiousness of dually infected individuals due to syphilis infection. In (2) , β s is the effective contact rate for the transmission of syphilis infection, while, β h is the effective contact rate for the transmission of HPV infection. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint Individuals with persistent HPV and exposed syphilis infection P I Individuals with persistent HPV and syphilis infection in late stage C Individuals infected with anal cancer C E Individuals infected with anal cancer and exposed to syphilis infection C I Individuals infected with anal cancer and syphilis infection in early stage C L Individuals infected with anal cancer and syphilis infection in late stage. R C Individuals who recovered from anal cancer R h Individuals who recovered from HPV 6 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint 2.1 Basic properties of the co-infection model (1) The basic dynamical properties of the model (1) will now be explored. Particularly, we establish the following positivity and invariance results. For the model (1) to be epidemiologically meaningful, it is important to prove that all its state variables are non-negative for all time. Following the same approach in Omame et al [29] , we establish the following results. Then the solutions (S, E, I, L, T, V h , H, P, C, R C , R h , H E , H I , H L , P E , P I , P L , C E , C I , C L ) of the model (1) are non-negative for all time t > 0. Therefore, D is positively invariant. Hence, no solution path can leave through any boundary of D and it is sufficient to consider the dynamics of the model (1) in D. Inside this region, the model is considered to be mathematically and epidemiologically well posed. It is instructive to analyze the sub-models of the co-infection model (1) , before analyzing the full model. The syphilis-only sub-model is (obtained by setting V h = H = P = C = R c = R h = H e = H i = H l = P e = P i = P l = C e = C i = C l = 0 in the model (1)) given by: where now, with N = S + E + I + L + T 7 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint 3.1.1 Basic reproduction number of the syphilis only sub-model The syphilis-only sub-model (3) has a DFE, obtained by setting the right-hand sides of the equations in the model (3) to zero, given by The basic reproduction number, using the next generation operator method [42] , is given by The local stability of the syphilis-only sub-model is analysed by the Jacobian matrix of the system (3) at ξ 0S , given by: where, The eigenvalues are −µ, −µ and the solution of the characteristic polynomial: where Applying the Routh-Hurwitz criterion, the cubic equation (7) will have roots with negative real parts if and only if α 1 > 0, α 3 > 0 and α 1 α 2 > α 3 . Clearly, α 1 > 0 and α 3 > 0 (if R 0S < 0). As a result, the disease-free equilibrium, ξ 0S is locally asymptotically stable if R 0S < 1. 3.1.3 Global asymptotic stability(GAS) of the disease-free equilibrium(DFE) ξ 0S of the syphilis-only sub-model We use the method presented in Castillo-Chavez et al [11] to investigate the global asymptotic stability of the disease free equilibrium of the syphilis-only sub-model. In this section, we list two conditions that if met, also guarantee the global asymptotic stability of the disease-free state. First, system (3) must be written in the form: 8 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . where W ∈ R m denotes (its components) the number of uninfected individuals and Q ∈ R n denotes (its components) the number of infected individuals including latent, infectious, etc. U 0 = (W * , 0) denotes the disease-free equilibrium of this system. The conditions (H1) and (H2) below must be met to guarantee local asymptotic stability: (H1): For dX dt = F (W, 0), W * is globally asymptotically stable (GAS), (H2): G(W, Q) = AI −Ĝ(W, Q)W, G(W, Q) ≥ 0 for (W, Q) ∈ Ω, where A = D 1 G(W * , 0) is an M-matrix (the off-diagonal elements of A are nonnegative) and Ω is the region where the model makes biological sense. If System (3) satisfies the above two conditions then the following theorem holds: The fixed point U 0S = (W * , 0) is a globally asymptotic stable (GAS) equilibrium of (3) provided that R 0 < 1 (LAS) and that assumptions (H1) and (H2) are satisfied where W denotes the number of non-infectious individuals and Q denotes the number of infected individuals It is clear from the above, that,Ĝ(W, Q) 0. Hence the DFE, U 0S , may not be globally asymptotically stable, suggesting the possibility of a backward bifurcation. This supports the backward bifurcation analysis in Section (3.1.5) where, substituting the right-hand sides of (10) into the force of infection (4) at steady states, we obtain the polynomial equation where, It is worthy of note, from the equation (11) , that a 1 > 0, while a 3 is greater than zero (less than zero) if R 0S < 1(> 1). Hence, the following result can be established: 9 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. Table 2 Theorem 3.2 The model (3) has i a unique endemic equilibrium if a 3 < 0 ⇐⇒ R 0S > 1; ii a unique endemic equilibrium if a 2 < 0 and a 1 = 0 or a 2 2 − 4a 3 a 1 = 0; iii two endemic equilibria if a 1 > 0, a 2 < 0 and a 2 2 − 4a 3 a 1 > 0 and R 0S < 1; iv no endemic equilibrium otherwise. It is also interesting to note that, setting the re-infection term ξ s = 0, reduces the quadratic (11) to a 2 λ * * s + a 3 = 0, resulting in no sign changes in the polynomial equation (11) , as a 2 > 0 and a 3 > 0 (for R 0S < 1). Hence no existence of an endemic equilibrium for R 0S < 1, ruling out the existence of backward bifurcation in the syphilis only sub-model (3) in the absence of re-infection of treated individuals. This is consitent with the proof in Theorem 3.1 Using the same approach as in Castillo-Chavez and Song [10] , we shall investigate the possibility of a backward bifurcation for the syphilis-only sub-model. The linearized system evaluate at the DFE is given as  where, The components of the right eigenvector are given as: The components of the left eigenvector are equally given by: CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . It can be observed from (14) , that setting the re-infection parameter, ξ s = 0, results in a < 0. Thus, re-infection induced the phenomenon of backward bifurcation in the Syphilis-only sub-model (3). This is consistent with the results obtained in the analyses in sections 3.1 and 3.1.4. An associated backward bifurcation diagram is depicted in Figure 2 ). 3.1.6 GAS of EEP of the Syphilis-only sub-model (3): special case(ξ s = δ s = 0) Now that the cause of the backward bifurcation in the syphilis-only sub-model (3) is removed, that is, ξ s = 0, we seek to prove the global asymptotic stability of the unique endemic equilibrium of the submodel, for a special case when disease induced death rates are neglibible, that is, δ s = 0. The Syphilis-only sub-model (3) has an endemic equilibrium given by It should be noted that setting Proof. Consider the syphilis-only sub-model (3) with (15) and ξ s = 0 andR 0 > 1, so that the associated unique endemic equilibrium exists. Also, consider the Lyapunov functional similar to the Goh-Volterra type considered by Ghosh et al [15] : The time derivative, using Leibniz rule for integration as illustrated in Adams [2] , is given bẏ Substituting the expressions for the derivatives,Ṡ,Ė,İandL from (3), into the Lyapunov derivative,L, and carrying out certain algebraic manipulations, we have thaṫ 11 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint (16) can be further simplified intȯ Thus,L ≤ 0 forR 0 > 1. Hence, L is a Lyapunov function in D and it follows from the La Salle's Invariance principle [19] , that every solution to the equations of the syphilis-only sub-model (3) with (15) and initial conditions in D D 0 approaches the associated unique endemic equilibrium ξ e , of the syphilis-only sub-model as t → ∞ forR 0S > 1. The epidemiological implication of Theorem 3.3 is that, if a previous infection with syphilis confers lifetime protection against re-infection, then syphilis infection will persist in the population, if the threshold quantity,R 0 > 1. The HPV-only sub-model is (obtained by setting (1)) given by: where now, The HPV-only sub-model (18) has a DFE, obtained by setting the right-hand sides of the equations in the model (18) to zero, given by The basic reproduction number, using the next generation operator method [42] is given by The local and global asymptotic stability analyses of the HPV-only sub-model (18) were well explored by Omame et al. [30] . 12 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . 4 Local asymptotic stability of disease-free equilibrium co-infection model (1) The epidemiological implication of Lemma 4.1 is that when R 0 < 1, a small influx of HPV or syphilisinfected individuals into the population will not generate large HPV or syphilis outbreaks, and the diseases will die out. 13 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . 4.1 Global asymptotic stability of the disease free equilirium of the full co-infection model (1) Using the same approach as in Section 3.1.3, we can establish the following result. Consider the conditions below: (H1): For dW dt = F (W, 0), W * is globally asymptotically stable (GAS), and Ω is the region where the model makes biological sense. If System 1 satisfies the above two conditions then the following theorem holds: 0) is a globally asymptotic stable (GAS) equilibrium of 1 provided that R 0 < 1 (LAS) and that assumptions (H1) and (H2) are satisfied The Proof follows as in section 3.1, wherê It is evident from the above, that,Ĝ(W, Q) 0, which means that the condition (H2) is not satisfied. Hence the DFE, U 0 = (W * , 0) may not be globally asymptotically stable, suggesting the possibility of a backward bifurcation. This supports the backward bifurcation analysis for the full-co-infection model (1) in the next section. We shall investigate the type of bifurcation the model (1) may undergo, using the Centre Manifold Theory as discussed in [10] . The following result can be obtained using the approach in [10] . then model (1) exhibits backward bifurcation at R 0 = 1. If a < 0, then the system (1) exhibits a forward bifurcation at R 0 = 1. Suppose represents any arbitrary endemic equilibrium of the model (that is, an endemic equilibrium in which at least one of the infected components is non-zero). To apply the Centre Manifold Theory, it is necessary 14 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint to carry out the following change of variables. Let Further, using the vector notation (1) can be re-written in the form as follows: x 10 = α 1 x 9 − (µ + λ s )x 10 where, Consider the case when R 0S = 1. Suppose, further, that β s is chosen as a bifurcation parameter. Solving for β s = β * s from R 0S = 1 gives CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint Evaluating the Jacobian of the system (20) at the DFE, J(ξ 0 ), and using the approach in [10] , we have that J(ξ 0 ) has a right eigenvector (associated with the simple zero eigenvalue of J(ξ 0 )) given by w = [ω 1 , ω 2 , ω 3 , ω 4 , ω 5 , ω 6 , ω 7 , ω 8 , ω 9 , ω 10 , ω 11 , ω 12 , ω 13 , ω 14 , ω 15 , ω 16 , ω 17 , ω 18 , ω 19 , ω 20 ] T where, where, K1 = σ1 + µ, K2 = σ2 + τ1 + µ, K3 = τ2 + µ + δs1, K4 = µ + δh1 + r1, K5 = µ + r2, K6 = µ + α1 + δc, K7 = σ3 + µ + r3 + δh2, K8 = σ4 + µ + τ3 + r4 + δh3, K9 = τ4 + r5 + µ + δs2 + δh4, K10 = r6 + µ + σ5, K11 = τ5 + r7 + µ + σ6, K12 = τ6 + r8 + µ + δs3, K13 = µ + δC + σ7 + α2, K14 = µ + δC + σ8 + α3 + τ7, 16 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint It follows from Theorem 4.1 in [10] , by computing the non-zero partial derivatives of F (x) (evaluated at the disease free equilibrium, (ξ0)) that the associated bifurcation coefficients defined by a and b, given by and Since the bifurcation coefficient b is positive, it follows from Theorem 4.1 in [10] that the model (1), or the transformed model (20) , will undergo a backward bifurcation if the backward bifurcation coefficient, a, given by (22) is positive. We apply Pontryagin's Maximum Principle, in this section, to determine the necessary conditions for the optimal control of the HPV-Syphilis co-infection model. We assume that the proportion of vaccinated individuals, φ for HPV and the Syphilis treatment rates u1, u2, u3, u4, u5, u6, u7, u8 are now time dependent and will therefore act as the control variables. Hence we have,Ṡ where, is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint For this, we consider the objective functional J u1, u2, u3, u4, u5, u6, u7, u8, φ = T 0 The parameters, W1 is the weight on the benefit and cost of implementing the optimal control measure, where W1 balances the cost factors as a result of the size and the relevance of the terms making up the objective functional. T is the final time. We seek to find an optimal control, φ * , u * 1 , u * 2 , u * 3 , u * 4 , u * 5 , u * 6 , u * 7 , and u * 8 such that where is the control set. The Pontryagin's Maximum Principle [33] gives the necessary conditions which an optimal control pair must satisfy. This principle transforms (23), (25) and (26) into a problem of minimizing a Hamiltonian, Ham, pointwisely with regards to the control functions, φ, u1, u2, u3, u4, u5, u6, u7 and u8 : Theorem 5.1 For an optimal control set φ, u1, u2, u3, u4, u5, u6, u7, u8 that minimizes J over U , there are adjoint variables, λ1, λ2, ..., λ20 satisfying is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint Furthermore, Proof of Theorem 5.1 Suppose U * = (φ * , u * 1 , u * 2 , . . . , u * 8 ) is an optimal control and S * , E * , I * , L * , T * , V * h , H * , P * , C * , R * C , R * h , H * I , H * L , P * E , P * I , P * L , C * E , C * I , C * L are the corresponding state solutions. Applying the Pontryagin's Maximum Principle [33] , there exist adjoint variables satisfying: with transversality conditions; We can determine the behaviour of the control by differentiating the Hamiltonian, Ham with respect to the controls (φ, u1, u2, u3, u4, u5, u6, u7, u8) at t. On the interior of the control set, where 0 < uj < 1 for all (j = 1, ..., 8) and 0 < φ < 1, we obtain Therefore, we have that [20] is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint We now simulate the optimal control model (23) numerically using the parameter estimates in Table 2 , so that the reproduction number, R0 = 1.80842 (unless otherwise stated), to assess the potential impact of various targeted control strategies on the transmission dynamics of HPV and SYphilis in the population. Demographic parameters relevant to the city of Rio de Janeiro in Brazil were chosen. Specifically, since the total population of sexually active susceptible individuals (15-64 years) in the State of Rio de Janeiro in Brazil are estimated to be 12,850,804, at disease-free equilibrium, Λ µ = 8, 946, 246 [9] . In Brazil, the life expectancy is estimated at 74 years [9] . Hence we have that µ = 0.0135, so that Λ = 120895 per year. Numerical simulations of the optimal control problem (23), adjoint equations (30) and characterizations of the control (33) are implemented by the Runge Kutta method using the forward backward sweep (carried out in MATLAB). The balancing factor w1 = 10 3 (The Centres for Disease Control cost per dose of the pediatric Gardasil 4, the Gardasil 9, and the Cervarix, respectively, is $121.03,$134.26 and $107.97 [12, 22] . Since the associated costs considered are combinations of price per dose, storage/ administration-related costs and transportation, etc, we will assume w1 = 10 3 ). Following the estimates for syphilis treatment in [18] , the balancing factors for syphilis-infected individuals are: w2 = w3 = 200, w4 = w5 = w6 = w7 = w8 = w9 = 220. Here, we assume that the cost of treatment for co-infected individuals is more than the cost of treatment for singly-infected individuals. Following the reported prevalence of Syphilis and HPV as well as the co-infection prevalence in [34, 39] and the demographic data obtained from [9] , we set the initial conditions at: Applying this strategy, we note in Figures 3 (a) and 3(b) that the total number of individuals singly infected with HPV and persistent HPV, respectively, is less than the total number when no control is applied. It can equally be observed that optimal HPV vaccination strategy has a high positive population level impact on the populations of infected individuals dually infected with HPV and Syphilis in early and late stages of infection, respectively, (Figures 4 (a) and 4(b) ), infected individuals dually infected with persistent HPV and Syphilis in early and late stages of infection,respectively ( Figure 5 (a) and 5(b)) and infected individuals dually infected with Cancer and Syphilis in early and late stages of infection, respectively (Figure 6 (a) and 6(b)). The results of the simulations agree with the epidemiological reports in [26, 25] , that prior HPV infections are risk factors for syphilis infection. Therefore, controlling HPV infections through vaccination can bring down the burden of co-infection of the two diseases. The simulations of the total number of dually infected individuals in the presence of Syphilis treatment controls are depicted in Figures 6a -6d) . Applying this control, we observe that the total number of individuals dually infected with HPV and Syphilis in early and late stages of infection respectively, is less than the total population when no control is applied as expected (Figure 6 (a) and 6(b)). Similarly, it is noticed from Figures (6 (c) and 6 (d), that the total number of individuals dually infected with persistent HPV and Syphilis in early and late stages of infection, respectively, is lesser when this control is applied. In addition high population level impact is noticed on the total number of individuals dually infected with Cancer and Syphilis in early and late stages of infection, respectively, as shown in (Figure 13 (a) and 13 (b) . This supports the epidemiological report in the introduction section that syphilis is a risk factor for HPV infection [34] . Hence, if we focus on syphilis treatment controls, it can significantly bring down the burden of the co-infection of HPV and syphilis in a population. The simulations equally agree with the findings in Tseng et al. [35] , that prior syphilis infection was associated with persistent HPV and increased susceptibility to anal cancer. As a result, treating syphilis infection in individuals dually 20 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint and 4d) , when optimal vaccination control only strategy (φ = 0) and syphilis treatment controls for singly-infected (u 1 = u 2 = 0) are administered. Here, β s = 7.0, β h = 2.0. All other parameters as in Table 2 21 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint Figure 5 : Plots of the co-infection cases for individuals dually infected with anal cancer and syphilis in early and late stages of infection, respectively (Figure 5a and 5b) , when optimal vaccination control only strategy (φ = 0) and syphilis treatment controls for singly-infected (u 1 = u 2 = 0) are administered. Here, β s = 7.0, β h = 2.0. All other parameters as in Table 2 infected with persistent HPV and anal cancer or focusing on syphilis treatment among individuals dually infected with syphilis and anal cancer, will significantly curb the mixed infections. The simulations of the optimal control model (23) in the presence of combined optimal HPV vaccination strategy and syphilis treatment controls are depicted in Figures 14, 15 and 16 . It is observed in Figure 14 that the combined control stratgy has a high population level impact on the populations of individuals singly infected with syphilis in early stage (Figure 14 (a) ), individuals singly infected with syphilis in late stage (Figure 14 (b) ), individuals singly infected with HPV ( Figure 14 (c) ) and individuals singly infected with persistent HPV (Figure 14 (d) ). In a similar manner, the total number of individuals dually infected with HPV and syphilis in early and late stages of infection, respectively (Figures 15 (a) and 15 (b)), total number of individuals dually infected with persistent HPV and syphilis in early and late stages of infection (Figures 15 (c) and 15 (d) ) and total number of individuals dually infected with anal cancer and syphilis in early and late stages of infection, respectively (Figures 16 (a) and 16(b) ) all recorded lesser population number when the universal strategy is implemented than no no control is administered. The cost-effectiveness analysis is used to evaluate the health interventions related benefits so as to justify the costs of the strategies [5] . This is obtained by comparing the differences among the health outcomes and costs of those interventions; achieved by computing the incremental cost-effectiveness ratio (ICER), which is defined as the cost per health outcome. It is given by: Difference in costs between strategies Difference in health effects between strategies . We calculated the total number of co-infection cases averted and the total cost of the strategies applied in Table 3 . This is equally presented in Figure 18 . The total number of co-infection cases prevented is obtained by calculating the total number of individuals when controls are implemented and the total number when there is no control applied. Similarly, we apply the cost functions 1 2 w1φ 2 , 1 2 w2u 2 1 , 1 2 w3u 2 2 , 1 2 w4u 2 3 , 1 2 w5u 2 4 , 1 2 w6u 2 5 , 1 2 w7u 2 6 , 1 2 w8u 2 7 , 1 2 w9u 2 8 , over time, to compute the total cost for the various strategies implemented. We compared the cost-effectiveness of strategy I (Optimal HPV vaccination strategy for sexually active susceptible individuals) and strategy II (syphilis treatment controls for singly infected individuals). is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint Figure 13 : Plots of the co-infection cases for individuals dually infected with anal cancer and syphilis in early and late stages of infection, respectively (Figures 13(a) and 13b) . Here, β s = 7.0, β h = 2.0, φ = u 1 = u 2 = 0, u 3 = u 4 = u 5 = u 6 = u 7 = u 8 = 0. All other parameters as in Table 2 Figure 16: Plots of the co-infection cases for individuals dually infected anal cancer and syphilis in early and late stages of infection, respectively (Figures 16a and 16(b) , when there is universal control strategy and when there is no control administered. Here, β s = 7.0, β h = 2.0. All other parameters as in Table 2 is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint Comparing strategy II and strategy III, we observe that ICER (III) is greater than ICER (II), showing that strategy III strongly dominated strategy II and is more expensive and less effective compared to strategy II. Therefore, strategy III is removed from the list of next alternative strategies and we re-calculate ICER for the remaining competing strategies II and IV, as shown in Table 5 Comparing strategy II and strategy IV, it can be observed that ICER (IV) is greater than and strongly dominates ICER (II), showing that strategy IV is more costly and less effective compared to strategy III. As a result, strategy II (the strategy that implements syphilis treatment controls for singly infected individuals) has the least ICER and is the most cost-effective of all the control strategies for the control of HPV and syphilis co-infections. This is clearly illustrated in Figure 18 , which also agrees with the results obtained from both ACER and ICER methods that strategy II is the most cost-effective strategy. In this work, we have developed and presented a co-infection model for HPV and syphilis with cost-effectiveness optimal control analysis. The full co-infection model was shown to undergo the phenomenon of backward bifurcation when a certain condition was satisfied. The global asymptotic stability of the disease-free equilibrium of the full model was shown not to exist, when the associated reproduction number was less than unity. The existence of endemic equilibrium of the syphilis-only sub-model was shown to exist and the global asymptotic stability of the disease-free and endemic equilibria of both the syphilis-only sub-model and HPV-only sub-model were established. The global asymptotic stability of disease-free equilirbium of the HPV-only sub-model was also proven. Numerical simulations of the optimal control model showed that: i. HPV vaccination control has a positive population level impact in reducing the burden of HPV and the co-infection cases in a population ii. Syphilis treatment controls for singly infected individuals not only help bring down the burden of syphilis infection, but also reduce the burden of the HPV and syphilis co-infections. iii. The control strategy which implements syphilis treatment for singly infected individuals is the most cost-effective of all the control strategies in reducing the burden of HPV and syphilis co-infections. . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted September 11, 2020. . https://doi.org/10.1101/2020.09.09.20191635 doi: medRxiv preprint Calculus: A Complete Course Frequency of human papillomavirus (HPV) subtypes 31, 33, 35, 39 and 45 among Yemeni women with cervical cancer Analysis of a risk-structured vaccination model for the dynamics of oncogenic and warts-causing HPV types Optimal control of a two-strain tuberculosis-HIV/AIDS co-infection model Bifurcation analysis of a mathematical model for TB-Dengue co-infection Population dynamics of a mathematical model for TB-Dengue co-infection Maternal syphilis : pathophysiology and treatment Dynamical models of tuberculosis and their applications On the computation of R0 and its role on global stability CDC Vaccine Price List, Vaccines for Children Program (VFC) Correlates of homosexual behaviour and the incidence of anal cancer Population dynamics of a mathematical model for syphilis Mathematical analysis of reinfection and relapse in malaria dynamics The natural history of syphilis: implications for the transition dynamics and control of infection Host immunity and synchronized epidemics of syphilis across the United States The Cost and Cost-Effectiveness of Scaling up Screening and Treatment of Syphilis in Pregnancy: A Model The Stability of Dynamical Systems Optimal Control Applied to Biological Models The impact of an imperfect vaccine and pap cytology screening on the transmission of Human Papillomavirus and occurrence of associated cervical dysplasia and cancer Optimal control with multiple human papillomavirus vaccines A new mathematical model of syphilis Mathematical analysis of the transmission dynamics of HIV syphilis Co-infection in the presence of treatment for syphilis Risk factors for syphilis in young women attending a family health program in Vitória, Brazil Mathematical Assessment of the Role of Early Latent Infections and Targeted Control Strategies on Syphilis Transmission Dynamics Analysis of a mathematical model for COVID-19 population dynamics in Lagos Mathematical analysis of a two-sex Human Papillomavirus (HPV) model Analysis of a co-infection model for HPV-TB A mathematical study of a model for HPV with two high risk strains Analysis of COVID-19 and comorbidity co-infection Model with Optimal Control The Mathematical Theory of Optimal Processes A syphilis Co-Infection Study in Human Papilloma Virus Patients Attended in the Sexually Transmitted Infection Ambulatory Clinic Risk factors for anal cancer: results of a population-based case:control study, Cancer Causes and Control A mathematical model of syphilis transmission in an MSM population Optimal Control against the Human Papillomavirus: Protection versus Eradication of the Infection American Society for Colspcopy and Cervical Pathology, and American Society for Clinical pathology screening guidelines for the prevention and early detection of cervical cancer HPV-16/18, and Treponema pallidum Infections in a Sample of Brazilian Men Who Have Sex with Men Deterministic and Stochastic Models of the Dynamics of Drug Resistant Tuberculosis Mathematical Model and Optimal Control of New-Castle Disease (ND) Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission Use of the nonavalent HPV vaccine in individuals previously fully or partially vaccinated with bivalent or quadrivalent HPV vaccines Eliminating congenital syphilis. www.who.int/repro ducti ve-healt h/stis/syphilis.html Prevalence and Related Risk Behaviors of HIV, syphilis, and Anal HPV Infection Among Men who have Sex with Men from Beijing