1. INTRODUCTION

Human papilloma virus (HPV) is the most common sexually transmitted infection that can affect the epithelial tissue, leading to the formation of lesion tissue - dysplasia and cancer [11], [12], [21], [28]. The significance of this medical issue has driven extensive re-search on HPV epidemic models over the past few decades. These models address the problem at three distinct levels: (i) social level (investigating the transmission of HPV between individuals and assessing the efficacy of HPV vaccination within populations, clinical management strategies) [6], [7], [26], [29], [32], (ii) molecular or tissue level which include mathematical study of cell-HPV population dynamics [1], [2], [3], [31], [36], [42], and (iii) microbiological study of HPV-induced oncogenesis [16], [30], [38] - [40]. This paper is focused on the analysis of HPV transmission in epithelial tissue of vaccinated organism which consists of susceptible and infected cells, lesion tissue – precancerous (dysplasia) and cancer cells. The purpose of this study is understanding the conditions of resistance in a vaccinated organism to HPV infection and conditions leading to disease, specifically related to lesion epithelial tissue. Thus, we study the interaction and dynamics of healthy (susceptible) cells, HPV infected cells and lesion cells with HPV and does not separate the lesion cells on precancerous and cancer cells in contrast to the earlier studies [1] - [3], [31].

The autonomous SILV epidemic model is based on the SIPCV epidemic models [1], [31] and considers the age-structured dynamics of susceptible, infective, lesion cell populations and dynamics of unstructured HPV population. In SILV epidemic model we assume that the immune response of vaccinated organism is tolerant with respect to its own cells and, as a consequence, death rates of infected and lesion cells do not depend from the HPV abundance [11], [21], [28]. Interaction strength between susceptible cells and HPV is a product of the Lotka-Voltera incidence rate and result in the growth of infective cells [3]. Interaction strength between infected and lesion cells is an HPV density-dependent sigmoid or S-function which describes the impact of HPV load on malignant transformation of cells. We assume that lesion cells do not apply pressure on the tissues of organism and have no effect on the proliferation and mortality of other cells. The efficacy of multivalent HPV vaccine is primary associated with the magnitude of immune response. On one hand, the humoral immune response (B-lymphocytes) leads to the increase on the level (titer) of antibodies in blood, which creates conditions for: (i) blocking HPV entry into host cell [27], [39], and (ii) destroying viruses that are freely moving in the inter cellular space [40]. In the first case, the strength of HPV infection reduces, and in the second case, the HPV death rate increases. On the other hand, vaccine stimulates cytotoxic activity and increases the number of T-cells which destruct the infected cells (HPV host cells) and, consequently, slow down the HPV replication. Thus, from mathematical point of view impact of the efficacy of HPV vaccine are taken into account implicitly in coefficients of SILV epidemic model: HPV infection rate, budding number of HPV virions, apoptosis rate of infected cells and HPV death rate [1] - [3]. HPV vaccine efficacy-dependent coefficients define conditions of existence and asymptotic stability of disease free an endemic equilibrium of system (DFE and EE, respectively).

Conditions of existence and asymptotic stability of non-trivial equilibria of age-structured epidemic models of population dynamics are studied in framework of the linear stability theory [3] - [5], [8], [9], [14], [20], [22] - [25], [34], [35], [41], [43]. In work [3] we studied the local asymptotic stability of EE, which is associated with the features of the course and possible treatment of an HPV induced disease. Although the influence of vaccination efficacy on the dynamical regimes of HPV - cellular populations is obvious, the implicit impact of HPV vaccine efficacy on the local asymptotical stability of DFE was not addressed in previous works. This paper is focused on the study of conditions of existence and local asymptotic stability of DFE of SILV epidemic model, which in turn is associated with the resistance of a vaccinated organism to the HPV induced disease. Asymptotic stability of DFE of epidemic SILV model means that epithelial tissue (population of susceptible cells) of vaccinated organism, being infected by HPV, eventually recovers and returns back to the stationary, healthy state. In this case the efficacy of HPV vaccine is sufficiently high and helps organism to overcome infectious disease and get rid the lesion tissue. The main indicator of the DFE asymptotic stability, basic reproduction number of infection \(R_0\), depends on the vaccine efficacy-dependent coefficients of the model. Thus, we can relate the vaccine efficacy with dynamical regime of system and resistance of vaccinated organism to HPV infection. On the other hand, asymptotic stability of EE means that HPV vaccine efficacy is sufficiently low, epithelial tissue of vaccinated organism is sick and the lesion tissue is asymptotically stable.

Overall, stability analysis of the autonomous age-structured SILV epidemic model helps us better understand the relationship between HPV vaccine efficacy-dependent model parameters and the asymptotic stability of DFE which is associated with vaccine-induced protection and resistance of the organism to infectious disease, and the asymptotic stability of EE which is associated with low efficacy of vaccine, infectious disease and cancer.

2. MODEL FORMULATION

We consider a SILV epidemic model that consists of susceptible cells (noninfected cells), infected cells (cells without changed by HPV morphology, CIN I and CIN II stages), lesion cells (dysplasia, CIN III stage and cancer cells) [11], [12], [21], [28] and human papilloma virus (HPV) that moves freely between cells. The age-specific densities of susceptible, infectious and lesion cells are denoted as S(a,t), I(a,t) and L(a,t). The quantity (size) of susceptible, infected and lesion cells subpopulations are denoted by \(N_S(t) = \int\limits_0^{a_d} S(a,t) da\),

\(N_I(t) = \int\limits_0^{a_d} I(a,t)da, N_L(t) = \int\limits_0^{a_d} L(a,t)da\), respectively, where cell's maximum lifespan is \(a_d\). The dynamics of cell subclasses (subpopulations) is described by the nonlinear age-structured model with age- dependent death (apoptosis) rates of susceptible \(d_s(a)\), infectious \(d_q(a) \geq d_s(a)\), and lesion cells \(d_l(a) \geq d_s(a)\) with the age reproductive windows of non-lesion cells \([a_r, a_m]\) and lesion cells \([a_c, a_g]\). Lesion cells differ from the susceptible cells in their lack of response to normal fertility control mechanism \((a_c < a_r, a_g < a_m)\) [11], [21], [28].

We assume that adaptive behavior of the HPV makes the natural immune response of host tolerant with respect to infectious and lesion cells that is their death rates do not depend from the HPV abundance [11], [12], [21], [28]. Host organism recognizes lesion cells (dysplasia and cancer) as its own and does not attempt to destroy them. Since viruses are not living things and cannot reproduce (multiply) until they enter a host cell, we use the unstructured model of HPV subpopulation dynamics. The dynamics of HPV quantity V(t)is governed by non-linear ODE with density-dependent death rate \(d_v(V(t)) = d_{v0}\sqrt{1 + d_{v1}V(t)}\) [3] (where \(d_{v0}, d_{v1} = const > 0\)). HPV death rate \(d_v(V(t))\) is an increasing function of V because increasing of HPV quantity changes the characteristics of intracellular space that result in the organism immune response through the activation of cell immunity (T-cells) and humoral immunity (B-lymphocytes) that leads to the elimination of HPV (i.e. monotone increasing of HPV death rate). On the other hand, the rate of change \(\frac{d(d_v(V))}{dV}\) is a decreasing function of V, approaching zero when \(V \to \infty\). That is the immune respond approaches the maximum value and \(d_v(V(t))\) changes slowly when HPV population gets larger.

