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

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.


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 highrisk (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 Figure 1). 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.
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 Figure 2, 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].
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 HPVinfected 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][17][18][19][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][22][23][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 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 simulate in silico the types of infection reported in the literature when HPV infection occurs.
The structure of this manuscript is as follows. In Section 2, 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.

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 Figure 1). Consequently, the cell dynamic model of the homeostatic 3 Computational and Mathematical Methods in Medicine stratified epithelium is composed of the stratum basale, stratum intermedium, and stratum corneum. The dynamic is visualized in Figure 3.
The model is given by B H ðtÞ, E H ðtÞ, and C H ðtÞ, which denote stratum basale, stratum intermedium, and stratum corneum formed by uninfected cells, respectively. The cells of B H ðtÞ proliferate (or lateral differentiation) at rate r B considering a logistic growth, with a carrying capacity K B ; these have suprabasal differentiation to E H ðtÞ at rate δ B and die at rate μ B . The dynamics of E H ðtÞ is generated by suprabasal differentiation of B H ðtÞ at a rate δ B . This population increases by cell proliferation, which is governed by logistic growth at a rate r E and a carrying capacity at rate K E . These decrease by suprabasal differentiation at a rate δ E , and we assume that there is no natural death. Finally, the dynamics of C H ðtÞ are generated by suprabasal differentiation of E H ðtÞ at rate δ E and shed from the epithelial tissue at a rate μ C . We arrive at the following mathematical model: Proof. At first, we will prove the positivity of B H ðtÞ. Let B H ð0Þ be the solution that satisfies the initial condition B H ðtÞ > 0. Assume that the solution is not always positive; i.e., there exists a t 0 ′ ∈ ℝ + , such as B H ðt 0 ′Þ < 0. By Bolzano's theorem, there exists a t 1 ∈ ð0, t 0 ′ Þ, such as B H ðt 1 Þ = 0. Let t 0 ∈ ℝ + be the initial time, such as B H ðt 0 Þ = 0, and then dB H ðtÞ/dt = 0. Note that if some t ≥ 0, B H ðtÞ = 0, then dB H ðtÞ/dt = 0. Then, any solution with B H ð0Þ = 0 will satisfy B H ðtÞ = 0 ∀t > 0. By uniqueness of solutions, we have that if B H ð0Þ > 0, then B H ðtÞ will remain positive ∀t > 0. Therefore, B H ðt 0 Þ = 0 leads to a contradiction. Hence, B H ðtÞ is nonnegative for all t > 0. Now, we will prove the positivity of E H ðtÞ. E H ðtÞ is not always positive; i.e., there exists t 0 ′ ∈ ℝ + , such as E H ðt 0 ′ Þ < 0. By Bolzano's theorem, there exists t 1 ∈ ð0, t 0 ′ Þ, such as E H ðt 1 Þ = 0. Let t 0 = min ft i | E H ðt i Þ = 0g. By the second equation (1), if E H ðt 0 Þ = 0, then dE H ðt 0 Þ/dt = δ b B H > 0 implies that E H is increasing at t = t 0 . Therefore, E H ðtÞ will be negative for values t < t 0 near to t 0 , that is, a contradiction. Finally, we will prove the positivity of C H . From the third equation (1), we obtain the following inequality dC H ð tÞ/dt ≥ −μ C C H ðtÞ. Integrating, we obtain the solution C H ≥ C H ð0Þe −μt ; therefore, C H ≥ 0. Proof. Let (B H ðtÞ, E H ðtÞ, C H ðtÞ) be any solution with nonnegative initial conditions. We define a function The time derivative along a solution of (2) is It follows that where η = min fμ B , δ E /n, μ C g. 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 equilibrium E 1 0 = ð0, 0, 0Þ always exists. We get the "epithelial cell homeostasis" equilibrium E 1 The following inequality r B > δ B + μ B is a biological condition for the maintenance of cell homeostasis of the epithelium.
The Jacobian matrix of (1) at E 1 0 is The eigenvalues of then JðE 1 0 Þ has one positive eigenvalue. Thus, the trivial equilibrium is unstable.
The Jacobian matrix for the equations of (1) at E 1 1 is Using the following identities The characteristic polynomial of JðE 1 1 Þ is The eigenvalues of Pðλ 1 1 Þ are λ 1 Clearly, all of the eigenvalues are negative. Consequently, the epithelial cell homeostasis equilibrium E 1 1 is locally asymptotically stable.
We then arrive at the following theorem. 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 the B H − E H subsystem of model (1), which is independent of the C H variables. We construct the following Volterra-type Lyapunov function [26] for the B H − E H subsystem The function UðtÞ is defined, continuous, and positive definite for all B H , E H > 0. Also, the global minimum UðB H , E H Þ = 0 occurs at ðB * H , E * H Þ, and therefore, U is a Lyapunov function. First, we calculate the time derivative of UðtÞ computed along solutions of the first two equations of the system (1), given by the expression Note that Using (14) and (15), we have The viral dynamics model includes that of uninfected cells, infected cells, and viral load. The viral dynamics are given by B I ðtÞ, E I ðtÞ, C I ðtÞ, and VðtÞ, which denote the cells of the stratum basale, stratum intermedium, stratum corneum infected, and the viral load, respectively. B H ðtÞ, E H ðtÞ, and C H ðtÞ are denoted in the same way as in the previous section. Figure 4 shows that the dynamics of B I ðtÞ 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 r * B , and a carrying capacity at rate K B . These decrease by suprabasal differentiation at rate δ * B and death cellular at rate μ * B . The dynamics of E I ðtÞ are generated by suprabasal differentiation of B I ðtÞ at a rate δ * B . This population increases by cell proliferation, governed by full logistic growth at a rate r * E , and a carrying capacity at rate K E . These decrease by suprabasal differentiation at a rate δ * E , and we assume that there is no natural death. The dynamics of C I ðtÞ are gener-ated by suprabasal differentiation of E I ðtÞ at rate δ * E . There is no proliferation of C I ðtÞ. This population decreases by desquamation at rate μ * C . Finally, VðtÞ increases by rupture of cell membrane of C I ðtÞ at rate σ, and they decline at rate γ. This is summarized in the following nonlinear ODE system: We consider the positivity and boundedness of the solutions of system (18).  Figure 4: Diagram of viral dynamics when a HPV infection occurs. 6 Computational and Mathematical Methods in Medicine The proof of Theorem 1 is very similar to the proof of Theorem 6. For this reason, we omit the proof.
The time derivative along a solution of (19) is It follows that where 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 Jacobian matrix of system (18) is where 3.1.1. Basic Reproductive Number R 0 . To compute the basic reproductive number R 0 , we use the next-generation operator introduced by van den Driessche and Watmough [29]. The Jacobian matrix J of this subsystem at the infectionfree 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 7 Computational and Mathematical Methods in Medicine given, respectively, by It follows that the basic reproductive number R 0 = ρðF V −1 Þ, where ρ is the spectral radius, is given by The parameter R 0 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.
taken as a new infection, we can get another basic reproductive number as follows: where When r * E = 0, it is easy to check that R 0 is equivalent to R 0 .

