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 (R 0) 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 R 0 = 1 and that in this case, the system presents a forward bifurcation, which shows the existence of an infected equilibrium point when R 0 > 1. We also study the solutions numerically (i.e., viral kinetic in silico) when R 0 > 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.

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 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.

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 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].

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 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.

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 Figure 1). 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 by B(t), E(t), and C(t), which denote stratum basale, stratum intermedium, and stratum corneum formed by uninfected cells, respectively. The cells of B(t) proliferate (or lateral differentiation) at rate r considering a logistic growth, with a carrying capacity K; these have suprabasal differentiation to E(t) at rate δ and die at rate μ. The dynamics of E(t) is generated by suprabasal differentiation of B(t) at a rate δ. This population increases by cell proliferation, which is governed by logistic growth at a rate r and a carrying capacity at rate K. These decrease by suprabasal differentiation at a rate δ, and we assume that there is no natural death. Finally, the dynamics of C(t) are generated by suprabasal differentiation of E(t) at rate δ and shed from the epithelial tissue at a rate μ. We arrive at the following mathematical model:

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

Let set Ω be Ω = {(B, E, C) ∈ ℝ+3 : B ≥ 0, E ≥ 0, C ≥ 0} ⊂ ℝ3. We consider the positivity and boundedness of the solutions of system (1).

Theorem 1 .

Given the initial conditions B(0) > 0, E(0) > 0, C(0) > 0, then all solutions of system (1) are positive.

Proof

At first, we will prove the positivity of B(t). Let B(0) be the solution that satisfies the initial condition B(t) > 0. Assume that the solution is not always positive; i.e., there exists a t0′ ∈ ℝ+, such as B(t0′) < 0. By Bolzano's theorem, there exists a t1 ∈ (0, t0′), such as B(t1) = 0. Let t0 ∈ ℝ+ be the initial time, such as B(t0) = 0, and then dB(t)/dt = 0. Note that if some t ≥ 0, B(t) = 0, then dB(t)/dt = 0. Then, any solution with B(0) = 0 will satisfy B(t) = 0∀t > 0. By uniqueness of solutions, we have that if B(0) > 0, then B(t) will remain positive ∀t > 0. Therefore, B(t0) = 0 leads to a contradiction. Hence, B(t) is nonnegative for all t > 0. Now, we will prove the positivity of E(t). E(t) is not always positive; i.e., there exists t0′ ∈ ℝ+, such as E(t0′) < 0. By Bolzano's theorem, there exists t1 ∈ (0, t0′), such as E(t1) = 0. Let t0 = min{t | E(t) = 0}. By the second equation (1), if E(t0) = 0, then dE(t0)/dt = δB > 0 implies that E is increasing at t = t0. Therefore, E(t) will be negative for values t < t0 near to t0, that is, a contradiction. Finally, we will prove the positivity of C. From the third equation (1), we obtain the following inequality dC(t)/dt ≥ −μC(t). Integrating, we obtain the solution C ≥ C(0)e−; therefore, C ≥ 0.

Theorem 2 .

Let (B(t), E(t), C(t)) be the solution of model (1) with the initial conditions B(0) > 0, E(0) > 0, and C(0) > 0. Then, B(t), E(t), and C(t) are all bounded for all t ≥ 0 at which the solution exists.