The interaction strength between susceptible and HPV is a product of the Lotka-Voltera incidence rate \(\alpha V(t)S(a,t)\) (where \(\alpha = const > 0\) is a rate of infection) [1], [3], [31], [35] and result in the growth of infectious cells, which move partially to the lesion subclass with rate \(\delta(V(t))I(a,t)\). In our previous works [1], [3] the progression rate from infected to precancerous cells (dysplasia) is \(\delta = const > 0\). However, the HPV viral load, or the amount of HPV DNA present in the intercellular space, influences the progression of infect-ed cells to lesion cells because: (i) high HPV load and HPV persistence in tissue increases the expression of oncoproteins, which accelerate the disruption of infected cell proliferation cycle control and enhance them transformation to lesion cells [17], [37], (ii) high HPV load overwhelms the both types of immune response (B-lymphocytes and T-cells), leading to persistent of HPV-infection and increased risk of transformation from infected to lesion cells [18]. Thus, the progression rate from infected to lesion cells \(\delta(V(t))\) should be an increasing function of HPV quantity \(\delta(V(t))\) with saturation which has the zero derivative at origin (slow start at V(t) = 0). The most suitable function satisfying these conditions is a sigmoid function (or S-function) \(\delta(V(t)) = \frac{\delta_d V^3(t)}{h_d^3 + V^3(t)}\) with saturation and half-saturation constants \(\delta_d > 0\), \(h_d > 0\), respectively.

Model considers also the partial apoptosis of infectious cells with rate \(n\int_{0}^{a_d}d_p(a)I(a,t)da\), where n is a mean budding number of HPV virions produced by one host cell (when viruses leave destroyed cells and ready to infect new susceptible cells), \(d_p(a)\) is an age- dependent death rates of infectious cells as a result of virus replication (HPV-induced apoptosis rate of infected cells) [3], [13]. These assumptions lead to the following age-structured epidemic model in domain \(Q = \{(a, t) | a \in (0, a_d), t \in (0, T)\}\):

\[\frac{\partial S(a,t)}{\partial t} + \frac{\partial S(a,t)}{\partial a} = -(d_s(a) + \alpha V(t))S(a,t), \tag{1}\]

\[\frac{\partial I(a,t)}{\partial t} + \frac{\partial I(a,t)}{\partial a} = -(d_q(a) + \delta(V(t)))I(a,t) + \alpha V(t)S(a,t), \tag{2}\]

\[\frac{\partial I(a,t)}{\partial t} + \frac{\partial I(a,t)}{\partial a} = -(d_q(a) + \delta(V(t)))I(a,t) + \alpha V(t)S(a,t), \qquad (2)\] \[\frac{\partial L(a,t)}{\partial t} + \frac{\partial L(a,t)}{\partial a} = -d_l(a)L(a,t) + \delta(V(t))I(a,t), \qquad (3)\]

\[\frac{\partial V(t)}{\partial t} = -d_v(V)V(t) + n \int_0^{a_d} d_p(a)I(a,t) da. \tag{4}\]

Equations (1) - (4) are completed by the boundary conditions and initial values:

\[S(0,t) = \int_{a_r}^{a_m} \beta_s(a)S(a,t)da, I(0,t) = \int_{a_r}^{a_m} \beta_s(a)I(a,t)da, L(0,t) = \int_{a_c}^{a_g} \beta_l(a)L(a,t)da,\] (5)

\[S(a,0) = S_0(a), I(a,0) = 0, L(a,0) = 0, V(0) = V_0.\] (6)

where \(\beta_s(a)\) is a fertility rate of the susceptible and infected cells, \(\beta_l(a)\) is a fertility rate of the lesion cells; \(S_0(a)\) is an initial density of susceptible cells, \(V_0\) is an initial value of HPV quantity. The model contains three sets of efficacy vaccine-dependent coefficients and parameters which are related implicitly to three factors of HPV vaccine impact on epithelial tissue. First coefficient is the rate of HPV infection \(\alpha\). Vaccine elicits a strong immune response that reduces the likelihood of HPV infection persistent and leads to a decrease in the \(\alpha\) coefficient. Second pair of coefficients are the HPV death rate constants \(d_{v0}\) and \(d_{v1}\). As it was mentioned above, vaccine-induced humoral immune response (B-lymphocytes) destroys viruses that are freely moving in the inter cellular space, increasing the HPV death rate and constants \(d_{v0}\) and \(d_{v1}\). Third parameter is a product \(nd_p(a)\) of the mean budding number of HPV virions produced by one host cell n and the HPV-induced apoptosis rate of infected cells \(d_p(a)\). HPV vaccine targets the viral oncoproteins which plays the key role in HPV ability to replicate and inhibits apoptosis in infected cells leading to decreasing of value of \(nd_p(a)\). We assume that all coefficients of System (1) - (6) are twice continuously differentiable functions and have the partial derivatives of the second order by all their arguments. In the next sections we analyse the existence and local asymptotic stability of all equilibria (stationary states) of System (1) - (6).

3. EQUILIBRIA OF SYSTEM (1) - (6)

It is easy to verify by direct substitution that trivial equilibrium (TE) (0,0,0,0) of System (1) – (6) always exists. The disease-free equilibrium (DFE) \((S_0^*(a),0,0,0)\) is a solution of stationary system (1) - (6):

\[S_0^*(a) = S_0^*(0) \exp\left(-\int_0^a d_s(s)ds\right) > 0.\] (7)

Integrating both sides of (7) with respect to a from 0 to \(a_d\), we can obtain another expression of DFE depending from the equilibrium value of susceptible cell subpopulation size \(N_{S_0}^* = \int_0^{a_d} S_0^*(a) da\):

\[S_0^*(a) = N_{S_0}^* \left( \int_0^{a_d} \exp\left(-\int_0^a d_s(s)ds\right) da \right)^{-1} \exp\left(-\int_0^a d_s(s)ds\right). \tag{8}\]

Plugging boundary condition (5) into (8), multiplying both sides of equation by \(\beta_s(a)\) and integrating them with respect to a from 0 to \(a_d\), after a little algebra, we arrive to the condition of existence of DFE:

\[R_S = \int_{a_r}^{a_m} \beta_s(a) \exp\left(-\int_0^a d_s(s)ds\right) da = 1.\] (9)

where \(R_S\) is a basic reproduction number of susceptible cells subpopulation. Equation (9) is the balance condition of stationary model which relates the susceptible cells death and fertility rates so that the number of newborn susceptible cells is equal to the number of dead susceptible cells of reproductive age. If Equation (9) holds and the initial value of the HPV population \(V_0=0\), the equilibrium density of susceptible cells is equal to the initial value \(S_0^*(a)=S_0(a)\) and the size of the susceptible sub-population does not change with time \(N_{S_0}^*=N_{S_0}=\int\limits_0^a S_0(a)da\). If \(V_0>0\) we cannot obtain the exact value of \(S_0^*(a)\), but we can assume that in this case the equilibrium \(S_0^*(a)\) depends also on the initial values \(S_0(a)\) and \(V_0\). This dynamical regime of subpopulations will be studied below in numerical experiments. Endemic equilibrium of system \((S^*(a),I^*(a),L^*(a),V^*)\), \(S^*(a)\geq 0\), \(I^*(a)\geq 0\), \(L^*(a)\geq 0\), \(V^*>0\) satisfies the stationary system:

\[\frac{dS^*(a)}{da} = -(d_s(a) + \alpha V^*)S^*(a), \tag{10}\]

\[\frac{dI^*(a)}{da} = -(d_q(a) + \delta(V^*))I^*(a) + \alpha V^* S^*(a), \tag{11}\]

