# Modeling and Mathematical Analysis of the Dynamics of HPV in Cervical Epithelial Cells: Transient, Acute, Latency, and Chronic Infections

Authors: Juan Carlos Sierra-Rojas; Ramón Reyes-Carreto; Cruz Vargas-De-León; Jorge Fernando Camacho
Journal: Computational and Mathematical Methods in Medicine
 (2022)

---

## Abstract

The aim of this paper is to model the dynamics of the human papillomavirus (HPV) in cervical epithelial cells. We developed a mathematical model of the epithelial cellular dynamics of the stratified epithelium of three (basale, intermedium, and corneum) stratums that is based on three ordinary differential equations. We determine the biological condition for the existence of the epithelial cell homeostasis equilibrium, and we obtain the necessary and sufficient conditions for its global stability using the method of Lyapunov functions and a theorem on limiting systems. We have also developed a mathematical model based on seven ordinary differential equations that describes the dynamics of HPV infection. We calculated the basic reproductive number (R0) of the infection using the next-generation operator method. We determine the existence and the local stability of the equilibrium point of the cellular homeostasis of the epithelium. We then give a sufficient condition for the global asymptotic stability of the epithelial cell homeostasis equilibrium using the Lyapunov function method. We proved that this equilibrium point is nonhyperbolic when R0=1 and that in this case, the system presents a forward bifurcation, which shows the existence of an infected equilibrium point when R0>1. We also study the solutions numerically (i.e., viral kinetic in silico) when R0>1. Finally, local sensitivity index was calculated to assess the influence of different parameters on basic reproductive number. Our model reproduces the transient, acute, latent, and chronic infections that have been reported in studies of the natural history of HPV.


---

## Body

## 1. Introduction

Human papillomavirus (HPV) is small, nonenveloped, icosahedral DNA viruses that have a diameter of 52-55 nm. HPV is one of the most common sexually transmitted diseases in the world, and it is the principal causative agent of cervical cancer (CC), which occurs in 99.7% of cases [1]. This is a large family of small viruses that are classified into low- and high-risk (HR) genotypes and which can cause abnormal cell proliferation, manifesting from epithelial warts to high-grade cervical intraepithelial neoplasia (CIN) [2].At the cellular level, HPV infects keratinocytes (i.e., cells that are the most predominant in the epidermis). Traditionally, the epidermis is segmented into distinct structural and functional compartments, which are called the cornified layer (stratum corneum), granular layer (stratum granulosum), spinous layer (stratum spinosum), and the deepest layer the basal layer (stratum basale) (see Figure1). Keratinocytes differentiate as they move through the cell layers, starting as basal keratinocytes. The keratinocytes produce more and more keratin, and they eventually undergo natural cell death and detachment (anoikis). The cornified keratinocytes that form the outermost layer of the epidermis are constantly shed off and replaced by new cells [3]. The differentiation can be lateral and suprabasal: in the first, other cells of the same stratum are produced; and in the second, other cells that change stratum are produced [4]. When a cell is infected by HPV, it retains this capacity for cell differentiation.Figure 1
Cervical stratified squamous epithelial cell architecture.The infection is produced by the exposure of the stratum basale cells together with a HPV particle when microtraumas occur in the epithelium during sexual activity [5, 6]. The virus only infects stratum basale cells because they contain the receptors that allow binding with the L1 protein in the capsid of the virus [7, 8].During the process of infection by HR genotypes, there is no viremic phase, replication is nonlytic, and levels of viral gene expression are kept low in the basal epithelium. This limits the innate immune response and adaptive is delayed, which favors the establishment of viral infection. Furthermore, the E6 and E7 oncoproteins interfere with the activation levels of type I interferon in the infected cell, which prevents the initiation of intracellular antiviral responses [9, 10].The HPV replicative cycle can last between 6 and 12 weeks, which considerably expands the viral genome. This results in new mature viral particles that can reach 1,000, which are released when rupture of the cell membrane occurs [11].As illustrated in Figure2, the release of viral loads from infected cells produces viral kinetics that allows us to classify the infection as transient, acute, latent, or chronic. The first is transient when viral genetic material is present in the host for a short period, and there is no infection in the cells of the stratum basale. The viral load is removed until its complete elimination during the following days. The second is acute when the infection is presented, replicates the viral genome, produces viral particles, and is eventually cleared. The third is latent when acute infections only appear to clear, but the viral genome remains in the infected cell without detectable activity. Finally, the fourth is chronic if genetic material is not eliminated during the acute infection but continues to increase the viral activity over time [12, 13].Figure 2
Theoretical HPV viral load kinetics. Figure adapted from Alizon et al. [12].The persistence of HR HPV infection is one of the principal causes of CIN I, II, or III and cancer. Although it is not clear whether the viral load has a causal effect on increasing the risk of cervical lesions and cancer in HPV-infected women, a high viral load is associated with the persistence of HPV infection. Thus, viral load may be a marker of HPV persistence [14].A retrospective natural history study of HPV infections through serial HPV viral load measurement in 261 untreated women with type-specific HPV DNA detected has shown a regression or decrease of a clonal cell population HPV-infected when the infection was latent [13]. This indicates that the kinetics of viral load can illustrate predictive scenarios on regression or progression of the type of viral infection presented by a HR HPV [13]. A protocol for a cohort study, called PAPCLEAR, has recently been reported, which was aimed at better understanding the course and natural history of cervical HPV infections in healthy, unvaccinated, and vaccinated young women [15]. This study will be relevant because of its impact in the clinic thanks to the possible integration of longitudinal data to mathematical models.A lot of literature has been published [16–20] on the mathematical modeling of viral infections, typically in human immunodeficiency virus, hepatitis B and C, and influenza. Although the literature reports studies on HPV modeling at the cellular level and viral kinetic scenarios [21–24], they do not show a direct relationship with the types of infection considering the different stratum of the squamous epithelium to show progression or regression of an infection by the HPV. Asih et al. [21] proposed a model that considers four compartments: susceptible cells, infected cells, precancerous cells, and cancer cells. They analyzed the local stability of the equilibrium points of the model and investigated the parameters that play an important role in the progression toward invasive cancer. Akimenko and Adi-Kusumo [22], and Sari and Adi-Kusumo [23] proposed two models of an age-structured subpopulations of susceptible, infected, precancerous, and cancer cells and unstructured subpopulation of HPV that aimed at gaining insight into the disease characteristics of cervical cancer. Verma et al. [24] developed a mathematical model of HIV/HPV coinfection of the oral mucosa. This model considered the cellular immune response and the basal and suprabasal layers of epithelial tissue but ignores viral transmission via suprabasal differentiation. They obtained simulations of the kinetics of an acute infection that tends to disappear over time and the kinetics of a chronic infection. Finally, Murall et al. [25] propose a viral dynamics model of the stratified epithelium of four layers. They simulated a scenario of a slow growing HR HPV infection that spontaneously regresses and another scenario where the infection is inoculated with few cells and the microabrasion repairs quickly. However, this mathematical model ignores the dynamics of the healthy stratified epithelium, and although it models the viral transmission via suprabasal differentiation, it does so with a linear term.Motivated by the above, we propose a deterministic model that includes several layers of the healthy and infected squamous epithelium without the presence of the immune system to simulatein silico the types of infection reported in the literature when HPV infection occurs.The structure of this manuscript is as follows. In Section2, we will describe a model to epithelial cellular dynamics of healthy tissue, and we will study the stability of its equilibrium points. In Section 3, we will describe the dynamics of HPV in the stratified epithelium of the cervix, and we will calculate the basic reproductive number for the viral infection, its equilibrium points, and their corresponding local and global stabilities. Additionally, we will show that one of these equilibrium points, the infection-free equilibrium, can be nonhyperbolic and that in this case, there is a bifurcation; this result shows the existence of an infected equilibrium point. Section 4 shows the typical values required in the model that describes the dynamics of HPV-infected tissue, and with them, it is estimated numerically the local bifurcation type that this model exhibits, the regions of existence and stability of its infected equilibrium point, and the simulation of different scenarios of interest that reproduce the transient, acute, latent, and chronic infections of the natural history of the HPV. In Section 5, we will perform a local sensitivity analysis of the basic reproductive number. This will be followed by a discussion of the results obtained and the conclusions, indicated in Sections 6 and 7, respectively.


## 2. Epithelial Cellular Dynamics

We assume that stratum spinosum and stratum granulosum cells have similar cell dynamics. Therefore, these cells belong to a stratum that we call stratum intermedium (see Figure1). Consequently, the cell dynamic model of the homeostatic stratified epithelium is composed of the stratum basale, stratum intermedium, and stratum corneum. The dynamic is visualized in Figure 3.Figure 3
Diagram of the dynamics of the homeostatic stratified epithelium.The model is given byBHt, EHt, and CHt, which denote stratum basale, stratum intermedium, and stratum corneum formed by uninfected cells, respectively. The cells of BHt proliferate (or lateral differentiation) at rate rB considering a logistic growth, with a carrying capacity KB; these have suprabasal differentiation to EHt at rate δB and die at rate μB. The dynamics of EHt is generated by suprabasal differentiation of BHt at a rate δB. This population increases by cell proliferation, which is governed by logistic growth at a rate rE and a carrying capacity at rate KE. These decrease by suprabasal differentiation at a rate δE, and we assume that there is no natural death. Finally, the dynamics of CHt are generated by suprabasal differentiation of EHt at rate δE and shed from the epithelial tissue at a rate μC. We arrive at the following mathematical model:
(1)dBHdt=rBBH1−BHKB−δBBH−μBBH,dEHdt=δBBH+rEEH1−EHKE−δEEH,dCHdt=δEEH−μCCH.


### 2.1. Positivity, Boundedness, Equilibria, and Their Local Stabilities

Let setΩP be ΩP=BH,EH,CH∈ℝ+3:BH≥0,EH≥0,CH≥0⊂ℝ3. We consider the positivity and boundedness of the solutions of system (1).Theorem 1.
Given the initial conditionsBH0>0,EH0>0,CH0>0, then all solutions of system (1) are positive.Proof.
At first, we will prove the positivity ofBHt. Let BH0 be the solution that satisfies the initial condition BHt>0. Assume that the solution is not always positive; i.e., there exists a t0′∈ℝ+, such as BHt0′<0. By Bolzano’s theorem, there exists a t1∈0,t0′, such as BHt1=0. Let t0∈ℝ+ be the initial time, such as BHt0=0, and then dBHt/dt=0. Note that if some t≥0, BHt=0, then dBHt/dt=0. Then, any solution with BH0=0 will satisfy BHt=0∀t>0. By uniqueness of solutions, we have that if BH0>0, then BHt will remain positive ∀t>0. Therefore, BHt0=0 leads to a contradiction. Hence, BHt is nonnegative for all t>0. Now, we will prove the positivity of EHt. EHt is not always positive; i.e., there exists t0′∈ℝ+, such as EHt0′<0. By Bolzano’s theorem, there exists t1∈0,t0′, such as EHt1=0. Let t0=minti∣EHti=0. By the second equation (1), if EHt0=0, then dEHt0/dt=δbBH>0 implies that EH is increasing at t=t0. Therefore, EHt will be negative for values t<t0 near to t0, that is, a contradiction. Finally, we will prove the positivity of CH. From the third equation (1), we obtain the following inequality dCHt/dt≥−μCCHt. Integrating, we obtain the solution CH≥CH0e−μt; therefore, CH≥0.Theorem 2.
LetBHt,EHt,CHt be the solution of model (1) with the initial conditions BH0>0, EH0>0, and CH0>0. Then, BHt, EHt, and CHt are all bounded for all t≥0 at which the solution exists.Proof.
Let (BHt, EHt, CHt) be any solution with nonnegative initial conditions. We define a function
(2)Bt=BHt+EHt+n−1nCHt,n≫1.
The time derivative along a solution of (2) is
(3)dBtdt=rBBH1−BHKB−μBBH+rEEH1−EHKE−δEnEH−μCn−1nCH=rBKB4+rEKE4−μBBH−δEnEH−μCn−1nCH−rBKBBH−KB22−rEKEEH−KE22≤rBKB4+rEKE4−μBBH−δEnEH−μCn−1nCH.
It follows that(4)dBtdt+ηBt≤rBKB+rEKE4,where η=minμB,δE/n,μC. Thus, limsupt⟶∞Bt≤rBKB+rEKE/4η. Therefore, BHt, EHt, and CHt are all bounded for all t≥0.Remark 3.
The total number of cervical epithelial cells is bounded by a weighted sum of the stratum basale and stratum intermedium carrying capacities, where the weights are the proliferation rates divided by four times of minimum of the death and differentiation rates. The above is biologically plausible that there is an upper bound in terms of the carrying capacities of the cells.
We can easily see that for all of the parameter values, the trivial equilibriumE01=0,0,0 always exists. We get the “epithelial cell homeostasis” equilibrium E11=BH∗,EH∗,CH∗, where
(5)BH∗=KB1−δB+μBrB,(6)EH∗=KE2rErE−δE+δE−rE2+4δBrEKBKE1−δB+μBrB,(7)CH∗=δEμCEH∗.The following inequalityrB>δB+μB is a biological condition for the maintenance of cell homeostasis of the epithelium.The Jacobian matrix of (1) at E01 is
(8)JE01=rB−δB−μB00δBrE−δE00δE−μC.The eigenvalues ofJE01 are λ1=rB−δB−μB, λ2=rE−δE, and λ3=−μC. By biological condition, rB>δB+μB, then JE01 has one positive eigenvalue. Thus, the trivial equilibrium is unstable.The Jacobian matrix for the equations of (1) at E11 is
(9)JE11=rB1−2BHKB−δB−μB00δBrE−2rEEHKE−δE00δE−μC.Using the following identitiesrB−rBBH∗/KB−δB+μB=0 and rE−rEEH∗/KE−δE=−δBBH∗/EH∗, we have
(10)JE11=−BH∗KB00δB−δBBH∗EH∗−rEEH∗KE00δE−μC.The characteristic polynomial ofJE11 is
(11)Pλ11=λ11−−BH∗KBλ11−−δBBH∗EH∗−rEEH∗KEλ11−−μC.The eigenvalues ofPλ11 are λ1,11=−BH∗/KB, λ1,21=−δBBH∗/EH∗−rEEH∗/KE, and λ1,31=−μC. Clearly, all of the eigenvalues are negative. Consequently, the epithelial cell homeostasis equilibrium E11 is locally asymptotically stable.We then arrive at the following theorem.Theorem 4.
Assume that the biological conditionrB>δB+μB is satisfied. The trivial equilibrium E01=0,0,0 always exists and is unstable. The epithelial cell homeostasis equilibrium E11=BH∗,EH∗,CH∗ always exists and is absolutely stable.