Local
Stability of E 2 0 . The Jacobian matrix (22) evaluated at the equilibrium point E 2 0 becomes By the biological condition, one of its eigenvalues is positive: while the other six are Thus, E 2 0 is unstable. This result can be summarized as follows.
The Jacobian matrix of system (18), evaluated in the infection-free equilibrium point E 2 1 , takes the form Using the following identities By biological condition, then J 22 ≤ 0 and J 44 < 0, respectively. On the other hand, the characteristic polynomial of (31) is given by where Note that in (43), a 4 = 0 when J 44 ≠ 0 and R 0 = 1; consequently, one of the roots of (39) is zero. Therefore, in this case the equilibrium point E 2 1 has a zero eigenvalue, that is, it is no-hyperbolic.
Three eigenvalues of characteristic polynomial (39) are λ 2 1,1 = −μ C < 0, λ 2 1,2 = J 33 < 0 and λ 2 1,3 = J 11 < 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 a 1 , a 2 , a 3 , and a 4 are positive and the relations hold. In order to show these relations, we will adopt the following notation: Note that quantities A > 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 We note that the coefficients a 1 , a 2 , and a 3 are positive, while if J 44 < 0 (or equivalently r * E /δ * E < r E /δ E ) and R 0 < 1 also, a 4 is positive.
Thus, the condition of Ruth-Hurwitz (44) can be written as First of all notice that from (46), A 2 = γ 2 + 2C + ðμ * C Þ 2 > 2C > C and B 2 = J 2 22 + 2D + J 2 44 > 2D > D: In this way, taking into account that A > 0 and B > 0, from the above inequalities, we get BA 2 > BC and AB 2 > AD: This result allows us to establish that Expanding the left hand side of (52), this can be 9 Computational and Mathematical Methods in Medicine rewritten as 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 2A D + 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 The left side of (55) can be rewritten as Note that the first factor on the right hand side of (56), taking (53) into account, takes the form Thus, the right hand side of (56) becomes Besides, the second term on the right hand side of (58) can be written as since E > 0, that is, Thus, from equations (56), (58), and (60), we have In this way, (55) is satisfied; that is, the inequality (45) holds.
In summary, by the biological condition r B > δ B + μ B , we find that in polynomial (39), its eigenvalues λ 2 1,2 = J 33 , λ 2 1,3 = J 11 , and λ 2 1,1 = −μ C are negative. Additionally, it has also been shown that the four coefficients a 1 ,a 2 , a 3 , and a 4 are positive and the two conditions a 1 a 2 > a 3 and a 1 a 2 a 3 > a 4 a 2 1 + a 2 3 are satisfied when r * B /ðδ * B + μ * B Þ ≤ r B /ðδ B + μ B Þ and r * E /δ * E < r E /δ E are fulfilled; in particular, a 4 is positive if it also holds that R 0 < 1. Thus, if these three conditions are met: and R 0 < 1, then the eigenvalues of the fourth order polynomial in (39) will have a negative real part, and consequently, the point E 2 1 will be asymptotically stable. On the other hand, by Descartes' rule of signs, if R 0 > 1, then a 4 < 0, and the full polynomial (39) will have a positive eigenvalue, in which case E 2 1 will be unstable. These results can also be summarized in the following theorem.
Theorem 11. Assume that the following conditions (18) is locally asymptotically stable for R 0 < 1 and unstable for R 0 > 1.
As already mentioned, E 2 1 is nonhyperbolic when R 0 = 1. It will be shown later that when this happens, a bifurcation occurs.