\[\frac{dL^*(a)}{\partial a} = -d_l(a)L^*(a) + \delta(V^*)I^*(a), \tag{12}\]

\[0 = -d_v(V^*)V^* + nd_p^* \int_0^{a_d} \tilde{d}_p(a)I^*(a)da,\] (13)

\[S^*(0) = \int_{a}^{a_m} \beta_s(a) S^*(a) da, I^*(0) = \int_{a}^{a_m} \beta_s(a) I^*(a) da, L^*(0) = \int_{a}^{a_g} \beta_l(a) L^*(a) da.\] (14)

where \(d_p^* = \max_{a \in [0, a_d]} (d_p(a))\), \(\tilde{d}_p(a) = d_p(a)/d_p^*\). Integrating (10) with respect to a, and using (14), we obtain the integral equation:

\[S^*(a) = A_S \exp\left(-\int_0^a (d_s(s)ds + \alpha V^*)ds\right),\tag{15}\]

\[A_S = \int_a^{a_m} \beta_s(a) S^*(a) da. \tag{16}\]

Multiplying (15) by \(\beta_s(a)\), integrating it with respect to a from \(a_r\) to \(a_m\), and using (16) we arrive to the transcendental equation for positive constant \(V^*\):

\[W(V^*) = \int_{a_r}^{a_m} \beta_s(a) \exp\left(-\int_0^a (d_s(s) + \alpha V^*) ds\right) da = 1,\] (17)

where \(W(V^*)\) is an auxiliary function. Since \(W(V^*)>0\) and \(\frac{dW(V^*)}{dV^*}<0\) for all \(V^*>0\), \(\lim_{V^*\to\infty}W(V^*)=0\), the unique positive root \(V^*>0\) of (17) exists if and only if the basic reproduction number of the population of susceptible cells \(W(0)=R_S>1\) (9). Integrating (11) with respect to a, and using (14), (15), (16), we arrive to the integral equation

\[I^{*}(a) = A_{I} \exp\left(-\int_{0}^{a} (d_{q}(s) + \delta(V^{*}))ds\right) + \alpha V^{*} \int_{0}^{a} S^{*}(z) \exp\left(-\int_{z}^{a} (d_{q}(s) + \delta(V^{*}))ds\right) dz, \quad (18)\]

\[A_I = \int_{a_n}^{a_m} \beta_s(a) I^*(a) da. \tag{19}\]

Multiplying (18) by \(\beta_s(a)\), integrating it with respect to a from \(a_r\) to \(a_m\), and using (19), we arrive to the linear equation for unknown constants \(A_I\) and \(A_S\):

\[(1 - R_I) A_I = \alpha V^* A_S B_S, \tag{20}\]

\[B_S = \int_{a_r}^{a_m} \beta_s(a) \int_0^a \exp\left(-\int_0^z (d_s(s) + \alpha V^*) ds\right) \exp\left(-\int_z^a (d_q(s) + \delta(V^*)) ds\right) dz da. \tag{21}\]

Since \(A_I > 0\) and \(A_s > 0\), the right side of (20) is positive, we obtain another restriction for coefficients of System (1-6):

\[R_{I} = \int_{a_{r}}^{a_{m}} \beta_{s}(a) \exp\left(-\int_{0}^{a} (d_{q}(s) + \delta(V^{*}))ds\right) da < 1.\] (22)

Equation (22) relates infected cell death and fertility rates with infected cell progression to lesion cells so that the number of newborn infected cells is less than the number of infected cells moved to lesion subpopulation and dead infected cells of reproductive age. From (20) we obtain constant \(A_I\):

\[A_I = \alpha V^* (1 - R_I)^{-1} B_S A_S. \tag{23}\]

Multiplying (18) by \(nd_p^*\tilde{d}_p(a)\), integrating it with respect to a from 0 to \(a_d\) and plugging it into (12), we arrive to the second linear equation for constants \(A_I\) and \(A_S\):

\[d_v(V^*)V^* = nd_p^*B_IA_I + \alpha nd_p^*V^*C_SA_S,\] (24)

\[B_{I} = \int_{0}^{a_{d}} \tilde{d}_{p}(a) \exp\left(-\int_{0}^{a} (d_{q}(s) + \delta(V^{*}))ds\right) da, \tag{25}\]

\[C_S = \int_0^{a_d} \tilde{d}_p(a) \int_0^a \exp\left(-\int_0^z (d_s(s) + \alpha V^*) ds\right) \exp\left(-\int_z^a (d_q(s) + \delta(V^*)) ds\right) dz da. \tag{26}\]

Plugging (23) into (24) yields

\[A_s = d_v(V^*) \left(\alpha n d_p^*\right)^{-1} \left( (1 - R_I)^{-1} B_I B_S + C_S \right)^{-1}.\] (27)

Integrating (12) with boundary condition (14), after a little algebra we obtain solution

\[L^{*}(a) = (1 - R_{L})^{-1} \delta(V^{*}) \int_{a_{c}}^{a_{g}} \left( \beta_{l}(a) \int_{0}^{a} I^{*}(z) \exp\left(-\int_{z}^{a} d_{l}(s) ds\right) dz \right) da \exp\left(-\int_{0}^{a} d_{l}(s) ds\right) + \delta(V^{*}) \int_{0}^{a} I^{*}(z) \exp\left(-\int_{z}^{a} d_{l}(s) ds\right) dz, \quad (28)\]

\[R_L = \int_{a_c}^{a_g} \beta_l(a) \exp\left(-\int_0^a d_l(s)ds\right) da < 1.\] (29)

Equation (29) relates the death rates of lesion cells and the fertility rate so that the number of newborn lesion cells is less than the number of dead lesion cells of reproductive age.

Thus, the basic reproduction numbers of susceptible, infected and lesion cells \(R_S > 1\) (Equation (9)), \(R_I < 1\) (Equation (22)) and \(R_L < 1\) (Equation (29)) together define the condition of existence of EE which is given in the following Theorem.

Theorem 3.1. System (1) - (6) posses

(i) unique disease-free equilibrium (DFE) \((S_0^*(a), 0, 0, 0)\) if \(R_S = 1\) (9); (ii) unique endemic equilibrium (EE) \((S^*(a), I^*(a), L^*(a), V^*)\) if \(R_S > 1\) (9), \(R_I < 1\) (22) and \(R_L < 1\) (29). Endemic equilibrium is defined by \(V^*\) (16), \(S^*(a)\) (14), (26), \(I^*(a)\) (17), (22), \(L^*(a)\) (28), (29).

4. LOCAL ASYMPTOTIC STABILITY OF DFE \((S_0^*(a), 0, 0, 0)\)

Linearizing System (1) - (6) at the DFE \((S_0^*(a),0,0,0)\), we arrive to the system for perturbations \(\bar{\psi}_s(a,t)=\psi_s(a)\exp(\lambda t)\) for \(S_0^*(a)\), \(\bar{\psi}_i(a,t)=\psi_i(a)\exp(\lambda t)\) for \(I^*(a)=0\), \(\bar{\psi}_l(a,t)=\psi_l(a)\exp(\lambda t)\) for \(L^*(a)=0\) and \(\psi_v(t)=\psi_v\exp(\lambda t)\) for \(V^*=0\):

\[\frac{d\psi_s(a)}{da} = -(\lambda + d_s(a))\psi_s(a) - \alpha S_0^*(a)\psi_v, \tag{30}\]