### 2.2. Global Stability of the Epithelial Cell Homeostasis Equilibrium

We use the method of Lyapunov functions and a theorem on limiting systems to analyze the global stability of the epithelial cell homeostasis equilibrium of the system (1).Consider theBH−EH subsystem of model (1), which is independent of the CH variables. We construct the following Volterra-type Lyapunov function [26] for the BH−EH subsystem
(12)Ut=2∫BH∗BH1−BH∗ηdη+rBδBKB∫EH∗EH1−EH∗ηdη.The functionUt is defined, continuous, and positive definite for all BH,EH>0. Also, the global minimum UBH,EH=0 occurs at BH∗,EH∗, and therefore, U is a Lyapunov function. First, we calculate the time derivative of Ut computed along solutions of the first two equations of the system (1), given by the expression
(13)dUtdt=21−BH∗BHdBHdt+rBδBKB1−EH∗EHdEHdt.Note that(14)1=BH∗KB+δB+μBrB,(15)1=δErE−δBrEBH∗EH∗+EH∗KE.Using (14) and (15), we have
(16)dUtdt=1−BH∗BHrBBH∗KBBH1−BHBH∗+1−BH∗BHrBBH∗KBBH1−BHBH∗+rBδBKB1−EE∗EHδBBH∗BHBH∗−EHEH∗+rEEHEH∗KE1−EHEH∗=rBBHBH∗KB2−BHBH∗−BH∗BH+rBδBKBrEEHEH∗KE2−EHEH∗−EH∗EH+rBBH∗KB3−BH∗BH−EHEH∗−BHEE∗BH∗EH<−rBKBBH−BH∗2−rBδBKBrEKEEH−EH∗2.BecausedUt/dt≤0, the Lyapunov stability theorem [27] implies that the BH∗,EH∗ equilibrium is globally asymptotically stable in ℝ+2 and limt⟶∞EHt=EH∗.From the third equation in (1), we can formally solve to obtain
(17)CHt=CH0+∫t0tδEEHτeμCτ−t0dτ/eμCt−t0.By L’Hospital’s rule, we obtainlimt⟶∞CHt=δEEH/μC=δEEH∗/μC=CH∗. An application of Lemma 1 [28] shows that the epithelial cell homeostasis equilibrium E11=BH∗,EH∗,CH∗ of model (1) is globally asymptotically stable in the interior of ΩP. We then have the following theorem.Theorem 5.
IfrB>δB+μB, then the epithelial cell homeostasis equilibrium E11=BH∗,EH∗,CH∗ of system (1) is globally asymptotically stable in the interior of ΩP.


## 2.1. Positivity, Boundedness, Equilibria, and Their Local Stabilities

Let setΩP be ΩP=BH,EH,CH∈ℝ+3:BH≥0,EH≥0,CH≥0⊂ℝ3. We consider the positivity and boundedness of the solutions of system (1).Theorem 1.
Given the initial conditionsBH0>0,EH0>0,CH0>0, then all solutions of system (1) are positive.Proof.
At first, we will prove the positivity ofBHt. Let BH0 be the solution that satisfies the initial condition BHt>0. Assume that the solution is not always positive; i.e., there exists a t0′∈ℝ+, such as BHt0′<0. By Bolzano’s theorem, there exists a t1∈0,t0′, such as BHt1=0. Let t0∈ℝ+ be the initial time, such as BHt0=0, and then dBHt/dt=0. Note that if some t≥0, BHt=0, then dBHt/dt=0. Then, any solution with BH0=0 will satisfy BHt=0∀t>0. By uniqueness of solutions, we have that if BH0>0, then BHt will remain positive ∀t>0. Therefore, BHt0=0 leads to a contradiction. Hence, BHt is nonnegative for all t>0. Now, we will prove the positivity of EHt. EHt is not always positive; i.e., there exists t0′∈ℝ+, such as EHt0′<0. By Bolzano’s theorem, there exists t1∈0,t0′, such as EHt1=0. Let t0=minti∣EHti=0. By the second equation (1), if EHt0=0, then dEHt0/dt=δbBH>0 implies that EH is increasing at t=t0. Therefore, EHt will be negative for values t<t0 near to t0, that is, a contradiction. Finally, we will prove the positivity of CH. From the third equation (1), we obtain the following inequality dCHt/dt≥−μCCHt. Integrating, we obtain the solution CH≥CH0e−μt; therefore, CH≥0.Theorem 2.
LetBHt,EHt,CHt be the solution of model (1) with the initial conditions BH0>0, EH0>0, and CH0>0. Then, BHt, EHt, and CHt are all bounded for all t≥0 at which the solution exists.Proof.
Let (BHt, EHt, CHt) be any solution with nonnegative initial conditions. We define a function
(2)Bt=BHt+EHt+n−1nCHt,n≫1.
The time derivative along a solution of (2) is
(3)dBtdt=rBBH1−BHKB−μBBH+rEEH1−EHKE−δEnEH−μCn−1nCH=rBKB4+rEKE4−μBBH−δEnEH−μCn−1nCH−rBKBBH−KB22−rEKEEH−KE22≤rBKB4+rEKE4−μBBH−δEnEH−μCn−1nCH.
It follows that(4)dBtdt+ηBt≤rBKB+rEKE4,where η=minμB,δE/n,μC. Thus, limsupt⟶∞Bt≤rBKB+rEKE/4η. Therefore, BHt, EHt, and CHt are all bounded for all t≥0.Remark 3.
The total number of cervical epithelial cells is bounded by a weighted sum of the stratum basale and stratum intermedium carrying capacities, where the weights are the proliferation rates divided by four times of minimum of the death and differentiation rates. The above is biologically plausible that there is an upper bound in terms of the carrying capacities of the cells.
We can easily see that for all of the parameter values, the trivial equilibriumE01=0,0,0 always exists. We get the “epithelial cell homeostasis” equilibrium E11=BH∗,EH∗,CH∗, where
(5)BH∗=KB1−δB+μBrB,(6)EH∗=KE2rErE−δE+δE−rE2+4δBrEKBKE1−δB+μBrB,(7)CH∗=δEμCEH∗.The following inequalityrB>δB+μB is a biological condition for the maintenance of cell homeostasis of the epithelium.The Jacobian matrix of (1) at E01 is
(8)JE01=rB−δB−μB00δBrE−δE00δE−μC.The eigenvalues ofJE01 are λ1=rB−δB−μB, λ2=rE−δE, and λ3=−μC. By biological condition, rB>δB+μB, then JE01 has one positive eigenvalue. Thus, the trivial equilibrium is unstable.The Jacobian matrix for the equations of (1) at E11 is
(9)JE11=rB1−2BHKB−δB−μB00δBrE−2rEEHKE−δE00δE−μC.Using the following identitiesrB−rBBH∗/KB−δB+μB=0 and rE−rEEH∗/KE−δE=−δBBH∗/EH∗, we have
(10)JE11=−BH∗KB00δB−δBBH∗EH∗−rEEH∗KE00δE−μC.The characteristic polynomial ofJE11 is
(11)Pλ11=λ11−−BH∗KBλ11−−δBBH∗EH∗−rEEH∗KEλ11−−μC.The eigenvalues ofPλ11 are λ1,11=−BH∗/KB, λ1,21=−δBBH∗/EH∗−rEEH∗/KE, and λ1,31=−μC. Clearly, all of the eigenvalues are negative. Consequently, the epithelial cell homeostasis equilibrium E11 is locally asymptotically stable.We then arrive at the following theorem.Theorem 4.
Assume that the biological conditionrB>δB+μB is satisfied. The trivial equilibrium E01=0,0,0 always exists and is unstable. The epithelial cell homeostasis equilibrium E11=BH∗,EH∗,CH∗ always exists and is absolutely stable.


## 2.2. Global Stability of the Epithelial Cell Homeostasis Equilibrium

We use the method of Lyapunov functions and a theorem on limiting systems to analyze the global stability of the epithelial cell homeostasis equilibrium of the system (1).Consider theBH−EH subsystem of model (1), which is independent of the CH variables. We construct the following Volterra-type Lyapunov function [26] for the BH−EH subsystem
(12)Ut=2∫BH∗BH1−BH∗ηdη+rBδBKB∫EH∗EH1−EH∗ηdη.The functionUt is defined, continuous, and positive definite for all BH,EH>0. Also, the global minimum UBH,EH=0 occurs at BH∗,EH∗, and therefore, U is a Lyapunov function. First, we calculate the time derivative of Ut computed along solutions of the first two equations of the system (1), given by the expression
(13)dUtdt=21−BH∗BHdBHdt+rBδBKB1−EH∗EHdEHdt.Note that(14)1=BH∗KB+δB+μBrB,(15)1=δErE−δBrEBH∗EH∗+EH∗KE.Using (14) and (15), we have
(16)dUtdt=1−BH∗BHrBBH∗KBBH1−BHBH∗+1−BH∗BHrBBH∗KBBH1−BHBH∗+rBδBKB1−EE∗EHδBBH∗BHBH∗−EHEH∗+rEEHEH∗KE1−EHEH∗=rBBHBH∗KB2−BHBH∗−BH∗BH+rBδBKBrEEHEH∗KE2−EHEH∗−EH∗EH+rBBH∗KB3−BH∗BH−EHEH∗−BHEE∗BH∗EH<−rBKBBH−BH∗2−rBδBKBrEKEEH−EH∗2.BecausedUt/dt≤0, the Lyapunov stability theorem [27] implies that the BH∗,EH∗ equilibrium is globally asymptotically stable in ℝ+2 and limt⟶∞EHt=EH∗.From the third equation in (1), we can formally solve to obtain
(17)CHt=CH0+∫t0tδEEHτeμCτ−t0dτ/eμCt−t0.By L’Hospital’s rule, we obtainlimt⟶∞CHt=δEEH/μC=δEEH∗/μC=CH∗. An application of Lemma 1 [28] shows that the epithelial cell homeostasis equilibrium E11=BH∗,EH∗,CH∗ of model (1) is globally asymptotically stable in the interior of ΩP. We then have the following theorem.Theorem 5.
IfrB>δB+μB, then the epithelial cell homeostasis equilibrium E11=BH∗,EH∗,CH∗ of system (1) is globally asymptotically stable in the interior of ΩP.


## 3. Epithelial Viral Dynamics of HPV

The viral dynamics model includes that of uninfected cells, infected cells, and viral load. The viral dynamics are given byBIt, EIt, CIt, and Vt, which denote the cells of the stratum basale, stratum intermedium, stratum corneum infected, and the viral load, respectively. BHt, EHt, and CHt are denoted in the same way as in the previous section.Figure4 shows that the dynamics of BIt results from contact between uninfected basal cells and HPV particle at rate β. This population increases by cell proliferation (or lateral differentiation), governed by full logistic growth at a rate rB∗, and a carrying capacity at rate KB. These decrease by suprabasal differentiation at rate δB∗ and death cellular at rate μB∗. The dynamics of EIt are generated by suprabasal differentiation of BIt at a rate δB∗. This population increases by cell proliferation, governed by full logistic growth at a rate rE∗, and a carrying capacity at rate KE. These decrease by suprabasal differentiation at a rate δE∗, and we assume that there is no natural death. The dynamics of CIt are generated by suprabasal differentiation of EIt at rate δE∗. There is no proliferation of CIt. This population decreases by desquamation at rate μC∗. Finally, Vt increases by rupture of cell membrane of CIt at rate σ, and they decline at rate γ. This is summarized in the following nonlinear ODE system:
(18)dBHdt=rBBH1−BH+BIKB−δBBH−μBBH−βBHV,dBIdt=βBHV+rB∗BI1−BH+BIKB−δB∗BI−μB∗BI,dEHdt=δBBH+rEEH1−EH+EIKE−δEEH,dEIdt=δB∗BI+rE∗EI1−EH+EIKE−δE∗EI,dCHdt=δEEH−μCCH,dCIdt=δE∗EI−μC∗CI,dVdt=σCI−γV.Figure 4
Diagram of viral dynamics when a HPV infection occurs.


### 3.1. Positivity, Boundedness, Equilibria, Basic Reproductive Number, and Local Stability

