Sensitivity analysis and optimal treatment control for a mathematical model of Human Papillomavirus infection

: Human papillomavirus (HPV) is one of the most common sexually transmitted viruses, and is a causal agent of cervical cancer. We aimed to develop a mathematical model of HPV natural history and qualitatively analyzed the stability of disease-free equilibrium, non-existence of limit cycle and existence of forward bifurcation. We performed sensitivity analysis to identify key epidemiological parameters. The Partial Rank Correlation Coefficient (PRCC) values for basic reproduction number shows that controlling contact rate plays an important role in disturbing equilibrium of HPV infection. Moreover, the increase of medical level is the most effective measure to prevent new HPV infections. Optimal treatment problem is solved and theoretical analysis is verified by numerical simulation.

relatively harmless, persistent infection with certain high-risk HPV genotypes can cause cervical precursor lesions and cervical carcinoma [1,2]. Cervical cancer is the most common malignancy reported among women in developing countries and is a leading cause of cancer deaths worldwide, with an estimated 233,000 deaths in 2000 [3]. By 2002, the number of women dying from cervical cancer had risen to 274,000 [4][5][6]. Nearly, all cervical cancers are attributable to genital infection with HPV [7,8]. Up to 70% of women in the general population will acquire genital HPV infection during their lifetime [9]. Most HPV related infections are cleared by treatment or natural healing of the host, while some are persistent and may develop into cervical cancer [10,11]. The progression from persistent infection to cervical cancer typically needs a period of 20 years or longer [12,13]. Fortunately, therapeutic HPV vaccines have made progress in clinical trials after long-term exploration and development. These vaccines include DNA vaccines, RNA replicon vaccines, peptide vaccines, vector vaccines, cell vaccines, and protein vaccines [14,15].
Mathematical models have been established to understand the dynamics of HPV related disease. In the modeling, some important variables specific for infectious disease are developed. For example, a variable representing how often an infectious individual contacts with a susceptible individual in a unit time is defined as the contact rate. The probability of each contact is set as  , indicating the ability of an infectious person to infect others (susceptible persons) in unit time. Considering that there are other members (e.g., immune persons   Rt), and that the number of members infected by one infected person in a unit time is limited, let  [16] formulated a mathematical model to investigate the impact of treatment and vaccination on HPV transmission dynamics. Elamin H and Elbasha et al [17] developed a nonlinear, deterministic and age-structured mathematical model. In the same year, Elamin H and Elbasha [18] studied the dynamic behaviors of a two-sex, deterministic model for assessing the potential impact of a prophylactic HPV vaccine with several properties. Mo'tassem and S Robert [19] developed an age-dependent two-sex mathematical model to describe the HPV vaccination program for a vaccine targeting HPV types 16 and 18 in both childhood and adult stages.
With the development of the mathematical model used to explain the spread of infectious diseases, optimization strategies are introduced into the model. Eihab B M et al [20] studied an optimal control model governed by a system of delay differential equations. The optimal vaccination strategy for a constrained time-varying SEIR (Susceptible, Exposed, Infected and Recovered) epidemic model was solved [21]. T K Kar and AshimBatabyal [22] formulated a nonlinear mathematical SIR epidemic model with a vaccination program. Yang and Tang et al [23] studied a mathematical model to explore the impact of vaccination and treatment on the transmission dynamics of tuberculosis (TB).
The purpose of this study was to incorporate both therapeutic HPV vaccines and partial immunity compartments based on the present models to investigate how they influence dynamic behavior and to further examine the optimal treatment control measures to minimize the number of asymptomatic patients and treatment costs. We formulate a mathematical model for HPV transmission, and the existence and stability of the equilibria are discussed. In the next section, we present a sensitivity analysis of the model through Partial Rank Correlation Coefficients to identify the key factors in the model. After that, we use optimal control and theory to minimize the number of asymptomatic patients and treatment costs; we also carry out some numerical simulations to highlight the effect of optimal treatment control on the model. Finally, we give concluding remarks.

Model formulation
Since persistent infection with high-risk HPV genotypes is the main cause of precursor cervical lesions and cervical carcinoma [24], we divided a cohort into three groups, namely, asymptomatic  Figure 1.
Combining all the aforementioned definitions and assumptions, it follows that the model for the transmission dynamics of HPV in a population is given by the following system of nonlinear differential equations: , , , , , , , , , , 0,

S E I I A R R S E I I A R S E I I
the region  is positively invariant. In this region, the model can be considered as been epidemiologically and mathematically well-posed [26]. Thus, each solution of the basic model (2.2) with initial conditions in  remains in  for all 0 t  .