Global Stability of the Infection-Free Equilibrium.
We obtained some conditions on global stability of the infection-free equilibrium of the system (18). (18) is globally asymptotically stable in Ω G .
Proof. We construct the following Lyapunov function for the system (18): where the following functions have been defined as Using After several calculations, we have Second, we calculate the time derivative of W 2 ðtÞ. Considering ð68Þ Taking into account that r B ð1 − B * H /K B Þ = δ B + μ B , we obtain Third, we calculate the time derivative of W 3 ðtÞ.
3.3. Existence of a Forward or Backward Bifurcation. It has already been mentioned that when R 0 = 1, in the characteristic equation (39), a 4 vanishes given rise to a zero eigenvalue and the infection-free equilibrium point E 2 1 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 which is in terms of the new variables and the new parameters We will consider that β is the bifurcation parameter and that when R 0 = 1, according to definition of R 0 , this parameter takes the value The Jacobian matrix of system (75) evaluated at the 12 Computational and Mathematical Methods in Medicine where the quantities J 11 , J 22 , J 33 , and J 44 , given by (33)- (36), are written as In this case, we know that zero is a simple eigenvalue of (79). A right eigenvector associated with zero eigenvalue is and the left eigenvector v, By algebraic calculations, it can be shown that the required nonzero second-order partial derivatives evaluated at E 2 1 , which are contained in the quantities a and b given in Theorem 4.1 of [31], are Thus, the quantities a and b take the form 13 Computational and Mathematical Methods in Medicine where w 1 , w 2 , w 3 , w 4 , and w 7 are the first, second, third, fourth, and seventh components of eigenvector (81), while v 2 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 v 2 > 0,w 7 > 0, and x * 1 > 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 point E 2 1 , when
It is important to note that, as a consequence of the previous result, when R 0 > 1, there is a family of asymptotically stable infected equilibrium points, which we will denote as

