Stability analysis of an age-structured model of cervical cancer cells and HPV dynamics

Stability analysis of an autonomous epidemic model of an age-structured sub-populations of susceptible, infected, precancerous and cancer cells and unstructured sub-population of human papilloma virus (HPV) (SIPCV epidemic model) aims to gain an insight into the features of cervical cancer disease. The model considers the immune functional response of organism to the virus population growing by the HPV-density dependent death rate, while the death rates of infected, precancerous and cancerous cells do not depend on the HPV quantity because the immune system of organism does not respond to its own cells. Interaction between susceptible cells and HPV is described by the LotkaVoltera incidence rate and leads to the growth of infected cells. Some of infected cells become precancerous cells, and the other apoptosis when viruses leave infected cells and are ready to infect new susceptible cells. Precancerous cells partially become cancer cells with the density-dependent saturated rate. Conditions of existence of the endemic equilibrium of system were obtained. It was proved that this equilibrium is always locally asymptotically stable whenever it exists. We obtained: (i) the conditions of cancer tumor localization (asymptotically stable dynamical regimes), (ii) outbreak of cancer cell population (that may correspond to metastasis), (iii) outbreak of dysplasia (precancerous cells) which induces the outbreak of cancer cells (that may correspond to metastasis). In cases (ii), (iii) the conditions of existence of endemic equilibrium do not hold. All cases are illustrated by numerical experiments.