Stability analysis
The model (2.2) has a disease-free equilibrium (DFE), given by , , , , , X E A I I S R  , using the notation [27], the model consists of nonnegative initial conditions together with the following system of equations: the matrices F and V , for the new infection terms and the remaining transfer terms, are, respectively, given by Consequently, we have the following conclusions Further, F is nonnegative, V is a nonsingular matrix and all eigenvalues of 4 J have positive real part. According to [27], all eigenvalues of FV  have negative real part. Therefore, the DFE of the model (2.2) is globally asymptotically stable if 0 1 R  .

Non-existence of limit cycles
Theorem 4. In the first quadrant, there is no limit cycle in system (2.2).
Proof We consider the Dulac function as   therefore, by the Dulac-Bendixson theorem [29], there is no periodic orbit for the system (2.2).

Existence of forward bifurcation
The epidemiological implications of backward bifurcation are that the effective disease control is only feasible if the basic reproduction number is reduced to values below another subthreshold less than unity [28]. Compared with backward bifurcation, the epidemiological significance of forward bifurcation means that the disease will disappear as long as the basic reproduction number is less than unity. Clearly, this phenomenon has important public health implications. Next, we analyze the existence and properties of bifurcation of the model. Let Setting the right-hand sides of model (2.2) to zero (at steady state) give by   * * * 1* 2* * * , , , , , It is instructive to characterize the type of bifurcation the model (2.2) may undergo, and we demonstrate the results in Appendix A. This phenomenon is illustrated by simulating model (2.2) with the set of parameter values listed in Table 2. The associated forward bifurcation diagram, depicted in Figure 2, shows that the model has a disease-free equilibrium (corresponding to 0 1 R  ) ( Figure 3) and an endemic equilibria (corresponding to 0 1 R  ) (Figure 4).

Efficacy of interventions
We used the parameter values given in Table 2 (unless otherwise stated) for our numerical simulations and changed the key parameter values to observe the corresponding effects on the disease outbreak.  To investigate the effect of contact rate on disease infection, we varied 1  , as shown in Figure 5.
The figure shows that as the contact rate of   Et with   St increases, the number of infected individuals increases, and consequently, the increased contact rate will result in an increase in the disease outbreak.
To examine the effect of treatment rate on disease infection, we plotted Figure 6. It demonstrates that increasing the parameter 1  and 2  decreases the endemic levels of asymptomatic patients and cancer-infected individuals. In particular, if we increased treatment rate 1  by 400% , the population of the asymptomatic patients and cancer-infected individuals would be decreased by 23% , and treatment rate 2  would be increased by 68% . Finally, the endemic levels of the cancer-infected individuals population would be decreased by 58% . Moreover, increasing treatment rate 2  by 81% would decrease the population of asymptomatic patients by 25% .

Sensitivity analysis of parameters in R 0
Latin hypercube sampling (LHS) is one of the Monte Carlo (MC) sampling methods that was first proposed by Mckay [34] in 1979. The advantage of LHS sampling is that the number of iterations is less than in other methods of random sampling, and the clustering phenomenon of sampling is avoided [35].
Model (2.2) has 14 parameters. In order to identify the factors associated with a given intervention that most strongly affect the spread of a new infection, following [36] we performed Latin Hypercube Sampling on the parameters that appear in 0 R . For the parameters in Eq. (2.7), Partial Rank Correlation Coefficients (PRCC) were calculated, and a total of 1000 simulations per LHS run were carried out. We chose a uniform distribution as the prior distribution when performing parameter sampling; the parameters of the model were set as input variables, and 0 R as the output variable. If the absolute value of the PRCC is larger, the influence of the parameter in 0 R is greater.
It is assumed that if the p value is greater than 0.05, the parameter is not significant for 0 R .