Viral Kinetic In Silico
In this section, the typical values of all the parameters involved in the HPV viral dynamics model given by the sys-tem (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 R 0 = 1 at the infection-free equilibrium point E 2 1 . The second is the determination, when R 0 > 1, of the regions of existence and stability of the infected equilibrium point E 2 2 . The third consists of the simulation of different scenarios of interest, such as the transient, acute, latent, and chronic infections visualized in Figure 2.

Typical Parameters of the HPV Viral Dynamics Model.
Most of the parameter values were obtained from the literature (see Table 1). 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 μm 2 of epithelial tissue, where the length is 7282:5 μm and the height is 355 μm. We estimate K B and K E 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: B H ð0Þ = 1000, B I ð0Þ = E I ð0Þ = C I ð0Þ = 0, E H ð0Þ = 617, C H ð 0Þ = 226, and Vð0Þ = 100. For parameters r * B and r * E , they were proposed with values close to the proliferation rates of uninfected stratum basale and stratum intermedium cells, respectively.

Existence of a Forward Bifurcation at Infection-Free Equilibrium Point.
To determine the type of bifurcation that occurs, when R 0 = 1, at the equilibrium point E 2 1 , were considered β = 10 −6 , σ = 10 3 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 × 1 0 −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 E 2 1 , when R 0 = 1, there is a forward bifurcation.

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 point E 2 2 of the model (18) when R 0 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 E 2 2 using the bisection method. We calculate the eigenvalues of the Jacobian matrix (22) evaluated at E 2 2 using the eig function in MATLAB [35], and the signs of the real part of the eigen-values 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 Figure 5, given the parameters σ and β. The green, blue, and red regions represent the values of both parameters that make E 2 2 not exist, or be locally asymptotically stable or unstable, respectively.

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 R 0 < 1 in this scenario is shown in Figure 6. In this case, the infection-free equilibrium point is locally asymptotically stable.
Some infections successfully establish, replicate the viral genome, produce viral particles, and are eventually cleared. This kinetics is known as acute infection. Figure 7 with σ = 1000 and β = 8:390 × 10 −4 such that R 0 > 1 shows that a ðσ, βÞ belongs to the instability region of E 2 1 and to the asymptotically stable region of E 2 2 , 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.

Computational and Mathematical Methods in Medicine
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 Figure 8 with σ = 1000 and β = 7:247 × 10 −7 and Figure 9 with σ = 1000 and β = 2:155 × 10 −6 such that R 0 > 1. The ordered pair ðσ, βÞ for Figure 8 belongs to the region of asymptotic stability of E 2 2 of Figure 5 and ðσ, βÞ for Figure 9 belongs to the region of instability of E 2 2 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.
Chronic infections are acute infections that not cleared and maintain viral activity over time. This is shown in Figure 10 with σ = 1000 and β = 1:447 × 10 −7 such that R 0 > 1. The ordered pair ðσ, βÞ belongs to the region of asymptotic stability of E 2 2 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.

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 Local sensitivity analysis is carried out in order to identify which parameters have the greatest influence on changes in the values of R 0 . As we have an explicit formula for R 0 , we derive an analytical expression for the sensitivity of R 0 , to each of the parameters described in R 0 . Figure 11 shows that all the parameters have either positive or negative effects on the R 0 . According to the sensitivity indices illustrated in Figure 11, we observe that the parameters, namely, r * E , β, σ, and δ * E , are the most positively sensitive parameters, respectively. This implies that decreasing the values of these parameters will decrease R 0 . Parameters μ * C , δ B , δ * B , and γ are the most negatively sensitive parameters. This implies that increasing the values of these parameters will decrease R 0 .

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 r B > δ B + μ B for the existence of the epithelial cell homeostasis equilibrium E 1 1 = ðB * H , E * H , C * H Þ, and consequently, we include cellular dynamics of the cervical epithelium. We have proven that trivial equilibrium E 1 0 = ð0, 0, 0Þ always exists and is unstable. Using the method of Lyapunov functions and a theorem on limiting systems, we have proven

18
Computational and Mathematical Methods in Medicine 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 R 0 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 E 2 0 and an epithelial cell homeostasis one E 2 1 . It was proved that E 2 0 always exists and is unstable, whereas if the biological condition r B > δ B + μ B is satisfied, then E 2 1 exists. If the following conditions r * B /ðδ * B + μ * B Þ ≤ r B /ðδ B + μ B Þ and r * E /δ * E < r E /δ E are satisfied, the E 2 1 is locally asymptotically stable if R 0 < 1 and unstable when R 0 > 1. We note that the stability conditions have the following biological interpretation: r * E /δ * E < r E /δ 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 stra-tum intermedium cells. A similar interpretation can be made of the condition r * Using an elegant construction of a Lyapunov function, we obtained conditions on parameters ( to prove the global asymptotic stability of the epithelial cell homeostasis equilibrium. We observe that under the following conditions on the parameters r B = r * B , , and δ * E > δ E , the conditions of the Theorem 11 continue to be satisfied. Furthermore, it was shown that when R 0 = 1, E 2 1 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 E 2 2 , 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 equilibrium E 2 2 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 Vð 0Þ = 100. We also performed simulations with different ini-tial values in Vð0Þ, specifically with Vð0Þ = 10 and Vð0Þ = 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 PAP-CLEAR 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 (Figure 11) should be considered together with simulated outputs of virus dynamics models and the possible

20
Computational and Mathematical Methods in Medicine 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.

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, when R 0 = 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 E 2 2 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 R 0 = 1, so the qualitative analysis of the local and global stabil-ities of these points for values of R 0 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 point E 2 2 . 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 (γ).

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this manuscript.

Acknowledgments
This work has been partially supported by the CONACYT-Ciencia de Frontera "Ecuaciones de evolución, sus estados estacionarios y su comportamiento asintótico con aplicaciones en Física y Biología" (project N. 684340). Sierra-Rojas is indebted to CONACYT for fellowship 2018-000068-02NACF-0586 that enabled him to pursue graduate studies for the degree of Maestría en Matemáticas Aplicadas. The authors also wish to thank Anayeli Cisneros-Ines for their constructive comments and suggestions. We dedicate this paper to the memory of Raúl Peralta-Rodriguez (1978-2020) and Mesuma K. Atakishiyeva (1939-2021), who had always been generous to us.