Let setΩG be ΩG=BH,BI,EH,EI,CH,CI,V∈ℝ+7:BH≥0,BI≥0,EH≥0,EI≥0,CH≥0,CI≥0,V≥0⊂ℝ7. We consider the positivity and boundedness of the solutions of system (18).Theorem 6.
Given the initial conditionsBH0>0, BI0≥0, EH0>0, EI0≥0, CH0>0, CH0≥0 and V0>0, then all solutions of system (18) are positive.The proof of Theorem1 is very similar to the proof of Theorem 6. For this reason, we omit the proof.Theorem 7.
LetBHt,BIt,EHt,EIt,CHt,CIt,Vt be the solution of model (18) with the initial conditions BH0>0, BI0≥0, EH0>0, EI0≥0, CH0>0, CH0≥0 and V0>0. Then, BHt, BIt, EHt, EIt, CHt, CIt, and Vt are all bounded for all t≥0 at which the solution exists.Proof.
LetBHt,BIt,EHt,EIt,CHt,CIt,Vt be any solution with nonnegative initial conditions. We define a function
(19)Bt=BHt+BIt+EHt+EIt+nn+1CHt+nn+1CIt+μC∗n−1σn+1Vt,n≫1.
The time derivative along a solution of (19) is
(20)dBtdt=rBBH1−BHKB+rB∗BI1−BIKB+rEEH1−EHKE+rE∗EI1−EIKE−rB+rB∗KBBHBI−rE+rE∗KEEHEI−μBBH−μB∗BI−δEn+1EH−δE∗n+1EI−μCnn+1CH−μC∗nnn+1CI−γμC∗n−1σn+1Vt=rB+rB∗KB4+rE+rE∗KE4−rBKBBH−KB22−rB∗KBBI−KB22−rEKEEH−KE22−rE∗KEEI−KE22−rB+rB∗KBBHBI−rE+rE∗KEEHEI−μBBH−μB∗BI−δEn+1EH−δE∗n+1EI−μCnn+1CH−μC∗nnn+1CI−γμC∗n−1σn+1Vt≤rB+rB∗KB+rE+rE∗KE4−μBBH−μB∗BI−δEn+1EH−δE∗n+1EI−μCnn+1CH−μC∗nnn+1CI−γμC∗n−1σn+1Vt.
It follows that(21)dBtdt+ηBt≤rB+rB∗KB+rE+rE∗KE4,where η=minμB,μB∗,δE/n+1,δE∗/n+1,μC,μC∗/n,γ. Thus, limsupt⟶∞Bt≤rB+rB∗KB+rE+rE∗KE/4η. Therefore, BHt, BIt, EHt, EIt, CHt, CIt, and Vt are all bounded for all t≥0.Remark 8.
The total number of cervical epithelial cells, both uninfected and infected, and viral load is bounded by a weighted sum of the stratum basale and stratum intermedium carrying capacities, where the weights are the proliferation rates of healthy and infected cells divided by four times of minimum of the death, differentiation, and viral clearance rates.The model (18) always have a trivial equilibrium E02=0,0,0,0,0,0,0, and a “epithelial cell homeostasis” equilibrium E12=BH∗,0,EH∗,0,CH∗,0,0 with the same coordinates BH∗, EH∗, and CH∗ given in (5), (6) and (7), respectively. We recall the “epithelial cell homeostasis” equilibrium as the infection-free equilibrium.The Jacobian matrix of system (18) is
(22)J=J11−rBBHKB0000−βBHβV−rB∗BIKBJ220000βBHδB0J33−rEEHKE0000δB∗−rE∗EIKEJ4400000δE0−μC00000δE∗0−μC∗000000σ−γ,where
(23)J11=rB−2rBKBBH−rBKBBI−δB−μB−βV,J22=rB∗−rB∗KBBH−2rB∗KBBI−δB∗−μB∗,J33=rE−2rEKEEH−rEKEEI−δE,J44=rE∗−rE∗KEEH−2rE∗KEEI−δE∗.


#### 3.1.1. Basic Reproductive NumberR0

To compute the basic reproductive numberR0, we use the next-generation operator introduced by van den Driessche and Watmough [29]. The Jacobian matrix J of this subsystem at the infection-free equilibrium is decomposed as J=F−V, where F is the viral transmission part and V describe the transition terms associated with the model (18). These quantities are given, respectively, by
(24)F≔rB∗1−BH∗KB00βBH∗000000000000,V≔δB∗+μB∗000−δB∗δE∗−rE∗1−EH∗KE000−δE∗μC∗000−σγ.It follows that the basic reproductive numberR0=ρFV−1, where ρ is the spectral radius, is given by
(25)R0=1δB∗+μB∗βσδB∗δE∗BH∗γμC∗δE∗−rE∗1−EH∗/KE+rB∗1−BH∗KB.The parameterR0 has an interesting biological meaning: it is the sum of average numbers of secondary infected cells produced by a single infected cell in a population of epithelial basal cells, by direct basal cell-to-HPV contact and infected basal cell division, respectively.Remark 9.
IfrE∗1−EH∗/KE is taken as a new infection, we can get another basic reproductive number as follows:
(26)R¯0=12a+b+c+12a2+b−c2+2ab+c,where
(27)a=βσδB∗BH∗γμC∗δB∗+μB∗,b=rB∗δB∗+μB∗1−BH∗KB,c=rE∗δE∗1−EH∗KE.
WhenrE∗=0, it is easy to check that R0 is equivalent to R¯0.


#### 3.1.2. Local Stability ofE02

The Jacobian matrix (22) evaluated at the equilibrium point E02 becomes
(28)JE02=rB−δB−μB0000000rB∗−δB∗−μB∗00000δB0rE−δE00000δB∗0rE∗−δE∗00000δB0−μC00000δB∗0−μC∗000000σ−γ.By the biological condition, one of its eigenvalues is positive:(29)λ0,12=rB−δB−μB>0,while the other six are
(30)λ0,22=rE−δE,λ0,32=rB∗−δB∗−μB∗,λ0,42=rE∗−δE∗,λ0,52=−μC,λ0,62=−μC∗,λ0,72=−γ.Thus,E02 is unstable. This result can be summarized as follows.Theorem 10.
Assume that the biological conditionrB>δB+μB is satisfied. The trivial equilibrium E02=0,0,0,0,0,0,0 always exists and is unstable.


#### 3.1.3. Local Stability ofE12

The Jacobian matrix of system (18), evaluated in the infection-free equilibrium point E12, takes the form
(31)JE12=J11−rBKBBH∗0000−βBH∗0J220000βBH∗δB0J33−rEKEEH∗0000δB∗0J4400000δE0−μC00000δE∗0−μC∗000000σ−γwhere
(32)J11=rB−δB−μB−2rBKBBH∗,J22=rB∗−δB∗−μB∗−rB∗KBBH∗,J33=rE−δE−2rEKEEH∗,J44=rE∗−δE∗−rE∗KEEH∗.Using the following identitiesrB−rBBH∗/KB−δB+μB=0 and rE−rEEH∗/KE−δE=−δBBH∗/EH∗, we have
(33)J11=−rB−δB−μB,(34)J22=rB∗−δB∗−μB∗−rB∗rBrB−δB−μB,(35)J33=−δBBH∗EH∗−rEKEEH∗,(36)J44=rE∗−δE∗−rE∗rErE−δE−rE∗rEδBBH∗EH∗.By biological condition,rB>δB+μB, then J11<0 and J33<0. Furthermore, if
(37)rB∗δB∗+μB∗≤rBδB+μB,(38)rE∗δE∗<rEδE,thenJ22≤0 and J44<0, respectively.On the other hand, the characteristic polynomial of (31) is given by
(39)Pλ12=μC+λ12J33−λ12J11−λ12λ124+a1λ123+a2λ122+a3λ12+a4,where
(40)a1≡μC∗+γ−J22+J44,(41)a2≡J22J44+γμC∗−γ+μC∗J22+J44,(42)a3≡μC∗+γJ22J44−γμC∗J22+J44,(43)a4≡γJ22J44μC∗−σβBH∗δB∗δE∗=γμC∗δB∗+μB∗J44R0−1.Note that in (43), a4=0 when J44≠0 and R0=1; consequently, one of the roots of (39) is zero. Therefore, in this case the equilibrium point E12 has a zero eigenvalue, that is, it is no-hyperbolic.Three eigenvalues of characteristic polynomial (39) are λ1,12=−μC<0, λ1,22=J33<0 and λ1,32=J11<0. To determine the sign of the other four, which are the roots of the quadratic equation in (39), we will use the Ruth-Hurwitz criterion. According to this, such polynomial has roots with negative real part if and only if all its coefficients a1, a2, a3, and a4 are positive and the relations
(44)a1a2>a3,(45)a1a2a3>a4a12+a32,hold. In order to show these relations, we will adopt the following notation:
(46)A≡γ+μC∗,(47)B≡−J22+J44,(48)C≡γμC∗,(49)D≡J22J44,(50)E≡σβBH∗δB∗δE∗.Note that quantitiesA>0 and C>0, since γ and μC∗ are positive parameters, and E is also positive. By the inequalities (37) and (38), we have that B is positive and D is nonnegative. Thus, in terms of this new notation, the quantities (40)–(43) can be rewritten as
(51)a1=A+B,a2=C+D+AB,a3=AD+BC,a4=CD−E.We note that the coefficientsa1, a2, and a3 are positive, while if J44<0 (or equivalently rE∗/δE∗<rE/δE) and R0<1 also, a4 is positive.Thus, the condition of Ruth-Hurwitz (44) can be written as
(52)A+BC+D+AB>AD+BC.First of all notice that from (46), A2=γ2+2C+μC∗2>2C>C and B2=J222+2D+J442>2D>D. In this way, taking into account that A>0 and B>0, from the above inequalities, we get BA2>BC and AB2>AD. This result allows us to establish that
(53)AB2+BA2>AD+BC.Expanding the left hand side of (52), this can be rewritten as
(54)AB2+A2B+AD+BD+AC+BC>2AD+2BC+BD+AC.In (54), AC>0 since A and C are positive definite, BC>0 since that B>0, and AD≥0 and BD≥0 since that D≥0. We would have on the right hand side of (54) that 2AD+2BC+BD+AC>2AD+2BC>AD+BC. Therefore, (52) is satisfied; that is, the inequality (44) holds.On the other hand, the last condition of Ruth-Hurwitz (45) can be written as
(55)A+BC+D+ABAD+BC>CD−EA+B2+AD+BC2.The left side of (55) can be rewritten as
(56)A+BC+D+ABAD+BC=AB2+A2B+C+DA+BAD+BC.Note that the first factor on the right hand side of (56), taking (53) into account, takes the form
(57)AB2+A2B+C+DA+B>AD+BC+C+DA+B.Thus, the right hand side of (56) becomes
(58)AB2+A2B+C+DA+BAD+BC>AD+BC2+C+DA+BAD+BC.Besides, the second term on the right hand side of (58) can be written as
(59)C+DA+BAD+BC=A2CD+A2D2+ABC2+2ABCD+ABD2+B2C2+B2CD=A+B2CD+BC2+AD2A+B>A+B2CD>A+B2CD−E,since E>0, that is,
(60)C+DA+BAD+BC>A+B2CD−E.Thus, from equations (56), (58), and (60), we have
(61)A+BC+D+ABAD+BC=AB2+A2B+C+DA+BAD+BC>AD+BC2+C+DA+BAD+BC>AD+BC2+A+B2CD−E.In this way, (55) is satisfied; that is, the inequality (45) holds.In summary, by the biological conditionrB>δB+μB, we find that in polynomial (39), its eigenvalues λ1,22=J33, λ1,32=J11, and λ1,12=−μC are negative. Additionally, it has also been shown that the four coefficients a1,a2, a3, and a4 are positive and the two conditions a1a2>a3 and a1a2a3>a4a12+a32 are satisfied when rB∗/δB∗+μB∗≤rB/δB+μB and rE∗/δE∗<rE/δE are fulfilled; in particular, a4 is positive if it also holds that R0<1. Thus, if these three conditions are met: rB∗/δB∗+μB∗≤rB/δB+μB, rE∗/δE∗<rE/δE, and R0<1, then the eigenvalues of the fourth order polynomial in (39) will have a negative real part, and consequently, the point E12 will be asymptotically stable. On the other hand, by Descartes’ rule of signs, if R0>1, then a4<0, and the full polynomial (39) will have a positive eigenvalue, in which case E12 will be unstable. These results can also be summarized in the following theorem.Theorem 11.
Assume that the following conditionsrB>δB+μB,rB∗/δB∗+μB∗≤rB/δB+μB, and rE∗/δE∗<rE/δE are satisfied. The infection-free equilibrium E12 of system (18) is locally asymptotically stable for R0<1 and unstable for R0>1.As already mentioned,E12 is nonhyperbolic when R0=1. It will be shown later that when this happens, a bifurcation occurs.


### 3.2. Global Stability of the Infection-Free Equilibrium