Introduction
Cervical cancer induced by sexually transmitted human papilloma virus (HPV) is the second most common cancer in women worldwide [1][2][3][4]. The widespread prevalence of this disease stimulated extensive studies of HPV-epidemic models over the last several years. The research in this field is carried out in two directions: (i) epidemiology and population science (models of spread of HPV through a population, HPV vaccination of people, etc.) [5,6]; (ii) epidemiology and biological sciences (cellular and tissue modelling, cell-HPV interaction studies, etc.) [7][8][9]. The early studies in the framework of the second scientific direction (ii) used the unstructured models of population dynamics for several compartments -subpopulations of biological cells and HPV. The dynamics of virus population can be efficiently described by the unstructured model on the basis of nonlinear ODE because viruses are non-living things which do not proliferate and can replicate only within a living host cell. But unstructured models are not suitable for the modelling of cell population dynamics because they do not describe the cell's life history and provide only a restricted description of their population dynamics in many applications. The urgent need to gain insight the complex biological processes for deep and accurate understanding of patterns of population dynamics and a variety of dynamical regimes of populations motivated the development and implementation of age-structured, and more generally, physiologically structured, cell population dynamics and tumor growth modelling [3,[10][11][12][13][14][15][16][17][18][19][20][21][22][23][24].
In this paper we study the new age-structured model of susceptible, infected, precancerous, cancer cells populations and unstructured model of human papilloma virus population (SIPCV epidemic model) dynamics which is continuation of the previous research described in works [8,25]. In this model we use the L. Hayflick limit theory [26] for modelling the proliferation in cell subpopulations. The model describes the life history of each cell (cell aging by L. Hayflick): birth, maturing up to the age when they can proliferate, division a limited number of times at the reproductive age, aging up to the final reproductive age and death. The cell division and mortality in subpopulations are described by the birth and death rates, respectively. The death rates of infected, precancerous and cancer cell subpopulations in our model do not depend on the HPV quantity since the immune system of organism does not respond to its own cells [1,3,4]. The death rate of HPV is considered as density-dependent function since immune system responds to the virus population growth [1,3,4]. Interaction between susceptible cells and HPV is described by the Lotka-Voltera incidence rate and leads to the growth of infected cells [25]. Some of infected cells become precancerous cells, and the other apoptosis when viruses leave infected cells and are ready to infect new susceptible cells [8,28]. Precancerous cells partially become cancer cells with the density-dependent saturated rate [8]. We assume that cancer cells do not apply pressure on the tissues of organism and have no effect on the proliferation and mortality of other cells and replication of HPV. This model is studied both theoretically and numerically. The existence theorem and explicit recurrent formula for the solution of the age-structured SIPCV model (like in work [25]) are beyond the aim and scope of this paper due to the complexity of the model and will be the subject of our further study.
The stability analysis of model is based on the study of conditions of existence of the positive (endemic) equilibrium of system and its asymptotical stability ( [13,[29][30][31][32][33][34][35][36][37]). We obtain the restrictions on the basic reproduction numbers of susceptible and cancer cell subpopulations, and coefficients of the system which guarantee existence of the endemic equilibrium. We prove that this equilibrium is always locally asymptotically stable whenever it exists. The numerical experiments confirm theoretical results and provide us with several topologically non-equivalent phase portraits of the model. Simulations reveal two asymptotically stable regimes with non-oscillating and oscillating dynamics in the vicinity of positive equilibrium and, at least, two unstable dynamical regimes (when the positive equilibrium does not exist). From the biological point of view the asymptotically stable dynamical regimes of cell-HPV population mean the localization of cancer without spreading it in organism. The unstable dynamical regimes of cell-HPV population mean the cancer metastasis. Overall, the theoretical and numerical analysis of autonomous age-structured SIPCV epidemic model help us better understand the features of HPV infectious and cancer disease.

Model
SIPCV epidemic model considers the biological tissue which consists of susceptible (noninfected), infected (without significant changing of morphology, CIN I and CIN II stages [1][2][3][4]27]), precancerous (with changed by virus morphology -dysplasia, but is differentiable yet, CIN III stage [2,3,27]), cancer (nondifferentiable) cells and human papilloma virus (HPV) that moves freely between cells. The agespecific densities of susceptible, infected, precancerous and cancer cells subclasses are denoted as , respectively. The dynamics of cell subpopulations is described by the autonomous nonlinear age-structured model with death rates of susceptible cells , infected cells , precancerous cells and cancer cells , with an age reproductive window of non-cancer cells [ , ] and cancer cells [ , ], < , < , (the age reproductive window of cancer cells is shifted in relation to the age reproductive window of non-cancer cells due to the abnormal program of cancer cell division when they divide before reaching the maturity of the reproductive window of non-cancer cells and therefore become nondifferentiable cells), the same birth rate of susceptible, infected and precancerous cells 0 , and the birth rate of cancer cells . Due to the adaptive behaviour of the HPV immune system of organism (both T-killers cells and humoral immunity) does not respond to the infected, precancerous and cancer cells, that is their death rates do not depend on the HPV quantity. Since viruses are non-living things which do not proliferate and replicate only within a living host cell, the dynamics of HPV quantity ( ) is described by the non-linear ODE with density-dependent death rate ( ). The latter describes the immune response of organism to the HPV population growth. The interaction between susceptible cells and HPV is a product of the Lotka-Voltera incidence rate ( ) ( , ) (where is an infection rate) and result in the growth of infected cells. Infected cells partially move to the precancerous subpopulation with rate ( , ) (where is a progression rate from infected to precancerous cells (dysplasia)) and partially apoptosis with rate ( ), when viruses leave infected cells and are ready to infect susceptible cells (where is a mean number of virions produced by an infectious cell, the death rate of infected cells . When the abundance of precancerous cells (dysplasia) increases from the small value, the risk of developing cancer cells increases from the small value too. In this case ( ) is directly proportional to the with coefficient (progression rate). Since ( ) is a fraction of precancerous cells which move to cancer subclass per unit of time, it is a bounded parameter which increases and eventually tends to the saturated level / < 1 with → ∞, where is a coefficient of saturation [8].
These assumptions lead to the following autonomous epidemic model where ( ) is an initial density of susceptible cells, 0 is an initial value of HPV quantity. We impose the following restrictions on the density-dependent HPV death rate and cell's birth and death rates [1][2][3][4]27,28]: , , , > 0, Equations (2.11) and (2.12) consider the positiveness and boundness of cell birth and death rates. The restrictions (2.11) mean that increasing of HPV quantity changes the characteristics of intracellular space that result in the organism immune response through the activation of cell immunity (T-killers) and humoral immunity (B-lymphocytes) that leads to the elimination of viruses (i.e. monotone increasing of their death rate). The constant 0 restricts the intrinsic density independent part of death rate.

Existence of the positive (endemic) equilibrium of autonomous system (2.1)-(2.10)
The trivial, disease free equilibrium (DFE) of the system (2.1)-(2.10) exists if an infection rate = 0. In this case biological tissue consists only from susceptible cells and its dynamics is defined by the basic reproduction number of susceptible cell subpopulation. Since we consider only positive infection rate > 0 (Eq (2.12)), the analysis of DFE is beyond of our study. where we assume that < ( * ) * . (2.21) Using Eq (2.18) we obtain the integral equation for the solution of Eq (2.13): * ( ) = 0 ∫ * ( ) (−( + * ) ).

(2.22)
Integrating Eq (2.22) with respect to form to , after a little algebra we arrive to the transcendental equation for determination of * : Integrating Eq (2.22) with respect to form 0 to we obtain the equation for the equilibrium quantity of susceptible cells: ).