\[\frac{d\psi_i(a)}{da} = -(\lambda + d_q(a))\psi_i(a) + \alpha S_0^*(a)\psi_v, \tag{31}\]

\[\frac{\partial \psi_l(a)}{\partial a} = -(\lambda + d_l(a))\psi_l(a),\tag{32}\]

\[0 = -(\lambda + d_{v0})\psi_v + nd_p^* \int_0^{a_d} \tilde{d}_p(a)\psi_i(a)da,\] (33)

\[\psi_s(0) = \int_{a_r}^{a_m} \beta_s(a)\psi_s(a)da, \psi_l(0) = \int_{a_r}^{a_m} \beta_s(a)\psi_l(a)da, \psi_l(0) = \int_{a_c}^{a_g} \beta_l(a)\psi_l(a)da.\] (34)

Integrating (30) and taking into account that \(R_S = 1\) (9) we have

\[\psi_s(a) = D_S \exp\left(-\int_0^a (\lambda + d_s(s))ds\right) - \alpha \psi_v \int_0^a S_0^*(z) \exp\left(-\int_z^a (\lambda + d_s(s))ds\right) dz, \tag{35}\]

\[D_S = -\alpha \psi_v \left( 1 - \int_{a_r}^{a_m} \beta_s(a) \exp\left( -\int_0^a (\lambda + d_s(s)) ds \right) da \right)^{-1}\] \[\times \int_{a_r}^{a_m} \beta_s(a) \int_0^a S_0^*(z) \exp\left( -\int_z^a (\lambda + d_s(s)) ds \right) dz da. \quad (36)\]

Linear Equation (36) relates unknown constants \(D_S\) and \(\psi_v\). Since \(R_S = 1\), from (36) it follows that \(\lambda \neq 0\) because in this case we have a trivial solution of System (30) - (34). Integrating (31) yields:

\[\psi_i(a) = D_I \exp\left(-\int_0^a (\lambda + d_q(s))ds\right) + \alpha \psi_v \int_0^a S_0^*(z) \exp\left(-\int_z^a (\lambda + d_q(s))ds\right) dz, \tag{37}\]

\[D_I = \int_a^{a_m} \beta_s(a)\psi_i(a)da. \tag{38}\]

Multiplying (37) by \(\beta_s(a)\), integrating it with respect to a from 0 to \(a_d\) after a little algebra we arrive to the linear equation for unknown \(D_I\) and \(\psi_v\):

\[\psi_v w_1(\lambda) - D_I w_2(\lambda) = 0, (39)\]

\[w_1(\lambda) = \alpha \int_{a_r}^{a_m} \left( \beta_s(a) \int_0^a S_0^*(z) \exp\left(-\int_z^a (\lambda + d_q(s)) ds\right) dz \right) da, \tag{40}\]

\[w_2(\lambda) = 1 - \int_{a_r}^{a_m} \beta_s(a) \exp\left(-\int_0^a (\lambda + d_q(s))ds\right) da.\] (41)

We assume that the expression in the parentheses of the second term of Equation (39) can't be equal to zero with some \(\lambda > 0\) because in this case we have a trivial solution of System (30) - (34). Multiplying (37) by \(nd_p^*\tilde{d}_p(a)\), integrating it with respect to a from 0 to \(a_d\) and plugging it into (30) we arrive to the another linear equation for \(D_I\) and \(\psi_v\):

\[\psi_v w_3(\lambda) - D_I w_4(\lambda) = 0, (42)\]

\[w_3(\lambda) = \lambda + d_{v0} - n\alpha d_p^* \int_0^{a_d} \left( \tilde{d}_p(a) \int_0^a S_0^*(z) \exp\left(-\int_z^a (\lambda + d_q(s)) ds\right) dz \right) da, \tag{43}\]

\[w_4(\lambda) = nd_p^* \int_0^{a_d} \tilde{d}_p(a) \exp\left(-\int_0^a (\lambda + d_q(s))ds\right) da. \tag{44}\]

Equating to zero determinant of System (39 - 40) we obtain the characteristic equation for \(\lambda\):

\[w_1(\lambda)w_4(\lambda) = w_2(\lambda)w_3(\lambda). \tag{45}\]

It is easy to verify that \(w_n(\lambda)\), n=1,...,4, are monotonic functions which possess the following properties for all \(\lambda > 0\):

\[w_1(\lambda) > 0, w_4(\lambda) > 0, \lim_{\lambda \to \infty} w_2(\lambda) = 1, \lim_{\lambda \to \infty} w_3(\lambda) = \infty,\] (46)

\[\frac{dw_1(\lambda)}{d\lambda} < 0, \frac{dw_4(\lambda)}{d\lambda} < 0, \frac{dw_2(\lambda)}{d\lambda} > 0, \frac{dw_3(\lambda)}{d\lambda} > 0.\] (47)

Since \(w_2(\lambda)\) and \(w_3(\lambda)\) are increasing functions and approach positive values at positive infinity, if \(w_2(0) \leq 0\) or/and \(w_3(0) \leq 0\) they always have a unique \(\lambda\)-intercept at \(\lambda \geq 0\). In this case the right side of characteristic equation (45) \((w_2(\lambda)w_3(\lambda))\) has at least one \(\lambda\)-intercept at \(\lambda \geq 0\), increases and approaches positive infinity when \(\lambda \to \infty\): \(\lim_{\lambda \to \infty} (w_2(\lambda)w_3(\lambda)) = \infty\). Next up, the positive decreasing function \((w_1(\lambda)w_4(\lambda))\) and the positive increasing from zero to positive infinity function \((w_2(\lambda)w_3(\lambda))\) always intersect, characteristic equation (45) always have at least one positive root \(\lambda^* > 0\). If \(w_2(0) > 0\) and \(w_3(0) > 0\) simultaneously, function \((w_2(\lambda)w_3(\lambda))\) does not have \(\lambda\)-intercept at \(\lambda \geq 0\), monotonically increases and approaches positive infinity when \(\lambda \to \infty\). Thus, if \(w_1(0)w_4(0) \geq w_2(0)w_3(0)\), positive decreasing function \((w_1(\lambda)w_4(\lambda))\) and the positive increasing function \((w_2(\lambda)w_3(\lambda))\) always intersect, characteristic equation (45) always have the unique non-negative root \(\lambda^* \geq 0\). If \(w_2(0) > 0\) and \(w_3(0) > 0\) simultaneously and \(w_1(0)w_4(0) < w_2(0)w_3(0)\), that is the basic reproduction number of infection [9], [10], [23]:

\[R_{0} = \frac{\alpha n d_{p}^{*}}{d_{v0}} \int_{a_{r}}^{a_{m}} \left( \beta_{s}(a) \int_{0}^{a} S_{0}^{*}(z) \exp\left(-\int_{z}^{a} d_{q}(s) ds\right) dz \right) da \int_{0}^{a_{d}} \tilde{d}_{p}(a) \exp\left(-\int_{0}^{a} d_{q}(s) ds\right) da\] \[\times \left(1 - \frac{\alpha n d_{p}^{*}}{d_{v0}} \int_{0}^{a_{d}} \left( \tilde{d}_{p}(a) \int_{0}^{a} S_{0}^{*}(z) \exp\left(-\int_{z}^{a} d_{q}(s) ds\right) dz \right) da \right)^{-1}\] \[\times \left(1 - \int_{a_{p}}^{a_{p}} \beta_{s}(a) \exp\left(-\int_{0}^{a} d_{q}(s) ds\right) da \right)^{-1} < 1, \quad (48)\] the positive decreasing function \((w_1(\lambda)w_4(\lambda))\) and the positive increasing function \((w_2(\lambda)w_3(\lambda))\) never intersect, characteristic equation (45) does not have the non-negative root \(\lambda^* \geq 0\).