We obtained some conditions on global stability of the infection-free equilibrium of the system (18).Theorem 12.
Assume thatrB>δB+μB, rB=rB∗, rE=rE∗, δB=δB∗, μB=μB∗, and δE∗≥δE. If R0≤1, then the infection-free equilibrium E12=BH∗,0,EH∗,0,CH∗,0,0 of system (18) is globally asymptotically stable in ΩG.Proof.
We construct the following Lyapunov function for the system (18):
(62)Wt=W1t+W2t+W3t,where the following functions have been defined as
(63)W1t=∫BH∗BH1−BH∗ηdη+BI+βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEEI+βσBH∗γμC∗CI+βBH∗γV,W2t=βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KE∫EH∗EH1−EH∗ηdη,W3t=βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KE∫CH∗CH1−CH∗ηdη.
The functionWt is defined, continuous, and positive definite for all BH,BI,EH,EI,CH,CI,V≥0. Also, the global minimum WBH,BI,EH,EI,CH,CI,V=0 occurs at E12=BH∗,0,EH∗,0,CH∗,0,0, and therefore, W is a Lyapunov function. First, we calculate the time derivative of W1t.
(64)dW1tdt=1−BH∗BHdBHdt+dBIdt+βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEdEIdt+βσBH∗γμC∗dCIdt+βBH∗γdVdt=1−BH∗BHrBBH1−BH+BIKB−δB+μBBH−βBHV+rBBI1−BH+BIKB−δB+μBBI+βBHV+βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδBBI+rEEI1−EH+EIKE−δE∗EI+βσBH∗γμC∗δE∗EI−μC∗CI+βBH∗γσCI−γV.
UsingrB−δB+μB=rB/KBBH∗, we have
(65)dW1tdt=−rBKBBH−BH∗BH−BH∗+BI−βBHV+βBH∗V+βBHV+rBBI1−BH∗KB−δB+μBBI−rBKBBIBH−BH∗+BI+βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδBBI+rEEI1−EH∗KE−δEEH−rEKEEIEH−EH∗+EI+βσBH∗γμC∗δE∗EI−μC∗CI+βBH∗γσCI−γV.
After several calculations, we have(66)dW1tdt=−rBKBBH−BH∗+BI2−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEIEH−EH∗+EI+BIβσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KE+rB1−BH∗KB−δB+μB=−rBKBBH−BH∗+BI2−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEIEH−EH∗+EI−δB+μB1−R0BI.
Second, we calculate the time derivative ofW2t.
(67)dW2tdt=1−EH∗EHdEHdt=1−EH∗EHδBBH+rEEH1−EH+EIKE−δEEH.
ConsideringrE=rE/KEEH∗+δE−δBBH∗/EH∗, we obtain
(68)dW2tdt=βσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEBH∗KE1+BH∗BH−EHEH∗−BH∗EHBHEH∗−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEH−EH∗EH−EH∗+EI=βσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEBH∗KE3−BHBH∗−EHEH∗−BH∗EHBHEH∗+βσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEBH∗KEBH∗BH+BHBH∗−2−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEH−EH∗EH−EH∗+EI.
WritingβσδBδE∗BH∗/γμC∗δE∗−rE1−EH∗/KE in terms of R0 and considering that βσδBδE∗BH∗/γμC∗δE∗−rE1−EH∗/KE=δB+μBR0−rB/δB+μB1−BH∗/KB, we find
(69)dW2tdt=βσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEBH∗KE3−BHBH∗−EHEH∗−BH∗EHBHEH∗+R0−rBδB+μB1−BH∗KBδB+μBrEBH∗KEBH∗BH+BHBH∗−2−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEH−EH∗EH−EH∗+EI=βσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEBH∗KE3−BHBH∗−EHEH∗−BH∗EHBHEH∗+1−R0δB+μBrEBH∗KE2−BH∗BH−BHBH∗+1−rBδB+μB1−BH∗KBδB+μBrEBH∗KEBH∗BH+BHBH∗−2−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEH−EH∗EH−EH∗+EI.
Taking into account thatrB1−BH∗/KB=δB+μB, we obtain
(70)dW2tdt=βσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEBH∗KE3−BHBH∗−EHEH∗−BH∗EHBHEH∗+1−R0δB+μBrEBH∗KE2−BH∗BH−BHBH∗−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEH−EH∗EH−EH∗+EI.
Third, we calculate the time derivative ofW3t.
(71)dW3tdt=1−CH∗CHdCHdt=1−CH∗CHδEEH−μCCH.
UsingδEEH∗=μCCH∗, we have
(72)dW3tdt=βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδEEH∗EHEH∗−CHCH∗−EHCH∗EH∗CH+1=βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδEEH∗3−EH∗EH−CHCH∗−EHCH∗EH∗CH+βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδEEH∗EHEH∗+EH∗EH−2,=βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδEEH∗3−EH∗EH−CHCH∗−EHCH∗EH∗CH+R0−rBδB+μB1−BH∗KBδEδB+μBEH∗δBEHEH∗+EH∗EH−2=βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδEEH∗3−EH∗EH−CHCH∗−EHCH∗EH∗CH+1−R0δEδB+μBEH∗δB2−EHEH∗−EH∗EH+rBδB+μB1−BH∗KB−1δEδB+μBEH∗δB2−EHEH∗−EH∗EH.
ConsideringrB1−BH∗/KB=δB+μB, we find
(73)dW3tdt=βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδEEH∗3−EH∗EH−CHCH∗−EHCH∗EH∗CH+1−R0δEδB+μBEH∗δB2−EHEH∗−EH∗EH.
Finally, considering expressions (66), (70), and (73), we obtain
(74)dWtdt=dW1tdt+dW2tdt+dW3tdt=−rBKBBH−BH∗+BI2−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEH−EH∗+EI2−δB+μB1−R0BI+βσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEBH∗KE3−BHBH∗−EHEH∗−BH∗EHBHEH∗+1−R0δB+μBrEBH∗KE2−BH∗BH−BHBH∗+βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδEEH∗3−EH∗EH−CHCH∗−EHCH∗EH∗CH+1−R0δEδB+μBEH∗δB2−EHEH∗−EH∗EH.
IfR0≤1, then dW/dt≤0. Note that dW/dt=0 if and only if BH=BH∗, BI=0, EH=EH∗, and EI=0 or R0=1, BH=BH∗, BI=0, EH=EH∗, and EI=0. Therefore, the largest compact invariant set in BHt,BIt,EHt,EIt,CHt,CIt,Vt:dW/dt=0 is the singleton E12. By the classical LaSalle invariance principle (Theorem 5.3 of [30]), E12 is globally asymptotically stable in ΩG if R0≤1.


### 3.3. Existence of a Forward or Backward Bifurcation

It has already been mentioned that whenR0=1, in the characteristic equation (39), a4 vanishes given rise to a zero eigenvalue and the infection-free equilibrium point E12 is nonhyperbolic. This suggests that in this case, a bifurcation occurs at that point. Indeed, this happens, as we will show below. For this purpose, we will use the Theorem 4.1 of [31]. According to this result, it is required to rewrite system (18) as
(75)dx1dt=a0x1−b0x1x2−βx1x7−b0x12≡f1,dx2dt=a1x2−b1x1x2+βx1x7−b1x22≡f2,dx3dt=δBx1+c0x3−d0x3x4−d0x32≡f3,dx4dt=x2δB∗+c1x4−d1x3x4−d1x42≡f4,dx5dt=δEx3−μCx5≡f5,dx6dt=δE∗x4−μC∗x6≡f6,dx7dt=σx6−γx7≡f7,which is in terms of the new variables
(76)x1≡BH,x2≡BI,x3≡EH,x4≡EI,x5≡CH,x6≡CI,x7≡V,and the new parameters
(77)a0≡rB−δB−μB,a1≡rB∗−δB∗−μB∗,b0≡rBKB,b1≡rB∗KB,c0≡rE−δE,c1≡rE∗−δE∗,d0≡rEKE,d1≡rE∗KE.We will consider thatβ is the bifurcation parameter and that when R0=1, according to definition of R0, this parameter takes the value
(78)β∗=γμC∗J22J44σδB∗δE∗x1∗.The Jacobian matrix of system (75) evaluated at the equilibrium point E12, when β=β∗, is
(79)JE12=J11−b0x1∗0000−β∗x1∗0J220000β∗x1∗δB0J33−d0x3∗0000δB∗0J4400000δE0−μC00000δE∗0−μC∗000000σ−γ,where the quantities J11, J22, J33, and J44, given by (33)–(36), are written as
(80)J11=a0−2b0x1∗,J22=a1−b1x1∗,J33=c0−2d0x3∗,J44=c1−d1x3∗.In this case, we know that zero is a simple eigenvalue of (79). A right eigenvector associated with zero eigenvalue is
(81)w=1μC∗+γJ44J22−γμC∗J22+J44γμC∗J44σδB∗δE∗J11J22−b0x1∗−J44δB∗μC∗δE∗γσγμC∗σδB∗δE∗−δBJ44J11J33J22−b0x1∗+δB∗d0x3∗J33μC∗δE∗γσδEμCγμC∗σδB∗δE∗−δBJ44J11J33J22−b0x1∗+δB∗d0x3∗J33γσ1,and the left eigenvectorv, satisfying v·w=1, is
(82)vT=σδB∗δE∗010−J22δB∗0β∗σx1∗γμC∗β∗x1∗γ.By algebraic calculations, it can be shown that the required nonzero second-order partial derivatives evaluated atE12, which are contained in the quantities a and b given in Theorem 4.1 of [31], are
(83)∂2f2E12∂x1∂x2=∂2f2E12∂x2∂x1=−b1,∂2f2E12∂x1∂x7=∂2f2E12∂x7∂x1=β∗,∂2f2E12∂x2∂x2=−2b1,∂2f4E12∂x3∂x4=∂2f4E12∂x4∂x3=−d1,∂2f4E12∂x4∂x4=−2d1,∂2f2E12∂β∂x7=∂2f2E12∂x7∂β=x1∗.Thus, the quantitiesa and b take the form
(84)a=2β∗v2w1w7−2b1v2w2w1+w2−2d1v4w4w3+w4,(85)b=v2w7x1∗,where w1, w2, w3, w4, and w7 are the first, second, third, fourth, and seventh components of eigenvector (81), while v2 is the second component of left eigenvector (82). Given that in (84) the first term on the right hand side of equality depends on β∗, the relationship that it has with the other terms will determine whether a>0 or a<0. On the other hand, b>0 since v2>0,w7>0, and x1∗>0. Therefore, taking into account the above and Theorem 4.1 of [31], the following theorem can be established.Theorem 13.
The infection-free equilibrium pointE12, when R0=1, presents a backward bifurcation if β∗<βR, while it has a forward bifurcation if β∗>βR, where
(86)βR≡b1v2w2w1+w2v2w1w7+d1v4w4w3+w4v2w1w7.It is important to note that, as a consequence of the previous result, whenR0>1, there is a family of asymptotically stable infected equilibrium points, which we will denote as E22=B¯H,B¯I,E¯H,E¯I,C¯H,C¯I,V¯, constituting the upper branch of these types of bifurcations.


## 3.1. Positivity, Boundedness, Equilibria, Basic Reproductive Number, and Local Stability

Let setΩG be ΩG=BH,BI,EH,EI,CH,CI,V∈ℝ+7:BH≥0,BI≥0,EH≥0,EI≥0,CH≥0,CI≥0,V≥0⊂ℝ7. We consider the positivity and boundedness of the solutions of system (18).Theorem 6.
Given the initial conditionsBH0>0, BI0≥0, EH0>0, EI0≥0, CH0>0, CH0≥0 and V0>0, then all solutions of system (18) are positive.The proof of Theorem1 is very similar to the proof of Theorem 6. For this reason, we omit the proof.Theorem 7.
LetBHt,BIt,EHt,EIt,CHt,CIt,Vt be the solution of model (18) with the initial conditions BH0>0, BI0≥0, EH0>0, EI0≥0, CH0>0, CH0≥0 and V0>0. Then, BHt, BIt, EHt, EIt, CHt, CIt, and Vt are all bounded for all t≥0 at which the solution exists.Proof.
LetBHt,BIt,EHt,EIt,CHt,CIt,Vt be any solution with nonnegative initial conditions. We define a function
(19)Bt=BHt+BIt+EHt+EIt+nn+1CHt+nn+1CIt+μC∗n−1σn+1Vt,n≫1.
The time derivative along a solution of (19) is
(20)dBtdt=rBBH1−BHKB+rB∗BI1−BIKB+rEEH1−EHKE+rE∗EI1−EIKE−rB+rB∗KBBHBI−rE+rE∗KEEHEI−μBBH−μB∗BI−δEn+1EH−δE∗n+1EI−μCnn+1CH−μC∗nnn+1CI−γμC∗n−1σn+1Vt=rB+rB∗KB4+rE+rE∗KE4−rBKBBH−KB22−rB∗KBBI−KB22−rEKEEH−KE22−rE∗KEEI−KE22−rB+rB∗KBBHBI−rE+rE∗KEEHEI−μBBH−μB∗BI−δEn+1EH−δE∗n+1EI−μCnn+1CH−μC∗nnn+1CI−γμC∗n−1σn+1Vt≤rB+rB∗KB+rE+rE∗KE4−μBBH−μB∗BI−δEn+1EH−δE∗n+1EI−μCnn+1CH−μC∗nnn+1CI−γμC∗n−1σn+1Vt.
It follows that(21)dBtdt+ηBt≤rB+rB∗KB+rE+rE∗KE4,where η=minμB,μB∗,δE/n+1,δE∗/n+1,μC,μC∗/n,γ. Thus, limsupt⟶∞Bt≤rB+rB∗KB+rE+rE∗KE/4η. Therefore, BHt, BIt, EHt, EIt, CHt, CIt, and Vt are all bounded for all t≥0.Remark 8.
The total number of cervical epithelial cells, both uninfected and infected, and viral load is bounded by a weighted sum of the stratum basale and stratum intermedium carrying capacities, where the weights are the proliferation rates of healthy and infected cells divided by four times of minimum of the death, differentiation, and viral clearance rates.The model (18) always have a trivial equilibrium E02=0,0,0,0,0,0,0, and a “epithelial cell homeostasis” equilibrium E12=BH∗,0,EH∗,0,CH∗,0,0 with the same coordinates BH∗, EH∗, and CH∗ given in (5), (6) and (7), respectively. We recall the “epithelial cell homeostasis” equilibrium as the infection-free equilibrium.The Jacobian matrix of system (18) is
(22)J=J11−rBBHKB0000−βBHβV−rB∗BIKBJ220000βBHδB0J33−rEEHKE0000δB∗−rE∗EIKEJ4400000δE0−μC00000δE∗0−μC∗000000σ−γ,where
(23)J11=rB−2rBKBBH−rBKBBI−δB−μB−βV,J22=rB∗−rB∗KBBH−2rB∗KBBI−δB∗−μB∗,J33=rE−2rEKEEH−rEKEEI−δE,J44=rE∗−rE∗KEEH−2rE∗KEEI−δE∗.