Sensitivity analysis of parameters in infectious variables
Time-varying sensitivity explains the influence of parameters on state variables when they change with time [35], thereby evaluating the influence and dependence of parameters on state variables in dynamic models. Given the different outbreak time periods of various diseases, one should pay attention to the effect of time period. For example, after being infected with HPV, the clinical manifestations are diverse, and individuals can even be asymptomatic. Therefore, the entire time period of disease prevention, outbreak, infection, and treatment should be considered.
Based on the pathogenesis of HPV, this section considers the sensitivity analysis of parameters in state variables from two aspects, namely, a continuous period of time and a single point in time.  In Figure 8, we show that the sensitivity of parameters in the early stages of disease outbreak changes significantly, especially before 500 t  . As can be seen from Figure 8 The influence of different parameters on state variables is also different at a certain time [37], so the influence of parameters on state variables at 4000 t  is studied, and the numerical simulation results can be seen in Table 4, Table 5, Figure 9, Figure 10, and Figure 11.       Table 4

 
Et is mainly affected by these parameters at this time.
From Figure 11 and Table 5, the monotonicity of parameter 1  , 2  , 1  , 2  , d , 1  and  is obvious, from Figure 9  has a positive correlation with   At, while 1  , 2  , d has a negative correlation with   At, which means that at this time, with the increase of 1  , 2  ,  , 1  , the number of   At will increase, but when 1  , 2  , d increases, the number of   At will decrease.

Optimal control of the extended model
The spread of infectious diseases can be controlled through reasonable treatment [16]. In this section, the optimal treatment problem will be addressed within the framework of the optimal control problem for constraints [21].

Formulation of the optimal treatment problem
We recorded 1  as u , and indicated the rate of people receiving treatment in   Et. Then an additional variable   Vt was introduced: where the value of   Vt at the initial time was set to 0 . According to the actual situation, the treatment capacity of a city or a country is limited, and the number of people who can be treated is limited as well [38], so the following formula can be established as is the time-dependent upper limit of the number of people that can be treated at each time instant. In addition, state restrictions can be introduced to keep the asymptomatic patients population at a low level, according to [39] we proposed the following constraints where max E is the upper limit of   Et.

T x S E I I A R 
, and the treatment rate   ut, was taken as the control variable. As to the cost function, the state variables and control variables are quadratic functions, which consider both the treatment rate and the number of asymptomatic patients. Based on the above description, the optimal treatment strategy can be formulated as the following problems with inequality constraints and free terminal state:

Analysis of optimal controls
The optimal control of the problem (P) is expressed as * u . Next, we use the maximum principle of Pontryagin [21] to derive the optimal system for problem (P). Prior to this, we converted the inequality constraints in the problem (P) into equality constraints with some non-negative parameter parameters    , with the transversality conditions From ( (4.12)

Numerical simulations
In this section, the symplectic pseudospectral method (SPM) developed in [40,41] is used to conduct numerical simulation. It owns structure-preserving property since the inherent Hamiltonian structure of optimal control problems is utilized. The local pseudospectral discretization scheme and the successive convexification technique are incorporated to make the algorithm fast and accurate. On one hand, the SPM can precisely quantify the energy variation for mechanical systems in optimal control and is thus widely used in vehicle trajectory planning [42][43][44]. On the other hand, it has sound stability for optimal control problems with long time span, which makes it an attractive method to solve policy-making problems for biological system [21].
When using the SPM in this section, the whole time span is divided into 10 regular sub-intervals and a 10-order approximation polynomial is used in each sub-interval. The parameters for the optimal control problem is set as Table 2 (unless otherwise stated in Table 6). The controlled solution is shown in Figure 12 along with the uncontrolled solution. It can be seen that the treatment strategy is effective, 1 I and A decreases rapidly and almost reaches zero at 40 t  , and the values of E and 2 I are also significantly reduced compared to before the addition of the treatment strategy. In addition, the values of S and R increased at the initial stage, and the effect before adding the treatment strategy is very obvious compared with the effect after the addition.
In the controlled case, the control variable   ut is continuously reduced and reaches zero at the terminal. It can be seen from Figure 12 (7) that the numerical solution of the control variable is in good agreement with the analytical solution given in (4.12), which verifies the characteristics of the optimal control.

Conclusion
For humans and some social animals, the standard incidence rate is more realistic than the bilinear incidence rate, so we used the standard incidence rate to indicate the ability of an infected person to infect others. Moreover, we assumed that the patients with persistent infection of HPV could be controlled and that mild infection could be cured, thus, those who are cured or have self-healed will gain immunity. Next, we established the model and analyzed its properties. It can be seen from the analysis results that controlling the value of the parameter in 0 R , making 0 R less than 1, the disease will be die out, and when 0 R greater than 1, the infection will experience an outbreak and form an endemic disease. The sensitivity analysis of the parameters in 0 R and state variables illustrates that the value of 0 R can be reduced or increased more efficiently, and the results can be used to control the spread of HPV. In the last section, the treatment rate is selected as the control variable, and the optimal treatment strategy is analyzed theoretically and simulated numerically. The results show that the appropriate treatment strategy can efficiently reduce the spread of the disease.