(2.24)
Substituting this expression in Eq (2.22) we obtain the solution * ( ) and equation for the equilibrium density of susceptible cells * : * ( ) = * (−( + * ) ), where is an auxiliary positive constant. Using the boundary condition (2.19) we can write the integral equation for the solution of Eq (2.14): Integrating Eq (2.27) with respect to form to , and using Eq (2.25), after a little algebra we arrive to the integral equation: Expressing ∫ * ( ) from Eq (2.28) and substituting it together with Eq (2.23) in Eq (2.27), after a little algebra we obtain the positive equilibrium * ( ): * ( ) = (−( + * ) ), HPV equilibrium * has to satisfy restriction * < (̃− ) −1 , (2.31) which guarantees positiveness of the auxiliary constant and equilibrium * ( ). Integrating * ( ) (Eq (2.29)) with respect to form 0 to and using Eq (2.30), we have * = * * Integrating Eq (2.34) with respect to form to , and using Eqs (2.23), (2.29), (2.30), after a little algebra we arrive to the integral equation: Since the left-hand side of Eq (2.35) must be positive the equilibrium quantity * and coefficients of the system have to satisfy the restriction:  Using the simple algebra, it is easy to verify that the equilibrium quantity of precancerous cells (2.41) is positive * > 0 and satisfies restriction (2.36) if the HPV equilibrium * and coefficients of the system satisfy condition * < −1 ( − + −1 ).
Integrating Eq (2.43) with respect to form to , after a little algebra we arrive to the equation: Since the expression in the right-hand side of Eq (2.44) is positive, integral ∫ * ( ) > 0 , if Integrating Eq (2.42) with respect to form 0 to we obtain the corresponding equilibrium quantity of cancer cells: (2.48) The results obtained above are summarized in the Theorem 2.1. Hence, the necessary and sufficient condition of existence of a positive root * > 0 of transcendental equation ( * ) = 1 (Eq (2.23)) is the condition (0) > 1, that is the basic reproduction number of susceptible cell population is bigger than one:  Integrating Eq (2.59) with respect to from to , and using Eqs (2.23), (2.63), we arrive to the equation: where +̃− − * > 0 for > 0 (see Eq (2.31)). We assume that ≠ * , because when = * function perturbations = 0 (Eq (2.64)), = 0 (Eq (2.63)) and we obtain the trivial solution.
Substituting ∫ ( ) in Eq (2.59) after integration and a little algebra we arrive to the solution of Eq (2.59):        In this case = 0 and we only obtain a trivial solution of system that is the characteristic number cannot take the trivial value.
Finally, if condition (2.84) does not hold and Eq (2.79) does not have a real positive root, we analyse existence of the complex characteristic roots with nonnegative real part. We seek the pure imaginary root = (where ≠ 0 is unknown real parameter) for which there exists the non-trivial solution of system (2.51)-(2.57). Substituting = in Eqs (2.51)-(2.55) and separating imaginary parts of equations yields:

Numerical experiments
The theoretical results obtained in Theorems 2.1 and 2.2 are illustrated by the numerical experiments. For simulation of the autonomous model (2.1-2.10), we use the explicit formulae of method of characteristics from [25]. The set of coefficients and initial values of autonomous system (2.1-2.10) used in simulations is given in Appendix A.
First, using the bisection method we find numerically the positive real root of Eq ( Since in case (iv) conditions (2.36), (2.42) of Theorem 2.1 do not hold, the positive equilibrium of system (2.1)-(2.10) does not exist. Graphs in Figures 3a, 3b, 3e show the oscillating dynamics of the quantity of susceptible, infected cell subpopulations and HPV subpopulation with following their convergence to some stationary values, like in the previous case (iii).
Since ′ ( ) > 0 in Figure 3c  In the last group of experiments, we study the behavior of solution when the basic reproduction number of cancer cell population < 1 and > 1. We consider the same set of the coefficients and constants as in the Experiment II but take the bigger values of birth rates of cancer cells (Appendix A) that corresponds to the bigger value of . In the first simulation = 0,99993, and the set of all coefficients and constants (except ) is the same as in case (iii). Since  Figure 2ac, e. The dynamical regime of the cancer cell subpopulation obtained in this experiment differs from the one obtained in case (iii) and is the same as in simulation of case (iv) (for < 1). The results of simulation are shown in the phase diagram in Figure 4b. Since ′ ( ) > 0, we observe the wiggle dynamical regime of the cancer cell subpopulation quantity (without oscillation) with subsequent unlimited exponential growth. Thus, the absence of the positive equilibrium in this case means the unstable dynamical regime of the system (2.1)-(2.10).

Discussion and conclusions
In this paper we study an autonomous epidemic model of age-structured population dynamics of susceptible, infected, precancerous and cancer cells and unstructured model of population dynamics of human papilloma virus (HPV) (SIPCV epidemic model). The model considers the problem of HPV propagation and cancer disease dynamics on the tissue level and includes the competitive system of initial-boundary value problems for semi-linear transport equations with non-local boundary conditions (renewal equations) and initial problem for nonlinear ODE. We carried out the stability analysis of this autonomous system, obtained the conditions of existence of the positive (endemic) equilibrium and proved that this equilibrium is always locally asymptotically stable whenever it exists. Theoretical analysis revealed two key parameters important for the study of cervical cancer diseasethe basic reproduction numbers of susceptible cell population and cancer cell population . The necessary and sufficient condition of existence of the positive equilibrium of system, as it was expected, imposes the restriction on the basic reproduction number of susceptible cell population: > 1, that is the healthy biological tissue must be growing for providing the sufficient environment for the HPV replication and propagation, and development of the cancer tissue. The particularity of the semi-linear transport equation of cancer cell population dynamics is that it does not impact on the dynamics of the other cell and HPV sub-populations since all other equations of the system do not depend from the cancer cell density or cancer cell quantity. This property of model is a consequence of our hypothesizes that immune system of organism is tolerant with respect to its own cervical cancer cells and cancer cells do not apply pressure on the tissues of organism and have no effect on the proliferation and mortality of the other cells and replication of HPV. Thus, when the positive equilibrium of system exists and the cancer cell population is not empty, their farther dynamics depends only on the basic reproduction number . If < 1 cancer cell population evolve eventually to the stationary state that means the localization of tumor tissue, and if > 1 cancer cell population grows infinitely that means the formation of metastasis in organism. The unlimited growth of cancer cell population can also appear even for < 1, when the positive equilibrium of system does not exist, but susceptible, infected cell subpopulations and HPV subpopulation eventually evolve to some stationary values and precancerous cell subpopulation (dysplasia) grows infinitely. In this case precancerous cells induce the cancer cell outbreak and subpopulation of cancer cells grows infinitely. Numerical experiments illustrate and confirm the theoretical results obtained in paper. When the basic reproduction number of susceptible cell subpopulation ≤ 1 (case (i)) the positive equilibrium of system does not exist, and simulation showed that all cell and HPV subpopulations eventually evolve to zero. This is a trivial solution of the model which is not valuable in epidemiological problem. When > 1, < 1, and all coefficients and constants of system (2.1)-(2.10) satisfy conditions (2.21), (2.31), (2.36), (2.42) of Theorem 2.1, the positive equilibrium of system exists. Simulation showed that in this case there are three types of dynamical asymptotically stable regimes: non-oscillating convergence of solution to the positive equilibrium (case (ii)), oscillating convergence of solution to the positive equilibrium (case (iii)), wiggle dynamics with subsequent exponential convergence of solution to the positive equilibrium (first simulation in case (v), = 0.99995). The common feature of all these dynamical regimes is that cancer cell population cannot grow infinitely, their dynamics is related to the dynamics of all other cells of tissue and the quantity of cancer cells eventually evolves to the stationary state. Results of the simulations in these cases confirm the theoretical conclusion about a localization of cancer tissue. When > 1, < 1, and coefficients and constants of system (2.1)-(2.10) do not satisfy condition (2.42) of Theorem 2.1, the positive equilibrium of system does not exist. Simulation showed that in this case the quantity of susceptible, infected cell subpopulations and HPV subpopulation oscillated and eventually evolved to some stationary values and the quantity of precancerous cell subpopulation (dysplasia) showed the wiggle dynamics with subsequent unlimited growth. In spite of the small value of the basic reproduction number < 1, cancer cell subpopulation also grew infinitely together with precancerous cell subpopulation that means the growth of dysplasia and cancer metastasis in organism.
And the last, but not least case when > 1, > 1, and all coefficients and constants of the system satisfy the conditions of Theorem 2.1, except . In this case (the second simulation in case (iv), = 1.00008), like in case (i), system does not have the positive equilibrium, but the dynamical regime of population is significantly different from case (i). The quantities of susceptible, infected, precancerous cells and HPV converge to the stationary states, except the cancer cell quantity. Cancer cell subpopulation exhibits the wiggle dynamics (without oscillation) with consequence unlimited exponential growth, that is the dynamics of cancer cells is not stable and leads to the formation of metastases in organism. Thus, the results of all numerical experiments confirm the conclusions of theoretical analysis.
Overall, the main result of this paper -development of the SIPCV age-structured epidemic model and stability analysis of autonomous system of this model, provide the theoretical instrument for the qualitative analysis of dynamical regimes of susceptible, infected, precancerous, cancerous cells and HPV populations that help us better understand the features of HPV infectious and cancer disease. In the different numerical experiments, in addition to the constants from the Table A1, we use the coefficients given in the Table A2.