### 3.1.1. Basic Reproductive NumberR0

To compute the basic reproductive numberR0, we use the next-generation operator introduced by van den Driessche and Watmough [29]. The Jacobian matrix J of this subsystem at the infection-free equilibrium is decomposed as J=F−V, where F is the viral transmission part and V describe the transition terms associated with the model (18). These quantities are given, respectively, by
(24)F≔rB∗1−BH∗KB00βBH∗000000000000,V≔δB∗+μB∗000−δB∗δE∗−rE∗1−EH∗KE000−δE∗μC∗000−σγ.It follows that the basic reproductive numberR0=ρFV−1, where ρ is the spectral radius, is given by
(25)R0=1δB∗+μB∗βσδB∗δE∗BH∗γμC∗δE∗−rE∗1−EH∗/KE+rB∗1−BH∗KB.The parameterR0 has an interesting biological meaning: it is the sum of average numbers of secondary infected cells produced by a single infected cell in a population of epithelial basal cells, by direct basal cell-to-HPV contact and infected basal cell division, respectively.Remark 9.
IfrE∗1−EH∗/KE is taken as a new infection, we can get another basic reproductive number as follows:
(26)R¯0=12a+b+c+12a2+b−c2+2ab+c,where
(27)a=βσδB∗BH∗γμC∗δB∗+μB∗,b=rB∗δB∗+μB∗1−BH∗KB,c=rE∗δE∗1−EH∗KE.
WhenrE∗=0, it is easy to check that R0 is equivalent to R¯0.


### 3.1.2. Local Stability ofE02

The Jacobian matrix (22) evaluated at the equilibrium point E02 becomes
(28)JE02=rB−δB−μB0000000rB∗−δB∗−μB∗00000δB0rE−δE00000δB∗0rE∗−δE∗00000δB0−μC00000δB∗0−μC∗000000σ−γ.By the biological condition, one of its eigenvalues is positive:(29)λ0,12=rB−δB−μB>0,while the other six are
(30)λ0,22=rE−δE,λ0,32=rB∗−δB∗−μB∗,λ0,42=rE∗−δE∗,λ0,52=−μC,λ0,62=−μC∗,λ0,72=−γ.Thus,E02 is unstable. This result can be summarized as follows.Theorem 10.
Assume that the biological conditionrB>δB+μB is satisfied. The trivial equilibrium E02=0,0,0,0,0,0,0 always exists and is unstable.


### 3.1.3. Local Stability ofE12

The Jacobian matrix of system (18), evaluated in the infection-free equilibrium point E12, takes the form
(31)JE12=J11−rBKBBH∗0000−βBH∗0J220000βBH∗δB0J33−rEKEEH∗0000δB∗0J4400000δE0−μC00000δE∗0−μC∗000000σ−γwhere
(32)J11=rB−δB−μB−2rBKBBH∗,J22=rB∗−δB∗−μB∗−rB∗KBBH∗,J33=rE−δE−2rEKEEH∗,J44=rE∗−δE∗−rE∗KEEH∗.Using the following identitiesrB−rBBH∗/KB−δB+μB=0 and rE−rEEH∗/KE−δE=−δBBH∗/EH∗, we have
(33)J11=−rB−δB−μB,(34)J22=rB∗−δB∗−μB∗−rB∗rBrB−δB−μB,(35)J33=−δBBH∗EH∗−rEKEEH∗,(36)J44=rE∗−δE∗−rE∗rErE−δE−rE∗rEδBBH∗EH∗.By biological condition,rB>δB+μB, then J11<0 and J33<0. Furthermore, if
(37)rB∗δB∗+μB∗≤rBδB+μB,(38)rE∗δE∗<rEδE,thenJ22≤0 and J44<0, respectively.On the other hand, the characteristic polynomial of (31) is given by
(39)Pλ12=μC+λ12J33−λ12J11−λ12λ124+a1λ123+a2λ122+a3λ12+a4,where
(40)a1≡μC∗+γ−J22+J44,(41)a2≡J22J44+γμC∗−γ+μC∗J22+J44,(42)a3≡μC∗+γJ22J44−γμC∗J22+J44,(43)a4≡γJ22J44μC∗−σβBH∗δB∗δE∗=γμC∗δB∗+μB∗J44R0−1.Note that in (43), a4=0 when J44≠0 and R0=1; consequently, one of the roots of (39) is zero. Therefore, in this case the equilibrium point E12 has a zero eigenvalue, that is, it is no-hyperbolic.Three eigenvalues of characteristic polynomial (39) are λ1,12=−μC<0, λ1,22=J33<0 and λ1,32=J11<0. To determine the sign of the other four, which are the roots of the quadratic equation in (39), we will use the Ruth-Hurwitz criterion. According to this, such polynomial has roots with negative real part if and only if all its coefficients a1, a2, a3, and a4 are positive and the relations
(44)a1a2>a3,(45)a1a2a3>a4a12+a32,hold. In order to show these relations, we will adopt the following notation:
(46)A≡γ+μC∗,(47)B≡−J22+J44,(48)C≡γμC∗,(49)D≡J22J44,(50)E≡σβBH∗δB∗δE∗.Note that quantitiesA>0 and C>0, since γ and μC∗ are positive parameters, and E is also positive. By the inequalities (37) and (38), we have that B is positive and D is nonnegative. Thus, in terms of this new notation, the quantities (40)–(43) can be rewritten as
(51)a1=A+B,a2=C+D+AB,a3=AD+BC,a4=CD−E.We note that the coefficientsa1, a2, and a3 are positive, while if J44<0 (or equivalently rE∗/δE∗<rE/δE) and R0<1 also, a4 is positive.Thus, the condition of Ruth-Hurwitz (44) can be written as
(52)A+BC+D+AB>AD+BC.First of all notice that from (46), A2=γ2+2C+μC∗2>2C>C and B2=J222+2D+J442>2D>D. In this way, taking into account that A>0 and B>0, from the above inequalities, we get BA2>BC and AB2>AD. This result allows us to establish that
(53)AB2+BA2>AD+BC.Expanding the left hand side of (52), this can be rewritten as
(54)AB2+A2B+AD+BD+AC+BC>2AD+2BC+BD+AC.In (54), AC>0 since A and C are positive definite, BC>0 since that B>0, and AD≥0 and BD≥0 since that D≥0. We would have on the right hand side of (54) that 2AD+2BC+BD+AC>2AD+2BC>AD+BC. Therefore, (52) is satisfied; that is, the inequality (44) holds.On the other hand, the last condition of Ruth-Hurwitz (45) can be written as
(55)A+BC+D+ABAD+BC>CD−EA+B2+AD+BC2.The left side of (55) can be rewritten as
(56)A+BC+D+ABAD+BC=AB2+A2B+C+DA+BAD+BC.Note that the first factor on the right hand side of (56), taking (53) into account, takes the form
(57)AB2+A2B+C+DA+B>AD+BC+C+DA+B.Thus, the right hand side of (56) becomes
(58)AB2+A2B+C+DA+BAD+BC>AD+BC2+C+DA+BAD+BC.Besides, the second term on the right hand side of (58) can be written as
(59)C+DA+BAD+BC=A2CD+A2D2+ABC2+2ABCD+ABD2+B2C2+B2CD=A+B2CD+BC2+AD2A+B>A+B2CD>A+B2CD−E,since E>0, that is,
(60)C+DA+BAD+BC>A+B2CD−E.Thus, from equations (56), (58), and (60), we have
(61)A+BC+D+ABAD+BC=AB2+A2B+C+DA+BAD+BC>AD+BC2+C+DA+BAD+BC>AD+BC2+A+B2CD−E.In this way, (55) is satisfied; that is, the inequality (45) holds.In summary, by the biological conditionrB>δB+μB, we find that in polynomial (39), its eigenvalues λ1,22=J33, λ1,32=J11, and λ1,12=−μC are negative. Additionally, it has also been shown that the four coefficients a1,a2, a3, and a4 are positive and the two conditions a1a2>a3 and a1a2a3>a4a12+a32 are satisfied when rB∗/δB∗+μB∗≤rB/δB+μB and rE∗/δE∗<rE/δE are fulfilled; in particular, a4 is positive if it also holds that R0<1. Thus, if these three conditions are met: rB∗/δB∗+μB∗≤rB/δB+μB, rE∗/δE∗<rE/δE, and R0<1, then the eigenvalues of the fourth order polynomial in (39) will have a negative real part, and consequently, the point E12 will be asymptotically stable. On the other hand, by Descartes’ rule of signs, if R0>1, then a4<0, and the full polynomial (39) will have a positive eigenvalue, in which case E12 will be unstable. These results can also be summarized in the following theorem.Theorem 11.
Assume that the following conditionsrB>δB+μB,rB∗/δB∗+μB∗≤rB/δB+μB, and rE∗/δE∗<rE/δE are satisfied. The infection-free equilibrium E12 of system (18) is locally asymptotically stable for R0<1 and unstable for R0>1.As already mentioned,E12 is nonhyperbolic when R0=1. It will be shown later that when this happens, a bifurcation occurs.


## 3.1.1. Basic Reproductive NumberR0

To compute the basic reproductive numberR0, we use the next-generation operator introduced by van den Driessche and Watmough [29]. The Jacobian matrix J of this subsystem at the infection-free equilibrium is decomposed as J=F−V, where F is the viral transmission part and V describe the transition terms associated with the model (18). These quantities are given, respectively, by
(24)F≔rB∗1−BH∗KB00βBH∗000000000000,V≔δB∗+μB∗000−δB∗δE∗−rE∗1−EH∗KE000−δE∗μC∗000−σγ.It follows that the basic reproductive numberR0=ρFV−1, where ρ is the spectral radius, is given by
(25)R0=1δB∗+μB∗βσδB∗δE∗BH∗γμC∗δE∗−rE∗1−EH∗/KE+rB∗1−BH∗KB.The parameterR0 has an interesting biological meaning: it is the sum of average numbers of secondary infected cells produced by a single infected cell in a population of epithelial basal cells, by direct basal cell-to-HPV contact and infected basal cell division, respectively.Remark 9.
IfrE∗1−EH∗/KE is taken as a new infection, we can get another basic reproductive number as follows:
(26)R¯0=12a+b+c+12a2+b−c2+2ab+c,where
(27)a=βσδB∗BH∗γμC∗δB∗+μB∗,b=rB∗δB∗+μB∗1−BH∗KB,c=rE∗δE∗1−EH∗KE.
WhenrE∗=0, it is easy to check that R0 is equivalent to R¯0.


## 3.1.2. Local Stability ofE02

The Jacobian matrix (22) evaluated at the equilibrium point E02 becomes
(28)JE02=rB−δB−μB0000000rB∗−δB∗−μB∗00000δB0rE−δE00000δB∗0rE∗−δE∗00000δB0−μC00000δB∗0−μC∗000000σ−γ.By the biological condition, one of its eigenvalues is positive:(29)λ0,12=rB−δB−μB>0,while the other six are
(30)λ0,22=rE−δE,λ0,32=rB∗−δB∗−μB∗,λ0,42=rE∗−δE∗,λ0,52=−μC,λ0,62=−μC∗,λ0,72=−γ.Thus,E02 is unstable. This result can be summarized as follows.Theorem 10.
Assume that the biological conditionrB>δB+μB is satisfied. The trivial equilibrium E02=0,0,0,0,0,0,0 always exists and is unstable.


## 3.1.3. Local Stability ofE12