Integrating Equation (32) we arrive to another characteristic equation for \(\lambda\):

\[\int_{a_{s}}^{a_{g}} \beta_{l}(a) \exp\left(-\int_{0}^{a} (\lambda + d_{l}(s))ds\right) da = 1.\] \[(49)\]

It is easy to verify that non-negative real root \(\lambda \ge 0\) of Equation (49) does not exist if \(R_L < 1\) (29) while it exists when \(R_L \ge 1\). The above results lead to the following Theorem.

Theorem 4.1. Let conditions of existence of the DFE \((S_0^*(a), 0, 0, 0)\) given in Theorem 3.1 (case (i)) hold. If \(w_2(0) > 0\) and \(w_3(0) > 0\) simultaneously (33), (41), the basic reproduction number of infection \(R_0 < 1\) (48) and \(R_L < 1\) (29), then DFE is locally asymptotically stable. Otherwise, if at least one of the following conditions holds:

• (i) \(w_2(0) \leq 0\) (41);
• (ii) \(w_3(0) \le 0\) (43);

(iii) \(w_2(0) > 0\) and \(w_3(0) > 0\) simultaneously and basic reproduction number of infection \(R_0 \ge 1\) (48); (iv) \(R_L \ge 1\) (29); then DFE \((S_0^*(a), 0, 0, 0)\) is unstable.

Remark 1. Using equation (10) in (48) yields the another expression of basic reproduction number of infection \(R_0\) depending from the parameter \(R_{vc} = \frac{\alpha N_{S_0}^* n d_p^*}{d_{v0}}\):

\[R_{0} = R_{vc} \left( \int_{0}^{a_{d}} \exp\left(-\int_{0}^{a} d_{s}(s)ds\right) da \right)^{-1} \int_{a_{r}}^{a_{m}} \left(\beta_{s}(a) \int_{0}^{a} \exp\left(-\int_{0}^{z} d_{s}(s)ds\right)\right) ds\] \[\times \exp\left(-\int_{z}^{a} d_{q}(s)ds\right) dz ds \int_{0}^{a_{d}} \tilde{d}_{p}(a) \exp\left(-\int_{0}^{a} d_{q}(s)ds\right) da \left(1 - R_{vc} \left(\int_{0}^{a_{d}} \exp\left(-\int_{0}^{a} d_{s}(s)ds\right) da\right)^{-1} \right)\] \[\times \int_{0}^{a_{d}} \left(\tilde{d}_{p}(a) \int_{0}^{a} \exp\left(-\int_{0}^{z} d_{s}(s)ds\right) \exp\left(-\int_{z}^{a} d_{q}(s)ds\right) dz\right) da \right)^{-1}\] \[\times \left(1 - \int_{a_{r}}^{a_{m}} \beta_{s}(a) \exp\left(-\int_{0}^{a} d_{q}(s)ds\right) da\right)^{-1} . (50)\]

Parameter \(R_{vc}\) is the complex model parameter which depends from all vaccine efficacy-dependent coefficients and parameters of the model. The higher the vaccine efficacy value the lower HPV virulence/virus activity \((\alpha)\), the lower HPV replication \((nd_p^*)\) and the higher the HPV death rate \((d_{v0})\) which means the decrease in parameter \(R_{vc}\). Derivative \(\frac{\partial R_0}{\partial R_{vc}} > 0\) for all \(R_{vc} > 0\). If \(w_2(0) > 0\) and \(w_3(0) > 0\) simultaneously and \(R_L < 1\) then the more vaccine efficacy the lower the values of \(R_{vc}\) and \(R_0\), and when \(R_0 < 1\) the DFE is asymptotically stable. We conclude that organism is resistant to HPV infection, remains healthy and can rid of the lesion tissue. In this case we can get the threshold value \(\bar{R}_{vc}\) of \(R_{vc}\):

\[0 < R_{vc} < \bar{R}_{vc} = \left(\int_{a_r}^{a_m} \left(\beta_s(a) \int_0^a \exp\left(-\int_0^z d_s(s)ds\right) \exp\left(-\int_z^a d_q(s)ds\right) dz\right) da\] \[\times \int_0^{a_d} \tilde{d}_p(a) \exp\left(-\int_0^a d_q(s)ds\right) da \left(1 - \int_{a_r}^{a_m} \beta_s(a) \exp\left(-\int_0^a d_q(s)ds\right) da\right)^{-1}\] \[+ \int_0^{a_d} \left(\tilde{d}_p(a) \int_0^a \exp\left(-\int_0^z d_s(s)ds\right) \exp\left(-\int_z^a d_q(s)ds\right) dz\right) da\right)^{-1} \left(\int_0^{a_d} \exp\left(-\int_0^a d_s(s)ds\right) da\right). \tag{51}\]

The criterion for local asymptotic stability of DFE, given by inequality (51), is based on a comparison of rate of HPV infection, HPV replication and HPV mortality (parameter \(R_{vc}\)) with fertility and death (apoptosis) rates of susceptible and infected cells (parameter \(\bar{R}_{vc}\)).

Remark 2. Since \(\frac{\partial R_0}{\partial R_{vc}} > 0\) for all \(R_{vc} > 0\), basic reproduction number of infection \(R_0\) is increasing function of \(N_{S_0}^*\). The DFE value of susceptible cell subpopulation size plays the role of the carrying capacity of susceptible cell environment for HPV subpopulation. If \(N_{S_0}^*\) is low enough, such that \(R_0(N_{S_0}^*) < 1\), carrying capacity of environment is insufficient for the development of HPV subpopulation. It means that starting from initial value \(V_0 > 0\) HPV subpopulation eventually disappears and cell population evolves to the asymptotically stable DFE. If \(N_{S_0}^*\) is large enough, such that \(R_0(N_{S_0}^*) \ge 1\) (case (iii) of Theorem 4.1), carrying capacity of susceptible cell environment is sufficient for development of HPV subpopulation and

DFE is unstable. From equation (51) we can get a threshold value of the carrying capacity of susceptible cell environment \(\hat{N}_{S_0}^*\) for which DFE is locally asymptotically stable:

\[N_{S_0}^* \le \hat{N}_{S_0}^* = \frac{R_{vc} d_{v0}}{\alpha n d_p^*}.\] (52)

Threshold value \(\hat{N}_{S_0}^*\) depends both from the vaccine efficacy-dependent constants \(\alpha\), \(nd_p^*\) and \(d_v\), and from fertility and death (apoptosis) rates of susceptible and infected cell populations. Since the exact value of \(N_{S_0}^*\) cannot be evaluated, we can assume that \(N_{S_0}^*\) depends from initial values \(S_0(a)\) and \(V_0\), and study the dynamical regimes of system defined by criteria (52) numerically.

5. Local asymptotic stability of EE \((S^*(a), I^*(a), L^*(a), V^*)\)

Linearizing System (1) – (6) at the endemic equilibrium (EE) \((S_{-}^*(a), I^*(a), L^*(a), V^*)\) we arrive to the system for perturbations \(\bar{\psi}_s(a,t) = \psi_s(a) \exp(\lambda t)\) for \(S^*(a) \ge 0\), \(\bar{\psi}_i(a,t) = \psi_i(a) \exp(\lambda t)\) for \(I^*(a) \ge 0\), \(\bar{\psi}_i(a,t) = \psi_i(a) \exp(\lambda t)\) for \(L^*(a) \ge 0\) and \(\bar{\psi}_v(t) = \psi_v \exp(\lambda t)\) for \(V^* > 0\):