Let (B(t), E(t), C(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{μ, δ/n, μ}. Thus, limsupB(t) ≤ (rK + rK)/4η. Therefore, B(t), E(t), and C(t) are all bounded for all t ≥ 0.

Remark 1 .

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 E01 = (0, 0, 0) always exists. We get the “epithelial cell homeostasis” equilibrium E11 = (B∗, E∗, C∗), where

The following inequality r > δ + μ is a biological condition for the maintenance of cell homeostasis of the epithelium.

The Jacobian matrix of (1) at E01 is

The eigenvalues of J(E01) are λ1 = r − δ − μ, λ2 = r − δ, and λ3 = −μ. By biological condition, r > δ + μ, then J(E01) has one positive eigenvalue. Thus, the trivial equilibrium is unstable.

The Jacobian matrix for the equations of (1) at E11 is

Using the following identities r − rB∗/K − (δ + μ) = 0 and r − rE∗/K − δ = −δB∗/E∗, we have

The characteristic polynomial of J(E11) is

The eigenvalues of P(λ11) are λ1,11 = −B∗/K, λ1,21 = −δB∗/E∗ − rE∗/K, and λ1,31 = −μ. 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 3 .

Assume that the biological condition r > δ + μ is satisfied. The trivial equilibrium E01 = (0, 0, 0) always exists and is unstable. The epithelial cell homeostasis equilibrium E11 = (B∗, E∗, C∗) 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 the B − E subsystem of model (1), which is independent of the C variables. We construct the following Volterra-type Lyapunov function [26] for the B − E subsystem

The function U(t) is defined, continuous, and positive definite for all B, E > 0. Also, the global minimum U(B, E) = 0 occurs at (B∗, E∗), 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

Because dU(t)/dt ≤ 0, the Lyapunov stability theorem [27] implies that the (B∗, E∗) equilibrium is globally asymptotically stable in ℝ+2 and limE(t) = E∗.

From the third equation in (1), we can formally solve to obtain

By L'Hospital's rule, we obtain limC(t) = δE/μ = δE∗/μ = C∗. An application of Lemma 1 [28] shows that the epithelial cell homeostasis equilibrium E11 = (B∗, E∗, C∗) of model (1) is globally asymptotically stable in the interior of Ω. We then have the following theorem.

Theorem 4 .

If r > δ + μ, then the epithelial cell homeostasis equilibrium E11 = (B∗, E∗, C∗) of system (1) is globally asymptotically stable in the interior of Ω.

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 by B(t), E(t), C(t), and V(t), which denote the cells of the stratum basale, stratum intermedium, stratum corneum infected, and the viral load, respectively. B(t), E(t), and C(t) are denoted in the same way as in the previous section.

Figure 4 shows that the dynamics of B(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∗, and a carrying capacity at rate K. These decrease by suprabasal differentiation at rate δ∗ and death cellular at rate μ∗. The dynamics of E(t) are generated by suprabasal differentiation of B(t) at a rate δ∗. This population increases by cell proliferation, governed by full logistic growth at a rate r∗, and a carrying capacity at rate K. These decrease by suprabasal differentiation at a rate δ∗, and we assume that there is no natural death. The dynamics of C(t) are generated by suprabasal differentiation of E(t) at rate δ∗. There is no proliferation of C(t). This population decreases by desquamation at rate μ∗. Finally, V(t) increases by rupture of cell membrane of C(t) at rate σ, and they decline at rate γ. This is summarized in the following nonlinear ODE system:

Figure 4

Diagram of viral dynamics when a HPV infection occurs.

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

Let set Ω be Ω = {(B, B, E, E, C, C, V) ∈ ℝ+7 : B ≥ 0, B ≥ 0, E ≥ 0, E ≥ 0, C ≥ 0, C ≥ 0, V ≥ 0} ⊂ ℝ7. We consider the positivity and boundedness of the solutions of system (18).

Theorem 5 .

Given the initial conditions B(0) > 0, B(0) ≥ 0, E(0) > 0, E(0) ≥ 0, C(0) > 0, C(0) ≥ 0 and V(0) > 0, then all solutions of system (18) are positive.

The proof of Theorem 1 is very similar to the proof of Theorem 6. For this reason, we omit the proof.

Theorem 6 .

Let (B(t), B(t), E(t), E(t), C(t), C(t), V(t)) be the solution of model (18) with the initial conditions B(0) > 0, B(0) ≥ 0, E(0) > 0, E(0) ≥ 0, C(0) > 0, C(0) ≥ 0 and V(0) > 0. Then, B(t), B(t), E(t), E(t), C(t), C(t), and V(t) are all bounded for all t ≥ 0 at which the solution exists.

Let (B(t), B(t), E(t), E(t), C(t), C(t), V(t)) be any solution with nonnegative initial conditions. We define a function

The time derivative along a solution of (19) is

It follows that where η = min{μ, μ∗, δ/(n + 1), δ∗/(n + 1), μ, μ∗/n, γ}. Thus, limsupB(t) ≤ ((r + r∗)K + (r + r∗)K)/4η. Therefore, B(t), B(t), E(t), E(t), C(t), C(t), and V(t) are all bounded for all t ≥ 0.

Remark 2 .

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 = (B∗, 0, E∗, 0, C∗, 0, 0) with the same coordinates B∗, E∗, and C∗ 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 where

3.1.1. Basic Reproductive Number R0

To compute the basic reproductive number R0, 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

It follows that the basic reproductive number R0 = ρ(FV−1), where ρ is the spectral radius, is given by

The parameter R0 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 3 .

If r∗(1 − E∗/K) is taken as a new infection, we can get another basic reproductive number as follows: where

When r∗ = 0, it is easy to check that R0 is equivalent to .

3.1.2. Local Stability of E02

The Jacobian matrix (22) evaluated at the equilibrium point E02 becomes

By the biological condition, one of its eigenvalues is positive: while the other six are

Thus, E02 is unstable. This result can be summarized as follows.

Theorem 7 .

Assume that the biological condition r > δ + μ is satisfied. The trivial equilibrium E02 = (0, 0, 0, 0, 0, 0, 0) always exists and is unstable.

3.1.3. Local Stability of E12

The Jacobian matrix of system (18), evaluated in the infection-free equilibrium point E12, takes the form where

Using the following identities r − rB∗/K − (δ + μ) = 0 and r − rE∗/K − δ = −δB∗/E∗, we have

By biological condition, r > δ + μ, then J11 < 0 and J33 < 0. Furthermore, if

then J22 ≤ 0 and J44 < 0, respectively.

On the other hand, the characteristic polynomial of (31) is given by where

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 = −μ < 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 hold. In order to show these relations, we will adopt the following notation:

Note that quantities A > 0 and C > 0, since γ and μ∗ 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 a1, a2, and a3 are positive, while if J44 < 0 (or equivalently r∗/δ∗ < r/δ) and R0 < 1 also, a4 is positive.

Thus, the condition of Ruth-Hurwitz (44) can be written as

First of all notice that from (46), A2 = γ2 + 2C + (μ∗)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

Expanding the left hand side of (52), this can be 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 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

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 > δ + μ, we find that in polynomial (39), its eigenvalues λ1,22 = J33, λ1,32 = J11, and λ1,12 = −μ 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 r∗/(δ∗ + μ∗) ≤ r/(δ + μ) and r∗/δ∗ < r/δ are fulfilled; in particular, a4 is positive if it also holds that R0 < 1. Thus, if these three conditions are met: r∗/(δ∗ + μ∗) ≤ r/(δ + μ), r∗/δ∗ < r/δ, 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 8 .

Assume that the following conditions r > δ + μ, r∗/(δ∗ + μ∗) ≤ r/(δ + μ), and r∗/δ∗ < r/δ 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 9 .

Assume that r > δ + μ, r = r∗, r = r∗, δ = δ∗, μ = μ∗, and δ∗ ≥ δ. If R0 ≤ 1, then the infection-free equilibrium E12 = (B∗, 0, E∗, 0, C∗, 0, 0) of system (18) is globally asymptotically stable in Ω.

We construct the following Lyapunov function for the system (18): where the following functions have been defined as

The function W(t) is defined, continuous, and positive definite for all B, B, E, E, C, C, V ≥ 0. Also, the global minimum W(B, B, E, E, C, C, V) = 0 occurs at E12 = (B∗, 0, E∗, 0, C∗, 0, 0), and therefore, W is a Lyapunov function. First, we calculate the time derivative of W1(t).

Using r − (δ + μ) = (r/K)B∗, we have

After several calculations, we have

Second, we calculate the time derivative of W2(t).

Considering r = (r/K)E∗ + δ − δ(B∗/E∗), we obtain

Writing βσδδ∗B∗/(γμ∗(δ∗ − r(1 − E∗/K))) in terms of R0 and considering that βσδδ∗B∗/(γμ∗(δ∗ − r(1 − E∗/K))) = (δ + μ)(R0 − (r/(δ + μ))(1 − B∗/K)), we find

Taking into account that r(1 − B∗/K) = δ + μ, we obtain

Third, we calculate the time derivative of W3(t).

Using δE∗ = μC∗, we have

Considering r(1 − B∗/K) = δ + μ, we find

Finally, considering expressions (66), (70), and (73), we obtain

If R0 ≤ 1, then dW/dt ≤ 0. Note that dW/dt = 0 if and only if B = B∗, B = 0, E = E∗, and E = 0 or R0 = 1, B = B∗, B = 0, E = E∗, and E = 0. Therefore, the largest compact invariant set in {(B(t), B(t), E(t), E(t), C(t), C(t), V(t)): dW/dt = 0} is the singleton {E12}. By the classical LaSalle invariance principle (Theorem 5.3 of [30]), E12 is globally asymptotically stable in Ω if R0 ≤ 1.

3.3. Existence of a Forward or Backward Bifurcation

It has already been mentioned that when R0 = 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 which is in terms of the new variables and the new parameters

We will consider that β is the bifurcation parameter and that when R0 = 1, according to definition of R0, this parameter takes the value

The Jacobian matrix of system (75) evaluated at the equilibrium point E12, when β = β∗, is where the quantities J11, J22, J33, and J44, 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, satisfying v · w = 1, is

By algebraic calculations, it can be shown that the required nonzero second-order partial derivatives evaluated at E12, 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 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 10 .

The infection-free equilibrium point E12, when R0 = 1, presents a backward bifurcation if β∗ < β, while it has a forward bifurcation if β∗ > β, where

It is important to note that, as a consequence of the previous result, when R0 > 1, there is a family of asymptotically stable infected equilibrium points, which we will denote as , constituting the upper branch of these types of bifurcations.

4. Viral Kinetic In 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 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 μm2 of epithelial tissue, where the length is 7282.5 μm and the height is 355 μm. We estimate K and K 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(0) = 1000, B(0) = E(0) = C(0) = 0, E(0) = 617, C(0) = 226, and V(0) = 100. For parameters r∗ and r∗, 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.

Parameter	Description	Value	Unit	Ranges	Ref
r B	Uninfected stratum basale cell proliferation rate	0.07	day−1	[0.03; 0.07]	[25]
r B ∗	Infected stratum basale cell proliferation rate	0.048	day−1	[−; −]	Fixed
r E	Uninfected stratum intermedium cell proliferation rate	0.039	day−1	[0.02; 1]	[25]
r E ∗	Infected stratum intermedium cell proliferation rate	0.04	day−1	[−; −]	Fixed
δ B	Uninfected cell differentiation rate from stratum basale cell to stratum intermedium	0.044	day−1	[−; −]	[33]
δ B ∗	Infected cell differentiation rate from stratum basale cell to stratum intermedium	0.05	day−1	[−; −]	Fixed
δ E	Uninfected cell differentiation rate from stratum intermedium cell to stratum corneum	0.099	day−1	[0.02; 1]	[25]
δ E ∗	Infected cell differentiation rate from stratum intermedium cell to stratum corneum	0.118	day−1	[0; 5]	[25]
μ B	Uninfected basal cell natural death rate	0		[−; −]	[34]
μ B ∗	Infected basal cell death rate	0		[−; −]	[34]
μ C	Desquamation rate of uninfected stratum corneum cells	0.27	day−1	[0.2; 1]	[25]
μ C ∗	Desquamation rate of infected stratum corneum cells	0.27	day−1	[−; −]	Fixed
K B	Stratum basale cell carrying capacity	2693	cells	[1443; 13465]	[32]
K E	Stratum intermedium cell carrying capacity	2114	cells	[553; 5010]	[32]
β	Viral transmission rate	Varies	virion−1 · day−1	[10−15; 10−5]	[25]
σ	Production rate of free virions	Varies	cell−1 · virion · day−1	[10; 103]	[11]
γ	Viral clearance rate	1.18	day−1	[0.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, when R0 = 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 β = 4.351 × 10−8. Since for the typical values of the parameters considered, we find that β∗ > β; 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 point E22 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 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 Figure 5, 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. Figure 7 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 Figure 8 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 Figure 10 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

Local sensitivity analysis is carried out in order to identify which parameters have the greatest influence on changes in the values of R0. As we have an explicit formula for R0, we derive an analytical expression for the sensitivity of R0, to each of the parameters described in R0.

Figure 11 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, r∗, β, σ, and δ∗, are the most positively sensitive parameters, respectively. This implies that decreasing the values of these parameters will decrease R0. Parameters μ∗, δ, δ∗, 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 number R0.

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 r > δ + μ for the existence of the epithelial cell homeostasis equilibrium E11 = (B∗, E∗, C∗), 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 r > δ + μ is satisfied, then E12 exists. If the following conditions r∗/(δ∗ + μ∗) ≤ r/(δ + μ) and r∗/δ∗ < r/δ 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: r∗/δ∗ < r/δ, 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 r∗/(δ∗ + μ∗) ≤ r/(δ + μ). Using an elegant construction of a Lyapunov function, we obtained conditions on parameters (r = r∗, r = r∗, δ = δ∗, μ = μ∗, and δ∗ ≥ δ) to prove the global asymptotic stability of the epithelial cell homeostasis equilibrium. We observe that under the following conditions on the parameters r = r∗, r = r∗, δ = δ∗, μ = μ∗, and δ∗ > δ, 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 equilibrium E22 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 initial 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 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 (Figure 11) 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 (δ∗) 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, when R0 = 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 point E22. 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 (γ).