The Jacobian matrix of system (18), evaluated in the infection-free equilibrium point E12, takes the form
(31)JE12=J11−rBKBBH∗0000−βBH∗0J220000βBH∗δB0J33−rEKEEH∗0000δB∗0J4400000δE0−μC00000δE∗0−μC∗000000σ−γwhere
(32)J11=rB−δB−μB−2rBKBBH∗,J22=rB∗−δB∗−μB∗−rB∗KBBH∗,J33=rE−δE−2rEKEEH∗,J44=rE∗−δE∗−rE∗KEEH∗.Using the following identitiesrB−rBBH∗/KB−δB+μB=0 and rE−rEEH∗/KE−δE=−δBBH∗/EH∗, we have
(33)J11=−rB−δB−μB,(34)J22=rB∗−δB∗−μB∗−rB∗rBrB−δB−μB,(35)J33=−δBBH∗EH∗−rEKEEH∗,(36)J44=rE∗−δE∗−rE∗rErE−δE−rE∗rEδBBH∗EH∗.By biological condition,rB>δB+μB, then J11<0 and J33<0. Furthermore, if
(37)rB∗δB∗+μB∗≤rBδB+μB,(38)rE∗δE∗<rEδE,thenJ22≤0 and J44<0, respectively.On the other hand, the characteristic polynomial of (31) is given by
(39)Pλ12=μC+λ12J33−λ12J11−λ12λ124+a1λ123+a2λ122+a3λ12+a4,where
(40)a1≡μC∗+γ−J22+J44,(41)a2≡J22J44+γμC∗−γ+μC∗J22+J44,(42)a3≡μC∗+γJ22J44−γμC∗J22+J44,(43)a4≡γJ22J44μC∗−σβBH∗δB∗δE∗=γμC∗δB∗+μB∗J44R0−1.Note that in (43), a4=0 when J44≠0 and R0=1; consequently, one of the roots of (39) is zero. Therefore, in this case the equilibrium point E12 has a zero eigenvalue, that is, it is no-hyperbolic.Three eigenvalues of characteristic polynomial (39) are λ1,12=−μC<0, λ1,22=J33<0 and λ1,32=J11<0. To determine the sign of the other four, which are the roots of the quadratic equation in (39), we will use the Ruth-Hurwitz criterion. According to this, such polynomial has roots with negative real part if and only if all its coefficients a1, a2, a3, and a4 are positive and the relations
(44)a1a2>a3,(45)a1a2a3>a4a12+a32,hold. In order to show these relations, we will adopt the following notation:
(46)A≡γ+μC∗,(47)B≡−J22+J44,(48)C≡γμC∗,(49)D≡J22J44,(50)E≡σβBH∗δB∗δE∗.Note that quantitiesA>0 and C>0, since γ and μC∗ are positive parameters, and E is also positive. By the inequalities (37) and (38), we have that B is positive and D is nonnegative. Thus, in terms of this new notation, the quantities (40)–(43) can be rewritten as
(51)a1=A+B,a2=C+D+AB,a3=AD+BC,a4=CD−E.We note that the coefficientsa1, a2, and a3 are positive, while if J44<0 (or equivalently rE∗/δE∗<rE/δE) and R0<1 also, a4 is positive.Thus, the condition of Ruth-Hurwitz (44) can be written as
(52)A+BC+D+AB>AD+BC.First of all notice that from (46), A2=γ2+2C+μC∗2>2C>C and B2=J222+2D+J442>2D>D. In this way, taking into account that A>0 and B>0, from the above inequalities, we get BA2>BC and AB2>AD. This result allows us to establish that
(53)AB2+BA2>AD+BC.Expanding the left hand side of (52), this can be rewritten as
(54)AB2+A2B+AD+BD+AC+BC>2AD+2BC+BD+AC.In (54), AC>0 since A and C are positive definite, BC>0 since that B>0, and AD≥0 and BD≥0 since that D≥0. We would have on the right hand side of (54) that 2AD+2BC+BD+AC>2AD+2BC>AD+BC. Therefore, (52) is satisfied; that is, the inequality (44) holds.On the other hand, the last condition of Ruth-Hurwitz (45) can be written as
(55)A+BC+D+ABAD+BC>CD−EA+B2+AD+BC2.The left side of (55) can be rewritten as
(56)A+BC+D+ABAD+BC=AB2+A2B+C+DA+BAD+BC.Note that the first factor on the right hand side of (56), taking (53) into account, takes the form
(57)AB2+A2B+C+DA+B>AD+BC+C+DA+B.Thus, the right hand side of (56) becomes
(58)AB2+A2B+C+DA+BAD+BC>AD+BC2+C+DA+BAD+BC.Besides, the second term on the right hand side of (58) can be written as
(59)C+DA+BAD+BC=A2CD+A2D2+ABC2+2ABCD+ABD2+B2C2+B2CD=A+B2CD+BC2+AD2A+B>A+B2CD>A+B2CD−E,since E>0, that is,
(60)C+DA+BAD+BC>A+B2CD−E.Thus, from equations (56), (58), and (60), we have
(61)A+BC+D+ABAD+BC=AB2+A2B+C+DA+BAD+BC>AD+BC2+C+DA+BAD+BC>AD+BC2+A+B2CD−E.In this way, (55) is satisfied; that is, the inequality (45) holds.In summary, by the biological conditionrB>δB+μB, we find that in polynomial (39), its eigenvalues λ1,22=J33, λ1,32=J11, and λ1,12=−μC are negative. Additionally, it has also been shown that the four coefficients a1,a2, a3, and a4 are positive and the two conditions a1a2>a3 and a1a2a3>a4a12+a32 are satisfied when rB∗/δB∗+μB∗≤rB/δB+μB and rE∗/δE∗<rE/δE are fulfilled; in particular, a4 is positive if it also holds that R0<1. Thus, if these three conditions are met: rB∗/δB∗+μB∗≤rB/δB+μB, rE∗/δE∗<rE/δE, and R0<1, then the eigenvalues of the fourth order polynomial in (39) will have a negative real part, and consequently, the point E12 will be asymptotically stable. On the other hand, by Descartes’ rule of signs, if R0>1, then a4<0, and the full polynomial (39) will have a positive eigenvalue, in which case E12 will be unstable. These results can also be summarized in the following theorem.Theorem 11.
Assume that the following conditionsrB>δB+μB,rB∗/δB∗+μB∗≤rB/δB+μB, and rE∗/δE∗<rE/δE are satisfied. The infection-free equilibrium E12 of system (18) is locally asymptotically stable for R0<1 and unstable for R0>1.As already mentioned,E12 is nonhyperbolic when R0=1. It will be shown later that when this happens, a bifurcation occurs.


## 3.2. Global Stability of the Infection-Free Equilibrium

We obtained some conditions on global stability of the infection-free equilibrium of the system (18).Theorem 12.
Assume thatrB>δB+μB, rB=rB∗, rE=rE∗, δB=δB∗, μB=μB∗, and δE∗≥δE. If R0≤1, then the infection-free equilibrium E12=BH∗,0,EH∗,0,CH∗,0,0 of system (18) is globally asymptotically stable in ΩG.Proof.
We construct the following Lyapunov function for the system (18):
(62)Wt=W1t+W2t+W3t,where the following functions have been defined as
(63)W1t=∫BH∗BH1−BH∗ηdη+BI+βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEEI+βσBH∗γμC∗CI+βBH∗γV,W2t=βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KE∫EH∗EH1−EH∗ηdη,W3t=βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KE∫CH∗CH1−CH∗ηdη.
The functionWt is defined, continuous, and positive definite for all BH,BI,EH,EI,CH,CI,V≥0. Also, the global minimum WBH,BI,EH,EI,CH,CI,V=0 occurs at E12=BH∗,0,EH∗,0,CH∗,0,0, and therefore, W is a Lyapunov function. First, we calculate the time derivative of W1t.
(64)dW1tdt=1−BH∗BHdBHdt+dBIdt+βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEdEIdt+βσBH∗γμC∗dCIdt+βBH∗γdVdt=1−BH∗BHrBBH1−BH+BIKB−δB+μBBH−βBHV+rBBI1−BH+BIKB−δB+μBBI+βBHV+βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδBBI+rEEI1−EH+EIKE−δE∗EI+βσBH∗γμC∗δE∗EI−μC∗CI+βBH∗γσCI−γV.
UsingrB−δB+μB=rB/KBBH∗, we have
(65)dW1tdt=−rBKBBH−BH∗BH−BH∗+BI−βBHV+βBH∗V+βBHV+rBBI1−BH∗KB−δB+μBBI−rBKBBIBH−BH∗+BI+βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδBBI+rEEI1−EH∗KE−δEEH−rEKEEIEH−EH∗+EI+βσBH∗γμC∗δE∗EI−μC∗CI+βBH∗γσCI−γV.
After several calculations, we have(66)dW1tdt=−rBKBBH−BH∗+BI2−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEIEH−EH∗+EI+BIβσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KE+rB1−BH∗KB−δB+μB=−rBKBBH−BH∗+BI2−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEIEH−EH∗+EI−δB+μB1−R0BI.
Second, we calculate the time derivative ofW2t.
(67)dW2tdt=1−EH∗EHdEHdt=1−EH∗EHδBBH+rEEH1−EH+EIKE−δEEH.
ConsideringrE=rE/KEEH∗+δE−δBBH∗/EH∗, we obtain
(68)dW2tdt=βσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEBH∗KE1+BH∗BH−EHEH∗−BH∗EHBHEH∗−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEH−EH∗EH−EH∗+EI=βσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEBH∗KE3−BHBH∗−EHEH∗−BH∗EHBHEH∗+βσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEBH∗KEBH∗BH+BHBH∗−2−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEH−EH∗EH−EH∗+EI.
WritingβσδBδE∗BH∗/γμC∗δE∗−rE1−EH∗/KE in terms of R0 and considering that βσδBδE∗BH∗/γμC∗δE∗−rE1−EH∗/KE=δB+μBR0−rB/δB+μB1−BH∗/KB, we find
(69)dW2tdt=βσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEBH∗KE3−BHBH∗−EHEH∗−BH∗EHBHEH∗+R0−rBδB+μB1−BH∗KBδB+μBrEBH∗KEBH∗BH+BHBH∗−2−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEH−EH∗EH−EH∗+EI=βσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEBH∗KE3−BHBH∗−EHEH∗−BH∗EHBHEH∗+1−R0δB+μBrEBH∗KE2−BH∗BH−BHBH∗+1−rBδB+μB1−BH∗KBδB+μBrEBH∗KEBH∗BH+BHBH∗−2−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEH−EH∗EH−EH∗+EI.
Taking into account thatrB1−BH∗/KB=δB+μB, we obtain
(70)dW2tdt=βσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEBH∗KE3−BHBH∗−EHEH∗−BH∗EHBHEH∗+1−R0δB+μBrEBH∗KE2−BH∗BH−BHBH∗−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEH−EH∗EH−EH∗+EI.
Third, we calculate the time derivative ofW3t.
(71)dW3tdt=1−CH∗CHdCHdt=1−CH∗CHδEEH−μCCH.
UsingδEEH∗=μCCH∗, we have
(72)dW3tdt=βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδEEH∗EHEH∗−CHCH∗−EHCH∗EH∗CH+1=βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδEEH∗3−EH∗EH−CHCH∗−EHCH∗EH∗CH+βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδEEH∗EHEH∗+EH∗EH−2,=βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδEEH∗3−EH∗EH−CHCH∗−EHCH∗EH∗CH+R0−rBδB+μB1−BH∗KBδEδB+μBEH∗δBEHEH∗+EH∗EH−2=βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδEEH∗3−EH∗EH−CHCH∗−EHCH∗EH∗CH+1−R0δEδB+μBEH∗δB2−EHEH∗−EH∗EH+rBδB+μB1−BH∗KB−1δEδB+μBEH∗δB2−EHEH∗−EH∗EH.
ConsideringrB1−BH∗/KB=δB+μB, we find
(73)dW3tdt=βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδEEH∗3−EH∗EH−CHCH∗−EHCH∗EH∗CH+1−R0δEδB+μBEH∗δB2−EHEH∗−EH∗EH.
Finally, considering expressions (66), (70), and (73), we obtain
(74)dWtdt=dW1tdt+dW2tdt+dW3tdt=−rBKBBH−BH∗+BI2−βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEKEEH−EH∗+EI2−δB+μB1−R0BI+βσδBδE∗BH∗γμC∗δE∗−rE1−EH∗/KErEBH∗KE3−BHBH∗−EHEH∗−BH∗EHBHEH∗+1−R0δB+μBrEBH∗KE2−BH∗BH−BHBH∗+βσδE∗BH∗γμC∗δE∗−rE1−EH∗/KEδEEH∗3−EH∗EH−CHCH∗−EHCH∗EH∗CH+1−R0δEδB+μBEH∗δB2−EHEH∗−EH∗EH.
IfR0≤1, then dW/dt≤0. Note that dW/dt=0 if and only if BH=BH∗, BI=0, EH=EH∗, and EI=0 or R0=1, BH=BH∗, BI=0, EH=EH∗, and EI=0. Therefore, the largest compact invariant set in BHt,BIt,EHt,EIt,CHt,CIt,Vt:dW/dt=0 is the singleton E12. By the classical LaSalle invariance principle (Theorem 5.3 of [30]), E12 is globally asymptotically stable in ΩG if R0≤1.


## 3.3. Existence of a Forward or Backward Bifurcation