\[\frac{d\psi_s(a)}{da} = -(\lambda + d_s(a) + \alpha V^*)\psi_s(a) - \alpha S^*(a)\psi_v, \tag{53}\]

\[\frac{d\psi_i(a)}{da} = -(\lambda + d_q(a) + \delta(V^*))\psi_i(a) + (\alpha S^*(a) + \delta'(V^*)I^*(a))\psi_v + \alpha V^*\psi_s(a), \tag{54}\]

\[\frac{\partial \psi_l(a)}{\partial a} = -(\lambda + d_l(a))\psi_l(a) + \delta(V^*)\psi_i(a) + \delta'(V^*)I^*(a)\psi_v,\tag{55}\]

\[0 = -(\lambda + \tilde{d}_v(V^*))\psi_v + nd_p^* \int_0^{a_d} \tilde{d}_p(a)\psi_i(a)da,\] (56)

\[\psi_s(0) = \int_{a_r}^{a_m} \beta_s(a) \psi_s(a) da, \psi_i(0) = \int_{a_r}^{a_m} \beta_s(a) \psi_i(a) da, \psi_l(0) = \int_{a_c}^{a_g} \beta_l(a) \psi_l(a) da.\] (57)

where \(\tilde{d}_v(V^*) = \frac{d_{v0}d_{v1}V^*}{2\sqrt{1+d_{v1}V^*}} + d_{v0}\sqrt{1+d_{v1}V^*} > 0\), \(\delta'(V^*) = \frac{3\delta_d V^{*2}}{(h_d^3+V^{*3})^2} > 0\). Integrating Equation (53) and using (9), (15) yields:

\[\psi_{s}(a) = \psi_{v} \left( D_{S} \exp\left( -\int_{0}^{a} (\lambda + d_{s}(s) + \alpha V^{*}) ds \right) - \alpha A_{S} \int_{0}^{a} \exp\left( -\int_{0}^{z} (d_{s}(s) + \alpha V^{*}) ds \right) \right)\] \[\times \exp\left( -\int_{z}^{a} (\lambda + d_{s}(s) + \alpha V^{*}) ds \right) dz = \psi_{v} \left( D_{S} \exp\left( -\int_{0}^{a} (\lambda + d_{s}(s) + \alpha V^{*}) ds \right) \right)\] \[-\frac{\alpha A_{S}}{\lambda} \left( \exp\left( -\int_{0}^{a} (d_{s}(s) + \alpha V^{*}) ds \right) - \exp\left( -\int_{0}^{a} (\lambda + d_{s}(s) + \alpha V^{*}) ds \right) \right), \quad (58)\]

\[D_{S} = (-\alpha) \left( 1 - \int_{a_{r}}^{a_{m}} \beta_{s}(a) \exp\left( -\int_{0}^{a} (\lambda + d_{s}(s) + \alpha V^{*}) ds \right) da \right)^{-1} \int_{a_{r}}^{a_{m}} \beta_{s}(a) \int_{0}^{a} S^{*}(z)\] \[\times \exp\left( -\int_{z}^{a} (\lambda + d_{s}(s) + \alpha V^{*}) ds \right) dz da = (-\alpha) \left( 1 - \int_{a_{r}}^{a_{m}} \beta_{s}(a) \exp\left( -\int_{0}^{a} (\lambda + d_{s}(s) + \alpha V^{*}) ds \right) da \right)^{-1}\] \[\times \frac{A_{S}}{\lambda} \left( \int_{a_{r}}^{a_{m}} \beta_{s}(a) \left( \exp\left( -\int_{0}^{a} (d_{s}(s) + \alpha V^{*}) ds \right) - \exp\left( -\int_{0}^{a} (\lambda + d_{s}(s) + \alpha V^{*}) ds \right) \right) da \right)\] \[= -\frac{\alpha A_{S}}{\lambda}. \quad (59)\]

Plugging (59) into Equation (58) yields (λ ̸= 0):

\[\psi_s(a) = \psi_v \left( -\frac{\alpha A_S}{\lambda} \right) \exp\left( -\int_0^a (d_s(s) + \alpha V^*) ds \right).\] (60)

Linear Equation (60) is only correct for λ > 0. Integrating Equation (54) yields:

\[\psi_i(a) = D_I \exp\left(-\int_0^a (\lambda + d_q(s) + \delta(V^*))ds\right) + \psi_v \frac{\tilde{d}_v(V^*)}{nd_p^* \left((1 - R_I)^{-1} B_I B_S + C_S\right)} f(a, \lambda),\] (61)

\[f(a,\lambda) = (\alpha A_s)^{-1} \int_0^a (\alpha S^*(y) + \delta'(V^*)I^*(y) + \alpha V^*\psi_s(y))\] \[\times \exp\left(-\int_y^a (\lambda + d_q(s) + \delta(V^*))ds\right) dy = \int_0^a \left(\left(1 - \frac{\alpha V^*}{\lambda}\right) \exp\left(-\int_0^y (d_s(s)ds + \alpha V^*)ds\right)\right)\] \[+ \frac{V^*\delta'(V^*)B_S}{(1 - R_I)} \exp\left(-\int_0^y (d_q(s)ds + \delta(V^*))ds\right) + V^*\delta'(V^*) \int_0^y \exp\left(-\int_0^z (d_s(s)ds + \alpha V^*)ds\right)\] \[\times \exp\left(-\int_z^y (d_q(s)ds + \delta(V^*))ds\right) dz \exp\left(-\int_y^a (\lambda + d_q(s) + \delta(V^*))ds\right) dy, \quad (62)\]

\[D_I = \int_{a_r}^{a_m} \beta_s(a)\psi_i(a)da, \tag{63}\]

Integrating Equation (55) yields:

\[\psi_{l}(a) = D_{L} \exp\left(-\int_{0}^{a} (\lambda + d_{l}(s))ds\right) + \psi_{v} \left(\frac{\tilde{d}_{v}(V^{*})\delta(V^{*})}{nd_{p}^{*}\left((1 - R_{I})^{-1}B_{I}B_{S} + C_{S}\right)} \int_{0}^{a} f(y, \lambda)\right)\] \[\times \exp\left(-\int_{y}^{a} (\lambda + d_{l}(s))ds\right) dy + \delta'(V^{*}) \int_{0}^{a} I^{*}(y) \exp\left(-\int_{y}^{a} (\lambda + d_{l}(s))ds\right) dy\] \[+ D_{I} \left(\delta(V^{*}) \int_{0}^{a} \exp\left(-\int_{0}^{y} (\lambda + d_{q}(s) + \delta(V^{*}))ds\right) \exp\left(-\int_{y}^{a} (\lambda + d_{l}(s))ds\right) dy\right). \quad (64)\]

Multiplying Equation (58) by \(\beta_l(a)\) and integrating it with respect to a from \(a_c\) to \(a_g\), after a little algebra we arrive to the linear equation for \(D_L\) which is defined through \(\psi_v\) and \(D_I\):

\[D_{L}\left(1 - \int_{a_{c}}^{a_{g}} \beta_{l}(a) \exp\left(-\int_{0}^{a} (\lambda + d_{l}(s))ds\right) da\right) = \psi_{v}\left(\frac{\tilde{d}_{v}(V^{*})\delta(V^{*})}{nd_{p}^{*}\left((1 - R_{I})^{-1}B_{I}B_{S} + C_{S}\right)} \int_{a_{c}}^{a_{g}} \beta_{l}(a)\right)\] \[\times \int_{0}^{a} f(y,\lambda) \exp\left(-\int_{y}^{a} (\lambda + d_{l}(s))ds\right) dy da + \delta'(V^{*}) \int_{a_{c}}^{a_{g}} \beta_{l}(a) \int_{0}^{a} I^{*}(y) \exp\left(-\int_{y}^{a} (\lambda + d_{l}(s))ds\right) dy da\right)\] \[+ D_{I}\left(\delta(V^{*}) \int_{a_{c}}^{a_{g}} \beta_{l}(a) \int_{0}^{a} \exp\left(-\int_{0}^{y} (\lambda + d_{q}(s) + \delta(V^{*}))ds\right) \exp\left(-\int_{y}^{a} (\lambda + d_{l}(s))ds\right) dy da\right). \tag{65}\]

Since \(R_L < 1\) (65) is correct for all \(\lambda > 0\). Multiplying (61) by \(\beta_s(a)\) and integrating it with respect to a from 0 to \(a_d\), after a little algebra we arrive to the linear equation for unknown \(\psi_v\) and \(D_I\):

\[\psi_v \left( \frac{\tilde{d}_v(V^*)}{nd_p^* \left( (1 - R_I)^{-1} B_I B_S + C_S \right)} \int_{a_r}^{a_m} \beta_s(a) f(a, \lambda) da \right) - D_I \left( 1 - \int_{a_r}^{a_m} \beta_s(a) \exp\left( - \int_0^a \left( \lambda + d_q(s) + \delta(V^*) \right) ds \right) da \right) = 0. \quad (66)\]

Since \(R_I < 1\) (64) is correct for all \(\lambda > 0\). Multiplying (61) by \(nd_p^* \tilde{d}_p(a)\), integrating it with respect to a from 0 to \(a_d\) and plugging it into Equation (56) we arrive to another linear equation for \(D_I\) and \(\psi_v\):

\[\psi_v \left( \lambda + d_v - \frac{\tilde{d}_v(V^*)}{\left( (1 - R_I)^{-1} B_I B_S + C_S \right)} \int_0^{a_d} \tilde{d}_p(a) f(a, \lambda) da \right) - D_I \left( n d_p^* \int_0^{a_d} \tilde{d}_p(a) \exp\left( - \int_0^a \left( \lambda + d_q(s) + \delta(V^*) \right) ds \right) da \right) = 0. \quad (67)\]

Equating to zero determinant of System (64), (65) after a little algebra we obtain the characteristic equation for \(\lambda\) (\(\lambda \neq 0\)):

\[\left(1 - \int_{a_r}^{a_m} \beta_s(a) \exp\left(-\int_0^a (\lambda + d_q(s) + \delta(V^*)) ds\right) da\right)^{-1}\] \[\times \int_{a_r}^{a_m} \beta_s(a) f(a, \lambda) da \int_0^{a_d} \tilde{d}_p(a) \exp\left(-\int_0^a (\lambda + d_q(s) + \delta(V^*)) ds\right) da\] \[+ \int_0^{a_d} \tilde{d}_p(a) f(a, \lambda) da = \left(\frac{\lambda + \tilde{d}_v(V^*)}{\tilde{d}_v(V^*)}\right) \left((1 - R_I)^{-1} B_S B_I + C_S\right). \quad (68)\]

Characteristic Equation (68) is pretty complex for theoretical analysis of existence of positive roots and will be studied numerically in the next section. We arrive to the theorem.

Theorem 5.1. Let conditions of existence of the EE \((S^*(a), I^*(a), L^*(a), V^*)\) given in Theorem 4.1 hold. If characteristic equation (68) has at least one real positive root \(\lambda > 0\) the EE is unstable, whereas it is locally asymptotically stable every time if Equation (68) does not have real positive root.

Remark 3. In contrast to characteristic equation (45) (DFE) characteristic equation of EE (68) does not depend from parameters of HPV replication n, \(d_p^*\), and depends from parameters of HPV activity (virulence) \(\alpha\) and HPV death rate \(d_v\). So, if vaccine does not sufficiently block HPV replication, DFE is unstable and system moves to the EE which means the formation of lesion tissue. In this case the medicine or therapeutic vaccine can reduce the HPV activity \(\alpha\), increase HPV death rate \(d_v\), decrease the basic reproduction number of lesion subclass \(R_L < 1\) and stabilizes the endemic equilibrium. The further treatment may include the surgical removal of lesion tissue.

6. SIMULATION

Three numerical experiments given in this section illustrate the theoretical results obtained in Theorems 3.1, 4.1, 5.1. We do not aim to cover all possible scenarios driven by the results obtained above, but we are going to model the dynamical regimes which exhibit the correctness of obtained theoretical results and may be interesting in practice. For simulation of the dynamical regimes of autonomous model (1) - (6), we use the explicit formulae of method of characteristics published in works [1], [2].

Experiment I. Locally asymptotically stable DFE. We consider vaccinated and HPV-infected (initial value \(V_0>0\)) organism. Parameters of system satisfy conditions of Theorems 3.1 (case (i)) and 2: \(R_S=1\), \(w_2(0)>0\) and \(w_3(0)>0\) simultaneously, \(R_L<1\), \(R_0<1\) (that is \(R_{vc}<\bar{R}_{vc}\), \(N_{S_0}\leq\hat{N}_{S_0}^*\)). In this case HPV vaccine efficacy is pretty high and HPV environmental carrying capacity is pretty small. It means that organism is resistant toward HPV infection and remains healthy. Results of numerical experiments are shown on Figures 1-4.

The size of susceptible cell population which can be treated as epithelial tissue declines from initial value to the stationary state \(N_{S_0}^* < N_{S_0}\) (Figure 1a) because the basic reproduction number of susceptible cell subpopulation \(R_S = 1\) (9) and initial value of HPV \(V_0 > 0\). The quantity of infected cells (Figure 2), lesion cells (Figure 3) and size of HPV subpopulation (Figure 4) evolve eventually to zero and system approaches the DFE. Numerical results of Experiment I confirm the results obtained in Theorems 3.1 and 4.1.

Experiment II. Unstable DFE. We consider vaccinated and HPV-infected (initial value \(V_0>0\)) organism. Parameters of system satisfy conditions of Theorems 3.1 (case (i)) and Theorem 4.1 (case (iii)): \(R_S=1\), \(w_2(0)>0\) and \(w_3(0)>0\) simultaneously, \(R_L<1\), \(R_0>1\) (that is \(R_{vc}>\bar{R}_{vc}\), \(N_{S_0}>>\hat{N}_{S_0}^*\)). In this case HPV vaccine efficacy is pretty high (the same as considered above in Experiment I) but HPV environmental carrying capacity is pretty large \(N_{S_0}>>\hat{N}_{S_0}^*\). In fact, we use here the same set of coefficients of system as in Experiment I (considered above) and just enlarge the initial value of susceptible cell subclass size \(N_{S_0}>>\hat{N}_{S_0}^*\). It means that we consider the larger epithelial tissue in comparison with previous experiment. Results of numerical experiments are shown on Figures \(S_0>0\). HPV subpopulation gets the larger environmental carrying capacity, virus replicates and destroys the epithelial tissue. Since the basic reproduction number of susceptible cell subpopulation \(R_S=1\) and the basic reproduction number of infection \(R_0>1\), susceptible

cells cannot proliferate in sufficient quantity in time to recover their subpopulation and the size of their sub-population evolves eventually to zero with time (Figure 5). Since the death rate of infected cells is dq(a) ≥ ds(a) and the basic reproduction number of lesion cell subpopulation R^L < 1 the size of infected, lesion cell subpopulations and, as a consequence, size of HPV subpopulation evolve eventually to zero with time (Figures 6 - 8). In this case the DFE is unstable and system approaches the trivial equilibrium.

Experiment III. Asymptotically stable EE. We consider vaccinated and HPV-infected (initial value V⁰ > 0) organism. Parameters of system satisfy conditions of Theorems 3.1 (case (ii)) and Theorem 5.1: R^S > 1, R^I < 1, R^L < 1 and characteristic equation (68) does not have real positive root. Since conditions of

Figure 7: Exp.II. Graph of \(N_L(t)\).

Figure 8: Exp.II. Graph of V(t).

existence of the DFE given in Theorem 3.1 (case (i)) do not hold the DFE does not exist, and if HPV vaccine efficacy is not sufficiently high system approaches the EE. In this case HPV infection causes lesion tissue - dysplasia (precancerous tissue) and cervical cancer. It means that organism is not resistant toward HPV infection and gets sick. Results of numerical experiments are shown on Figures 9-12. The size of susceptible cell population (epithelial tissue) declines from initial value to the stationary state \(N_{S_0}^* < N_{S_0}\) (Figure 9) that is a result of balance between the susceptible cells proliferation \((R_S > 1)\) and their infection by HPV. After a transient process the quantity of infected (Figure 10) and lesion (Figure 11) cells and size of HPV subpopulation (Figure 12) evolve eventually to positive equilibrium, and system approaches the EE. Numerical results of Experiment III confirm results obtained in Theorems 3.1 and 5.1.

Figure 9: Exp.III. Graph of \(N_S(t)\).

Figure 10: Exp.III. Graph of \(N_I(t)\).

Figure 11: Exp.III. Graph of \(N_L(t)\).

Figure 12: Exp.III. Graph of V(t).

7. CONCLUSION

In this paper we studied the epidemic model of age-structured sub-populations of susceptible, infectious, lesion cells and unstructured population of human papilloma virus (HPV) (SILV epidemic model). The model considers HPV transmission and epithelial tissue cells (susceptible, infected and lesion cells) population dynamics of vaccinated organism. Stability analysis of autonomous system revealed the conditions of existence of unique disease-free and unique endemic equilibria (DFE and EE, respectively) which do not contain efficacy vaccine-dependent coefficients and parameters of model. In particular, condition of existence of DFE depends on the basic reproduction number of susceptible cells \(R_S\) which contains only fertility and death (apoptosis) rates of susceptible cells. Condition of existence of EE depends on the basic reproduction numbers of susceptible, infected and lesion cells \(R_S\), \(R_I\) and \(R_L\), respectively, which depend on the fertility and death (apoptosis) rates of susceptible, infected and lesion cells and the progression rate from infected to lesion cells \(\delta\). The main indicator of DFE asymptotic stability – basic reproduction number of infection \(R_0\) is increasing function of \(R_{vc} = \frac{\alpha N_{S_0}^* n d_p^*}{d_{v0}}\). Parameter \(R_{vc}\) contains the vaccine efficacy-dependent constants and DFE value of susceptible population quantity \(N_{S_0}^* = \int\limits_{a}^{a_d} S_0^*(a) da\). The higher the vaccine efficacy the lower HPV virulence/virus activity \((\alpha)\), the lower HPV replication \((nd_p^*)\) and the higher the HPV death rate \((d_{v0})\) which means the decrease in parameters \(R_{vc}\) and \(R_0\). If coefficients of the system satisfy conditions of Theorem 4.1 and \(R_0 < 1\) the DFE is asymptotically stable that is vaccine efficacy is sufficiently high, vaccinated organism is resistant to HPV infection, and can rid of the lesion tissue.

On the other hand, the basic reproduction number of infection \(R_0\) (50) is increasing function of \(N_{S_0}^*\). Since we can't evaluate the exact value of \(N_{S_0}^*\), we can just assume that \(N_{S_0}^*\) is proportional to the initial values \(N_{S_0}=\int\limits_0^{a_d}S_0(a)da\) and \(V_0\), and study the asymptotic stability of DFE with different initial values numerically. In Experiments I and II, two systems differed from each other only in their initial values \(N_{S_0}\)(Figures 1 - 4 and 5 - 8) and have the same vaccine efficacy-dependent parameters. In the first case system with smaller initial value \(N_{S_0}\) and smaller \(R_{vc}\) (for which \(R_0 < 1\)) satisfied conditions of asymptotic stability of DFE of Theorem 4.1 and approached the asymptotically stable DFE (Figures 1 - 4). Quantity of infected, lesion cells and HPV approached zero, vaccinated organism was resistant to HPV infection and got rid of the lesion tissue and HPV. In the second case system with bigger initial value \(N_{S_0}\) and bigger \(R_{vc}\) (for which \(R_0 > 1\)) did not satisfy conditions of asymptotic stability of DFE of Theorem 4.1 and approached the trivial equilibrium (Figures 5 - 8). Quantity of all susceptible, infected, lesion cells and HPV approached zero, vaccinated organism was not resistant to HPV infection and infection destroyed epithelial tissue layer. We considered here a particular type of unstable DFE when system approached the trivial equilibrium with small basic reproduction number of lesion cells population \(R_L < 1\), because we used the same coefficients of system as in the first case (satisfying conditions of Theorem 4.1) except only the initial value. We can conclude that asymptotic stability of DFE (i.e. resistance of vaccinated organism to HPV infection) depends from the initial quantity of susceptible cells of epithelial tissue layer. The tissue of smaller size can be more persistent to HPV infection than the tissue of larger size. Thus, HPV vaccine efficacy should always be sufficiently high to reduce \(R_{vc}\) to a value where \(R_0 < 1\) for the maxi-mum possible epithelial tissue size. Thus, the analysis of asymptotical stability of DFE of SILV epidemic model helps to better understand the features of HPV vaccination and develop recommendations for practical medicine when studying the efficacy of new vaccines in control groups in clinical trials [15], [19]. In particular, we can conclude that HPV vaccine efficacy should always be sufficiently high to reduce \(R_{vc}\) to a value where \(R_0 < 1\) for the maximum possible

The local asymptotic stability of EE means that vaccinated organism is not resistant to HPV infection which leads to viral disease and formation of lesion tissue. The indicator of EE asymptotic stability does not depend from vaccine efficacy-dependent set of parameters and depends on the fertility and death (apoptosis) rates of susceptible, infected and lesion cells, the progression rate from infected to lesion cells, HPV death rate and its partial derivative. Asymptotically stable EE studied in the Experiment III, where system satisfied conditions of asymptotic stability of EE of Theorem 5.1. Quantity of susceptible, infected, lesion cells and HPV approached positive equilibrium (Figures 10 - 12), which means that vaccinated organism is not resistant to HPV infection and gets sick, which leads to the formation of lesion tissue.

Overall, we can conclude that stability analysis of age-structured SILV epidemic model with vaccination

provides the theoretical instrument for study the dynamics of susceptible, infected, lesion cells and HPV populations that help us better understand the impact of vaccination on HPV infectious disease.

ACKNOWLEDGEMENT

I would like to acknowledge Faculty of Sciences, University of Manitoba, Winnipeg, Canada, for the organizational support of this study.