Abstract
The aim of this paper is to model the dynamics of the human papillomavirus (HPV) in cervical epithelial cells. We developed a mathematical model of the epithelial cellular dynamics of the stratified epithelium of three (basale, intermedium, and corneum) stratums that is based on three ordinary differential equations. We determine the biological condition for the existence of the epithelial cell homeostasis equilibrium, and we obtain the necessary and sufficient conditions for its global stability using the method of Lyapunov functions and a theorem on limiting systems. We have also developed a mathematical model based on seven ordinary differential equations that describes the dynamics of HPV infection. We calculated the basic reproductive number () 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 and that in this case, the system presents a forward bifurcation, which shows the existence of an infected equilibrium point when . We also study the solutions numerically (i.e., viral kinetic in silico) when . 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.
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 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.
The model is given by , , and , which denote stratum basale, stratum intermedium, and stratum corneum formed by uninfected cells, respectively. The cells of proliferate (or lateral differentiation) at rate considering a logistic growth, with a carrying capacity ; these have suprabasal differentiation to at rate and die at rate . The dynamics of is generated by suprabasal differentiation of at a rate . This population increases by cell proliferation, which is governed by logistic growth at a rate and a carrying capacity at rate . These decrease by suprabasal differentiation at a rate , and we assume that there is no natural death. Finally, the dynamics of are generated by suprabasal differentiation of 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 We consider the positivity and boundedness of the solutions of system (1).
Theorem 1. Given the initial conditions , then all solutions of system (1) are positive.
Proof. At first, we will prove the positivity of . Let be the solution that satisfies the initial condition . Assume that the solution is not always positive; i.e., there exists a , such as . By Bolzano’s theorem, there exists a , such as . Let be the initial time, such as , and then . Note that if some , , then . Then, any solution with will satisfy . By uniqueness of solutions, we have that if , then will remain positive . Therefore, leads to a contradiction. Hence, is nonnegative for all . Now, we will prove the positivity of . is not always positive; i.e., there exists , such as . By Bolzano’s theorem, there exists , such as . Let . By the second equation (1), if , then implies that is increasing at . Therefore, will be negative for values near to , that is, a contradiction. Finally, we will prove the positivity of . From the third equation (1), we obtain the following inequality . Integrating, we obtain the solution ; therefore, .
Theorem 2. Let be the solution of model (1) with the initial conditions , , and . Then, , , and are all bounded for all at which the solution exists.
Proof. Let (, , ) be any solution with nonnegative initial conditions. We define a function The time derivative along a solution of (2) is It follows that where . Thus, . Therefore, , , and are all bounded for all .
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 always exists. We get the “epithelial cell homeostasis” equilibrium , where
The following inequality is a biological condition for the maintenance of cell homeostasis of the epithelium.
The Jacobian matrix of (1) at is
The eigenvalues of are , , and . By biological condition, , then has one positive eigenvalue. Thus, the trivial equilibrium is unstable.
The Jacobian matrix for the equations of (1) at is
Using the following identities and , we have
The characteristic polynomial of is
The eigenvalues of are , , and . Clearly, all of the eigenvalues are negative. Consequently, the epithelial cell homeostasis equilibrium is locally asymptotically stable.
We then arrive at the following theorem.
Theorem 4. Assume that the biological condition is satisfied. The trivial equilibrium always exists and is unstable. The epithelial cell homeostasis equilibrium 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 subsystem of model (1), which is independent of the variables. We construct the following Volterra-type Lyapunov function [26] for the subsystem
The function is defined, continuous, and positive definite for all . Also, the global minimum occurs at , and therefore, is a Lyapunov function. First, we calculate the time derivative of computed along solutions of the first two equations of the system (1), given by the expression
Note that
Because , the Lyapunov stability theorem [27] implies that the equilibrium is globally asymptotically stable in and .
From the third equation in (1), we can formally solve to obtain
By L’Hospital’s rule, we obtain . An application of Lemma 1 [28] shows that the epithelial cell homeostasis equilibrium of model (1) is globally asymptotically stable in the interior of . We then have the following theorem.
Theorem 5. If , then the epithelial cell homeostasis equilibrium 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 , , , and , which denote the cells of the stratum basale, stratum intermedium, stratum corneum infected, and the viral load, respectively. , , and are denoted in the same way as in the previous section.
Figure 4 shows that the dynamics of 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 , and a carrying capacity at rate . These decrease by suprabasal differentiation at rate and death cellular at rate . The dynamics of are generated by suprabasal differentiation of at a rate . This population increases by cell proliferation, governed by full logistic growth at a rate , and a carrying capacity at rate . These decrease by suprabasal differentiation at a rate , and we assume that there is no natural death. The dynamics of are generated by suprabasal differentiation of at rate . There is no proliferation of . This population decreases by desquamation at rate . Finally, increases by rupture of cell membrane of at rate , and they decline at rate . This is summarized in the following nonlinear ODE system:
3.1. Positivity, Boundedness, Equilibria, Basic Reproductive Number, and Local Stability
Let set be We consider the positivity and boundedness of the solutions of system (18).
Theorem 6. Given the initial conditions , , , , , and , 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 7. Let be the solution of model (18) with the initial conditions , , , , , and . Then, , , , , , , and are all bounded for all at which the solution exists.
Proof. Let be any solution with nonnegative initial conditions. We define a function The time derivative along a solution of (19) is It follows that where . Thus, . Therefore, , , , , , , and are all bounded for all .
Remark 8. The total number of cervical epithelial cells, both uninfected and infected, and viral load is bounded by a weighted sum of the stratum basale and stratum intermedium carrying capacities, where the weights are the proliferation rates of healthy and infected cells divided by four times of minimum of the death, differentiation, and viral clearance rates.
The model (18) always have a trivial equilibrium , and a “epithelial cell homeostasis” equilibrium with the same coordinates , , and 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
To compute the basic reproductive number , we use the next-generation operator introduced by van den Driessche and Watmough [29]. The Jacobian matrix of this subsystem at the infection-free equilibrium is decomposed as , where is the viral transmission part and describe the transition terms associated with the model (18). These quantities are given, respectively, by
It follows that the basic reproductive number , where is the spectral radius, is given by
The parameter has an interesting biological meaning: it is the sum of average numbers of secondary infected cells produced by a single infected cell in a population of epithelial basal cells, by direct basal cell-to-HPV contact and infected basal cell division, respectively.
Remark 9. If is taken as a new infection, we can get another basic reproductive number as follows: where When , it is easy to check that is equivalent to .
3.1.2. Local Stability of
The Jacobian matrix (22) evaluated at the equilibrium point becomes
By the biological condition, one of its eigenvalues is positive: while the other six are
Thus, is unstable. This result can be summarized as follows.
Theorem 10. Assume that the biological condition is satisfied. The trivial equilibrium always exists and is unstable.
3.1.3. Local Stability of
The Jacobian matrix of system (18), evaluated in the infection-free equilibrium point , takes the form where
Using the following identities and , we have
By biological condition, , then and . Furthermore, if
then and , respectively.
On the other hand, the characteristic polynomial of (31) is given by where
Note that in (43), when and ; consequently, one of the roots of (39) is zero. Therefore, in this case the equilibrium point has a zero eigenvalue, that is, it is no-hyperbolic.
Three eigenvalues of characteristic polynomial (39) are , and . 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 , , , and are positive and the relations hold. In order to show these relations, we will adopt the following notation:
Note that quantities and , since and are positive parameters, and is also positive. By the inequalities (37) and (38), we have that is positive and is nonnegative. Thus, in terms of this new notation, the quantities (40)–(43) can be rewritten as
We note that the coefficients , , and are positive, while if (or equivalently ) and also, is positive.
Thus, the condition of Ruth-Hurwitz (44) can be written as
First of all notice that from (46), and In this way, taking into account that and , from the above inequalities, we get and This result allows us to establish that
Expanding the left hand side of (52), this can be rewritten as
In (54), since and are positive definite, since that and and since that . We would have on the right hand side of (54) that 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 , 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 , we find that in polynomial (39), its eigenvalues , , and are negative. Additionally, it has also been shown that the four coefficients , , and are positive and the two conditions and are satisfied when and are fulfilled; in particular, is positive if it also holds that . Thus, if these three conditions are met: , , and , then the eigenvalues of the fourth order polynomial in (39) will have a negative real part, and consequently, the point will be asymptotically stable. On the other hand, by Descartes’ rule of signs, if , then , and the full polynomial (39) will have a positive eigenvalue, in which case will be unstable. These results can also be summarized in the following theorem.
Theorem 11. Assume that the following conditions , and are satisfied. The infection-free equilibrium of system (18) is locally asymptotically stable for and unstable for .
As already mentioned, is nonhyperbolic when . It will be shown later that when this happens, a bifurcation occurs.
3.2. Global Stability of the Infection-Free Equilibrium
We obtained some conditions on global stability of the infection-free equilibrium of the system (18).
Theorem 12. Assume that , , , , , and . If , then the infection-free equilibrium of system (18) is globally asymptotically stable in .
Proof. We construct the following Lyapunov function for the system (18): where the following functions have been defined as The function is defined, continuous, and positive definite for all . Also, the global minimum occurs at , and therefore, is a Lyapunov function. First, we calculate the time derivative of . Using we have After several calculations, we have Second, we calculate the time derivative of . Considering we obtain Writing in terms of and considering that , we find Taking into account that we obtain Third, we calculate the time derivative of . Using , we have Considering we find Finally, considering expressions (66), (70), and (73), we obtain If , then . Note that if and only if , , , and or , , , , and . Therefore, the largest compact invariant set in is the singleton . By the classical LaSalle invariance principle (Theorem 5.3 of [30]), is globally asymptotically stable in if .
3.3. Existence of a Forward or Backward Bifurcation
It has already been mentioned that when , in the characteristic equation (39), vanishes given rise to a zero eigenvalue and the infection-free equilibrium point 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 , according to definition of , this parameter takes the value
The Jacobian matrix of system (75) evaluated at the equilibrium point , when , is where the quantities , , , and , 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 , satisfying , is
By algebraic calculations, it can be shown that the required nonzero second-order partial derivatives evaluated at , which are contained in the quantities and given in Theorem 4.1 of [31], are
Thus, the quantities and take the form where , , , , and are the first, second, third, fourth, and seventh components of eigenvector (81), while 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 or . On the other hand, since , and . 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 , when , 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 , 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 at the infection-free equilibrium point . The second is the determination, when , of the regions of existence and stability of the infected equilibrium point . 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 of epithelial tissue, where the length is and the height is . We estimate and 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: , , , , and . For parameters and , they were proposed with values close to the proliferation rates of uninfected stratum basale and stratum intermedium cells, respectively.
4.2. Existence of a Forward Bifurcation at Infection-Free Equilibrium Point
To determine the type of bifurcation that occurs, when , at the equilibrium point , were considered , and all other required quantities from Table 1. With these values it is found, by (85), that , which is positive. Furthermore, it is found, according to (78) and (86), respectively, that and . Since for the typical values of the parameters considered, we find that ; consequently, according to Theorem 13, at the infection-free equilibrium point , when , 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 of the model (18) when 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 using the bisection method. We calculate the eigenvalues of the Jacobian matrix (22) evaluated at 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 not exist, or be locally asymptotically stable or unstable, respectively.
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. and such that 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 and such that shows that a belongs to the instability region of and to the asymptotically stable region of , 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 days, and finally, the infection is resolved.
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 and and Figure 9 with and such that . The ordered pair for Figure 8 belongs to the region of asymptotic stability of of Figure 5 and for Figure 9 belongs to the region of instability of 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 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 and such that . The ordered pair belongs to the region of asymptotic stability of 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.
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, , that depends differentiable on a parameter, , 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 . As we have an explicit formula for , we derive an analytical expression for the sensitivity of , to each of the parameters described in .
Figure 11 shows that all the parameters have either positive or negative effects on the . According to the sensitivity indices illustrated in Figure 11, we observe that the parameters, namely, , , , and , are the most positively sensitive parameters, respectively. This implies that decreasing the values of these parameters will decrease . Parameters , , , and are the most negatively sensitive parameters. This implies that increasing the values of these parameters will decrease .
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 for the existence of the epithelial cell homeostasis equilibrium , and consequently, we include cellular dynamics of the cervical epithelium. We have proven that trivial equilibrium 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 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 and an epithelial cell homeostasis one . It was proved that always exists and is unstable, whereas if the biological condition is satisfied, then exists. If the following conditions and are satisfied, the is locally asymptotically stable if and unstable when . We note that the stability conditions have the following biological interpretation: , 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 . Using an elegant construction of a Lyapunov function, we obtained conditions on parameters (, , , , and ) to prove the global asymptotic stability of the epithelial cell homeostasis equilibrium. We observe that under the following conditions on the parameters , , , , and , the conditions of the Theorem 11 continue to be satisfied. Furthermore, it was shown that when , 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 , 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 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 . We also performed simulations with different initial values in , specifically with and , 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 , 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 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 , so the qualitative analysis of the local and global stabilities of these points for values of 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 . 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.