It has already been mentioned that whenR0=1, in the characteristic equation (39), a4 vanishes given rise to a zero eigenvalue and the infection-free equilibrium point E12 is nonhyperbolic. This suggests that in this case, a bifurcation occurs at that point. Indeed, this happens, as we will show below. For this purpose, we will use the Theorem 4.1 of [31]. According to this result, it is required to rewrite system (18) as
(75)dx1dt=a0x1−b0x1x2−βx1x7−b0x12≡f1,dx2dt=a1x2−b1x1x2+βx1x7−b1x22≡f2,dx3dt=δBx1+c0x3−d0x3x4−d0x32≡f3,dx4dt=x2δB∗+c1x4−d1x3x4−d1x42≡f4,dx5dt=δEx3−μCx5≡f5,dx6dt=δE∗x4−μC∗x6≡f6,dx7dt=σx6−γx7≡f7,which is in terms of the new variables
(76)x1≡BH,x2≡BI,x3≡EH,x4≡EI,x5≡CH,x6≡CI,x7≡V,and the new parameters
(77)a0≡rB−δB−μB,a1≡rB∗−δB∗−μB∗,b0≡rBKB,b1≡rB∗KB,c0≡rE−δE,c1≡rE∗−δE∗,d0≡rEKE,d1≡rE∗KE.We will consider thatβ is the bifurcation parameter and that when R0=1, according to definition of R0, this parameter takes the value
(78)β∗=γμC∗J22J44σδB∗δE∗x1∗.The Jacobian matrix of system (75) evaluated at the equilibrium point E12, when β=β∗, is
(79)JE12=J11−b0x1∗0000−β∗x1∗0J220000β∗x1∗δB0J33−d0x3∗0000δB∗0J4400000δE0−μC00000δE∗0−μC∗000000σ−γ,where the quantities J11, J22, J33, and J44, given by (33)–(36), are written as
(80)J11=a0−2b0x1∗,J22=a1−b1x1∗,J33=c0−2d0x3∗,J44=c1−d1x3∗.In this case, we know that zero is a simple eigenvalue of (79). A right eigenvector associated with zero eigenvalue is
(81)w=1μC∗+γJ44J22−γμC∗J22+J44γμC∗J44σδB∗δE∗J11J22−b0x1∗−J44δB∗μC∗δE∗γσγμC∗σδB∗δE∗−δBJ44J11J33J22−b0x1∗+δB∗d0x3∗J33μC∗δE∗γσδEμCγμC∗σδB∗δE∗−δBJ44J11J33J22−b0x1∗+δB∗d0x3∗J33γσ1,and the left eigenvectorv, satisfying v·w=1, is
(82)vT=σδB∗δE∗010−J22δB∗0β∗σx1∗γμC∗β∗x1∗γ.By algebraic calculations, it can be shown that the required nonzero second-order partial derivatives evaluated atE12, which are contained in the quantities a and b given in Theorem 4.1 of [31], are
(83)∂2f2E12∂x1∂x2=∂2f2E12∂x2∂x1=−b1,∂2f2E12∂x1∂x7=∂2f2E12∂x7∂x1=β∗,∂2f2E12∂x2∂x2=−2b1,∂2f4E12∂x3∂x4=∂2f4E12∂x4∂x3=−d1,∂2f4E12∂x4∂x4=−2d1,∂2f2E12∂β∂x7=∂2f2E12∂x7∂β=x1∗.Thus, the quantitiesa and b take the form
(84)a=2β∗v2w1w7−2b1v2w2w1+w2−2d1v4w4w3+w4,(85)b=v2w7x1∗,where w1, w2, w3, w4, and w7 are the first, second, third, fourth, and seventh components of eigenvector (81), while v2 is the second component of left eigenvector (82). Given that in (84) the first term on the right hand side of equality depends on β∗, the relationship that it has with the other terms will determine whether a>0 or a<0. On the other hand, b>0 since v2>0,w7>0, and x1∗>0. Therefore, taking into account the above and Theorem 4.1 of [31], the following theorem can be established.Theorem 13.
The infection-free equilibrium pointE12, when R0=1, presents a backward bifurcation if β∗<βR, while it has a forward bifurcation if β∗>βR, where
(86)βR≡b1v2w2w1+w2v2w1w7+d1v4w4w3+w4v2w1w7.It is important to note that, as a consequence of the previous result, whenR0>1, there is a family of asymptotically stable infected equilibrium points, which we will denote as E22=B¯H,B¯I,E¯H,E¯I,C¯H,C¯I,V¯, constituting the upper branch of these types of bifurcations.


## 4. Viral KineticIn Silico

In this section, the typical values of all the parameters involved in the HPV viral dynamics model given by the system (18) are indicated. With these values, three relevant features of this model are determined numerically. The first one is the type of local bifurcation that this system exhibits when R0=1 at the infection-free equilibrium point E12. The second is the determination, when R0>1, of the regions of existence and stability of the infected equilibrium point E22. The third consists of the simulation of different scenarios of interest, such as the transient, acute, latent, and chronic infections visualized in Figure 2.


### 4.1. Typical Parameters of the HPV Viral Dynamics Model

Most of the parameter values were obtained from the literature (see Table1). The initial conditions of the number of cells in each stratum were calculated from the cell count in the micrographs of Figure 2 in Walker et al. [32]. We present simulations on a rectangular area of 2585288μm2 of epithelial tissue, where the length is 7282.5μm and the height is 355μm. We estimate KB and KE assuming that the tissue is divided into four equal parts, all having the same area, and each is divided by the area that a cell occupies according to its stratum. Consequently, the initial conditions of the viral kinetic in silico are the following: BH0=1000, BI0=EI0=CI0=0, EH0=617, CH0=226, and V0=100. For parameters rB∗ and rE∗, they were proposed with values close to the proliferation rates of uninfected stratum basale and stratum intermedium cells, respectively.Table 1
Parameters of epithelial cellular dynamics of healthy tissue and dynamics of HPV infection.
ParameterDescriptionValueUnitRangesRefrBUninfected stratum basale cell proliferation rate0.07day−10.03;0.07[25]rB∗Infected stratum basale cell proliferation rate0.048day−1−;−FixedrEUninfected stratum intermedium cell proliferation rate0.039day−10.02;1[25]rE∗Infected stratum intermedium cell proliferation rate0.04day−1−;−FixedδBUninfected cell differentiation rate from stratum basale cell to stratum intermedium0.044day−1−;−[33]δB∗Infected cell differentiation rate from stratum basale cell to stratum intermedium0.05day−1−;−FixedδEUninfected cell differentiation rate from stratum intermedium cell to stratum corneum0.099day−10.02;1[25]δE∗Infected cell differentiation rate from stratum intermedium cell to stratum corneum0.118day−10;5[25]μBUninfected basal cell natural death rate0−;−[34]μB∗Infected basal cell death rate0−;−[34]μCDesquamation rate of uninfected stratum corneum cells0.27day−10.2;1[25]μC∗Desquamation rate of infected stratum corneum cells0.27day−1−;−FixedKBStratum basale cell carrying capacity2693cells1443;13465[32]KEStratum intermedium cell carrying capacity2114cells553;5010[32]βViral transmission rateVariesvirion−1·day−110−15;10−5[25]σProduction rate of free virionsVariescell−1·virion·day−110;103[11]γViral clearance rate1.18day−10.2;3[25]−;− denotes ranges of values not evidenced in the literature.


### 4.2. Existence of a Forward Bifurcation at Infection-Free Equilibrium Point

To determine the type of bifurcation that occurs, whenR0=1, at the equilibrium point E12, were considered β=10−6, σ=103 and all other required quantities from Table 1. With these values it is found, by (85), that b=152,994, which is positive. Furthermore, it is found, according to (78) and (86), respectively, that β∗=9.941×10−8 and βR=4.351×10−8. Since for the typical values of the parameters considered, we find that β∗>βR; consequently, according to Theorem 13, at the infection-free equilibrium point E12, when R0=1, there is a forward bifurcation.


### 4.3. Regions of Existence and Stability of the Infected Equilibrium Point

A numerical scheme was implemented to show the regions of existence and local stability of the infected equilibrium pointE22 of the model (18) when R0 is greater to unity. All of the parameters were set, except the rate of viral transmission β and the production rate of free virions σ. Because an infected cell produces up to 1000 viral particles [11], we take this restriction for the choice of the parameter space σ. For each ordered pair σ,β, we determine numerically the basic reproductive number value and the coordinate V¯ of the infected equilibrium point E22 using the bisection method. We calculate the eigenvalues of the Jacobian matrix (22) evaluated at E22 using the eig function in MATLAB [35], and the signs of the real part of the eigenvalues were identified for the stability classification of the hyperbolic equilibrium point.The regions of existence and stability of the infected equilibrium point are shown in Figure5, given the parameters σ and β. The green, blue, and red regions represent the values of both parameters that make E22 not exist, or be locally asymptotically stable or unstable, respectively.Figure 5
Stability regions of infected equilibrium associated with the parametersσ,β.


### 4.4. Simulation of Different Scenarios of Interest

When some HPV genotype is in a host, there is no infection in the cells of the stratum basale and the viral load is removed until its complete elimination during the following days, without altering the homeostasis of the epithelium. Similar dynamics are consistent with transient type infections.σ=1000 and β=9.447×10−8 such that R0<1 in this scenario is shown in Figure 6. In this case, the infection-free equilibrium point is locally asymptotically stable.Figure 6
Dynamics of the model (18) when there is an infection generated by 100 virions with σ=1000 and β=9.447×10−8, such as R0<1.Some infections successfully establish, replicate the viral genome, produce viral particles, and are eventually cleared. This kinetics is known as acute infection. Figure7 with σ=1000 and β=8.390×10−4 such that R0>1 shows that a σ,β belongs to the instability region of E12 and to the asymptotically stable region of E22, as shown in Figure 5. In this scenario, the homeostasis of the epithelial cells is altered, and the viral load describes a unimodal curve for 1500 days, and finally, the infection is resolved.Figure 7
Dynamics of the model (18) when there is an infection generated by 100 virions with σ=1000 and β=8.390×10−4, such as R0>1.Latent infections are another form of kinetics, where the acute infections only appear to clear, but the viral genome remains in the infected cell without detectable activity. This is illustrated in Figure8 with σ=1000 and β=7.247×10−7 and Figure 9 with σ=1000 and β=2.155×10−6 such that R0>1. The ordered pair σ,β for Figure 8 belongs to the region of asymptotic stability of E22 of Figure 5 and σ,β for Figure 9 belongs to the region of instability of E22 of Figure 5. According to Figure 8, a rapid growth of viral load is shown but without completely infecting the epithelium during the first 250 days after infection. After this period, the viral load tends to decrease up to 1000 days postinfection. From this time on, the dynamics show a similar pattern, which is comparable with damped dynamics until the infection is stabilized throughout the epithelium. This behavior shows the dynamics of a latent infection that can turn into a chronic infection. Likewise, Figure 9 shows the rapid growth of viral load during the first 100 days and decreases until about 2000 days after infection. Although a similar pattern holds over time, the solution does not show damped dynamics as in Figure 8. In this scenario, the resolution of the infection does not occur because the model (18) does not consider mechanisms of prevention or elimination of the virus, the epithelium in its entirety is infected rapidly, and there can be no proliferation of healthy basal cells.Figure 8
Dynamics of the model (18) when there is an infection generated by 100 virions with σ=1000 and β=7.247×10−7, such as R0>1.Figure 9
Dynamics of the model (18) when there is an infection generated by 100 virions with σ=1000 and β=2.155×10−6, such as R0>1.Chronic infections are acute infections that not cleared and maintain viral activity over time. This is shown in Figure10 with σ=1000 and β=1.447×10−7 such that R0>1. The ordered pair σ,β belongs to the region of asymptotic stability of E22 of Figure 5. After infection of the epithelium, viral load production increases monotonically during the first 1500 days until equilibrium is reached, but the infection of the cells does not occur in all stratified epithelial tissue.Figure 10
Dynamics of the model (18) when there is an infection generated by 100 virions with σ=1000 and β=1.447×10−7, such as R0>1.


## 4.1. Typical Parameters of the HPV Viral Dynamics Model

Most of the parameter values were obtained from the literature (see Table1). The initial conditions of the number of cells in each stratum were calculated from the cell count in the micrographs of Figure 2 in Walker et al. [32]. We present simulations on a rectangular area of 2585288μm2 of epithelial tissue, where the length is 7282.5μm and the height is 355μm. We estimate KB and KE assuming that the tissue is divided into four equal parts, all having the same area, and each is divided by the area that a cell occupies according to its stratum. Consequently, the initial conditions of the viral kinetic in silico are the following: BH0=1000, BI0=EI0=CI0=0, EH0=617, CH0=226, and V0=100. For parameters rB∗ and rE∗, they were proposed with values close to the proliferation rates of uninfected stratum basale and stratum intermedium cells, respectively.Table 1
Parameters of epithelial cellular dynamics of healthy tissue and dynamics of HPV infection.
ParameterDescriptionValueUnitRangesRefrBUninfected stratum basale cell proliferation rate0.07day−10.03;0.07[25]rB∗Infected stratum basale cell proliferation rate0.048day−1−;−FixedrEUninfected stratum intermedium cell proliferation rate0.039day−10.02;1[25]rE∗Infected stratum intermedium cell proliferation rate0.04day−1−;−FixedδBUninfected cell differentiation rate from stratum basale cell to stratum intermedium0.044day−1−;−[33]δB∗Infected cell differentiation rate from stratum basale cell to stratum intermedium0.05day−1−;−FixedδEUninfected cell differentiation rate from stratum intermedium cell to stratum corneum0.099day−10.02;1[25]δE∗Infected cell differentiation rate from stratum intermedium cell to stratum corneum0.118day−10;5[25]μBUninfected basal cell natural death rate0−;−[34]μB∗Infected basal cell death rate0−;−[34]μCDesquamation rate of uninfected stratum corneum cells0.27day−10.2;1[25]μC∗Desquamation rate of infected stratum corneum cells0.27day−1−;−FixedKBStratum basale cell carrying capacity2693cells1443;13465[32]KEStratum intermedium cell carrying capacity2114cells553;5010[32]βViral transmission rateVariesvirion−1·day−110−15;10−5[25]σProduction rate of free virionsVariescell−1·virion·day−110;103[11]γViral clearance rate1.18day−10.2;3[25]−;− denotes ranges of values not evidenced in the literature.


## 4.2. Existence of a Forward Bifurcation at Infection-Free Equilibrium Point

To determine the type of bifurcation that occurs, whenR0=1, at the equilibrium point E12, were considered β=10−6, σ=103 and all other required quantities from Table 1. With these values it is found, by (85), that b=152,994, which is positive. Furthermore, it is found, according to (78) and (86), respectively, that β∗=9.941×10−8 and βR=4.351×10−8. Since for the typical values of the parameters considered, we find that β∗>βR; consequently, according to Theorem 13, at the infection-free equilibrium point E12, when R0=1, there is a forward bifurcation.


## 4.3. Regions of Existence and Stability of the Infected Equilibrium Point

A numerical scheme was implemented to show the regions of existence and local stability of the infected equilibrium pointE22 of the model (18) when R0 is greater to unity. All of the parameters were set, except the rate of viral transmission β and the production rate of free virions σ. Because an infected cell produces up to 1000 viral particles [11], we take this restriction for the choice of the parameter space σ. For each ordered pair σ,β, we determine numerically the basic reproductive number value and the coordinate V¯ of the infected equilibrium point E22 using the bisection method. We calculate the eigenvalues of the Jacobian matrix (22) evaluated at E22 using the eig function in MATLAB [35], and the signs of the real part of the eigenvalues were identified for the stability classification of the hyperbolic equilibrium point.The regions of existence and stability of the infected equilibrium point are shown in Figure5, given the parameters σ and β. The green, blue, and red regions represent the values of both parameters that make E22 not exist, or be locally asymptotically stable or unstable, respectively.Figure 5
Stability regions of infected equilibrium associated with the parametersσ,β.


## 4.4. Simulation of Different Scenarios of Interest

When some HPV genotype is in a host, there is no infection in the cells of the stratum basale and the viral load is removed until its complete elimination during the following days, without altering the homeostasis of the epithelium. Similar dynamics are consistent with transient type infections.σ=1000 and β=9.447×10−8 such that R0<1 in this scenario is shown in Figure 6. In this case, the infection-free equilibrium point is locally asymptotically stable.Figure 6
Dynamics of the model (18) when there is an infection generated by 100 virions with σ=1000 and β=9.447×10−8, such as R0<1.Some infections successfully establish, replicate the viral genome, produce viral particles, and are eventually cleared. This kinetics is known as acute infection. Figure7 with σ=1000 and β=8.390×10−4 such that R0>1 shows that a σ,β belongs to the instability region of E12 and to the asymptotically stable region of E22, as shown in Figure 5. In this scenario, the homeostasis of the epithelial cells is altered, and the viral load describes a unimodal curve for 1500 days, and finally, the infection is resolved.Figure 7
Dynamics of the model (18) when there is an infection generated by 100 virions with σ=1000 and β=8.390×10−4, such as R0>1.Latent infections are another form of kinetics, where the acute infections only appear to clear, but the viral genome remains in the infected cell without detectable activity. This is illustrated in Figure8 with σ=1000 and β=7.247×10−7 and Figure 9 with σ=1000 and β=2.155×10−6 such that R0>1. The ordered pair σ,β for Figure 8 belongs to the region of asymptotic stability of E22 of Figure 5 and σ,β for Figure 9 belongs to the region of instability of E22 of Figure 5. According to Figure 8, a rapid growth of viral load is shown but without completely infecting the epithelium during the first 250 days after infection. After this period, the viral load tends to decrease up to 1000 days postinfection. From this time on, the dynamics show a similar pattern, which is comparable with damped dynamics until the infection is stabilized throughout the epithelium. This behavior shows the dynamics of a latent infection that can turn into a chronic infection. Likewise, Figure 9 shows the rapid growth of viral load during the first 100 days and decreases until about 2000 days after infection. Although a similar pattern holds over time, the solution does not show damped dynamics as in Figure 8. In this scenario, the resolution of the infection does not occur because the model (18) does not consider mechanisms of prevention or elimination of the virus, the epithelium in its entirety is infected rapidly, and there can be no proliferation of healthy basal cells.Figure 8
Dynamics of the model (18) when there is an infection generated by 100 virions with σ=1000 and β=7.247×10−7, such as R0>1.Figure 9
Dynamics of the model (18) when there is an infection generated by 100 virions with σ=1000 and β=2.155×10−6, such as R0>1.Chronic infections are acute infections that not cleared and maintain viral activity over time. This is shown in Figure10 with σ=1000 and β=1.447×10−7 such that R0>1. The ordered pair σ,β belongs to the region of asymptotic stability of E22 of Figure 5. After infection of the epithelium, viral load production increases monotonically during the first 1500 days until equilibrium is reached, but the infection of the cells does not occur in all stratified epithelial tissue.Figure 10
Dynamics of the model (18) when there is an infection generated by 100 virions with σ=1000 and β=1.447×10−7, such as R0>1.


## 5. Local Sensitivity Analysis

In the context of viral infections, sensitivity analysis can provide valuable information on how HPV viral kinetics are with respect to changes in model parameters. Local sensitivity analysis is a classic method that studies the impact of small perturbations on the model outputs. The normalized forward sensitivity index of a variable,u, that depends differentiable on a parameter, p, is defined as
(87)Ψpu=pu∂u∂p.Local sensitivity analysis is carried out in order to identify which parameters have the greatest influence on changes in the values ofR0. As we have an explicit formula for R0, we derive an analytical expression for the sensitivity of R0,
(88)ΨpR0=pR0∂R0∂p,to each of the parameters described in R0.Figure11 shows that all the parameters have either positive or negative effects on the R0. According to the sensitivity indices illustrated in Figure 11, we observe that the parameters, namely, rE∗, β, σ, and δE∗, are the most positively sensitive parameters, respectively. This implies that decreasing the values of these parameters will decrease R0. Parameters μC∗, δB, δB∗, and γ are the most negatively sensitive parameters. This implies that increasing the values of these parameters will decrease R0.Figure 11
Local sensitivity analysis of the basic reproduction numberR0.


## 6. Discussion

Several approaches have been reported in the literature for the mathematical modeling of the process of viral infection by HPV. For example, Verma et al. [24] proposed a model of HIV/HPV coinfection in oral mucosal epithelial cells with anti-HPV immune response. They consider the basal and suprabasal layers of oral mucosal but ignore viral transmission via suprabasal differentiation, which is very relevant in the persistence of the virus. In their simulations, the authors reproduce acute and chronic infections. Murall et al. [25] developed a viral dynamics model of the cervical stratified epithelium of cornified, granular, spinous, and basal layers. In their simulations, they reproduce a scenario where the infection is inoculated with a few cells and the microabrasion repairs quickly, and another scenario of slow growing HR HPV infection that spontaneously regresses. However, this model ignores the epithelial cellular dynamics of healthy tissue, and it models the viral transmission via suprabasal differentiation with a linear term. Motivated by this review, we proposed a mathematical model that includes several layers of the healthy and infected squamous epithelium without the presence of the immune response. We modeled the viral transmission via suprabasal differentiation with a full logistic term.We start by proposing a model of the epithelial cellular dynamics of the cervix stratified epithelium of three (basale, intermedium, and corneum) stratums, given by the system (1), where the stratum intermedium is formed by granular and spinous layers. This assumption considers that the dynamics of basal and suprabasal differentiation are homogeneous between the two stratums. In addition, we have a simpler mathematical model that is easier to deal with from qualitative analysis. We determine biological condition rB>δB+μB for the existence of the epithelial cell homeostasis equilibrium E11=BH∗,EH∗,CH∗, and consequently, we include cellular dynamics of the cervical epithelium. We have proven that trivial equilibrium E01=0,0,0 always exists and is unstable. Using the method of Lyapunov functions and a theorem on limiting systems, we have proven that the epithelial cell homeostasis equilibrium is globally asymptotically stable when the biological condition is satisfied.Later, we formulated an extended model, given by system (18), in which the infection by the HPV virus was taken into account. For this system, the basic reproductive number R0 of the infection was calculated. This allowed us to design the different scenarios of the theoretical viral loads. The model also has a trivial equilibrium point E02 and an epithelial cell homeostasis one E12. It was proved that E02 always exists and is unstable, whereas if the biological condition rB>δB+μB is satisfied, then E12 exists. If the following conditions rB∗/δB∗+μB∗≤rB/δB+μB and rE∗/δE∗<rE/δE are satisfied, the E12 is locally asymptotically stable if R0<1 and unstable when R0>1. We note that the stability conditions have the following biological interpretation: rE∗/δE∗<rE/δE, the ratio of proliferation rate and differentiation rate of the infected stratum intermedium cell is less than ratio of proliferation rate and differentiation rate of the uninfected stratum intermedium cells. A similar interpretation can be made of the condition rB∗/δB∗+μB∗≤rB/δB+μB. Using an elegant construction of a Lyapunov function, we obtained conditions on parameters (rB=rB∗, rE=rE∗, δB=δB∗, μB=μB∗, and δE∗≥δE) to prove the global asymptotic stability of the epithelial cell homeostasis equilibrium. We observe that under the following conditions on the parameters rB=rB∗, rE=rE∗, δB=δB∗, μB=μB∗, and δE∗>δE, the conditions of the Theorem 11 continue to be satisfied. Furthermore, it was shown that when R0=1, E12 is nonhyperbolic and that in this case, it experiences a forward bifurcation. This last result shows the existence of a family of asymptotically stable infected equilibrium points, denoted as E22, which bifurcate from the nonhyperbolic equilibrium point and are located in the upper branch of the forward bifurcation.We numerically study the local stability of the infected equilibriumE22 varying the virus-to-cell transmission rate and the viral production rate, and we obtained the regions of stability that are shown in Figure 5. We have reproduced the viral kinetics in silico that have been proposed by Alizon et al. [12] as the transient, acute, latent, and chronic infections. In the transient infection (Figure 6), the viral load is removed until its complete elimination during the following days, without altering the homeostasis of the epithelium. In acute infection, the viral load describes a unimodal curve for fifteen hundred days, and the homeostasis of the epithelial cells is altered; finally, the infection is resolved (Figure 7). In latent infection, the viral load shows behavior of damped and self-sustaining oscillations (Figures 8 and 9), and the resolution of the infection does not occur, and the epithelium in its entirety is infected rapidly. In the chronic infection, the viral load production increases monotonically during the first thousand days until equilibrium is reached, but the infection of the cells does not occur in all stratified epithelial tissue (Figure 10). All of the simulations were performed with the initial condition of viral inoculation of V0=100. We also performed simulations with different initial values in V0, specifically with V0=10 and V0=1000, and we obtained the same scenarios for each case.As described by Alizon et al. [12], transient, acute, latent, and chronic infections differ in terms of viral activity, such as viral and cellular gene expression patterns, effects on cell replication dynamics, or induced local immunosuppression. Currently, protocols are being developed, such as the PAPCLEAR study [15], that allow adequate monitoring to characterize the stages of infection in healthy young women, particularly in the detection of viral genetic material associated with latent or chronic infections. The PAPCLEAR study will be relevant because of its impact in understanding the natural history of cervical HPV infections as the possible integration of longitudinal data to mathematical models.Sensitivity analysis provides insights on possible strategies for the control of a viral infection. The results of the local sensitivity analysis (Figure11) should be considered together with simulated outputs of virus dynamics models and the possible interventions to be carried out. In our model of the dynamics of HPV infection, the viral transmission rate (β) and viral production rate (σ) could be reduced by developing antiviral therapies targeting inhibition of new infections and viral replication, respectively, which is still under investigation [36, 37]. The viral clearance rate (γ) is traditionally increased by the antibody immune response induced by vaccines. Vaccine-induced antibody levels have been reported to be stable over time, which is associated with high long-term protection [38]. Differentiation rate of infected basal cell (δB∗) could be decreased by cytotoxic T cell immune response, but there is no direct evidence that viral antigen-specific immune effector mechanisms are responsible for virus elimination [39].Our model (18) suffers from a limitation because it does not consider the four strata of the stratified epithelium, but the model illustrates all the theoretical viral kinetic scenarios proposed by Alizon et al. [12]. Therefore, the innate and cellular immune response can be studied in the future work where it is possible to reproduce the kinetics obtained in the retrospective study [13] as the kinetics of latency with regression or progression of the infection.


## 7. Conclusions

In summary, in this research work, the dynamics of an HPV model were studied through the use of qualitative and numerical analyses. Regarding the first, the positivity and boundary conditions of their solutions were determined, and their main equilibrium points were found, as well as their local and global stabilities. In addition, it was shown that, whenR0=1, the model has a nonhyperbolic equilibrium point which, for typical values of the parameters involved, gives rise to a forward bifurcation; consequently, there is a family of asymptotically stable infected equilibrium points E22 that branch off from the nonhyperbolic point. It is worth mentioning that this last result is merely local and only valid in the neighborhood of R0=1, so the qualitative analysis of the local and global stabilities of these points for values of R0 much greater than one is still an open problem.Through numerical analysis of the model solutions, we study some features of the intricate behavior that they exhibit around the infected equilibrium pointE22. Specifically, our simulated results (i.e., viral kinetic in silico) suggested that the dynamics of HPV model reproduce the transient, acute, latent, and chronic infections that have been reported in studies of the natural history of HPV.Finally, using local sensitivity analysis, we found some parameters that could be controlled to remove HPV infection in epithelial tissue, as the viral transmission rate (β), viral production rate (σ), and viral clearance rate (γ).



---

*Source: 8650071-2022-08-23.xml*