Mathematical analysis of a SIPC age-structured model of cervical cancer

Eminugroho Ratna Sari; Fajar Adi-Kusumo; Lina Aryati; Eminugroho Ratna Sari; Fajar Adi-Kusumo; Lina Aryati

doi:10.3934/mbe.2022281

Mathematical Biosciences and Engineering

2022, Volume 19, Issue 6: 6013-6039. doi: 10.3934/mbe.2022281

Previous Article Next Article

Research article Special Issues

Mathematical analysis of a SIPC age-structured model of cervical cancer

1.
Department of Mathematics, Universitas Gadjah Mada, Yogyakarta 55281, Indonesia
2.
Department of Mathematics Education, Universitas Negeri Yogyakarta, Yogyakarta 55281, Indonesia

Academic Editor: Sandra Saliani

Received: 22 February 2022 Revised: 26 March 2022 Accepted: 31 March 2022 Published: 12 April 2022

Human Papillomavirus (HPV), which is the main causal factor of cervical cancer, infects normal cervical cells on the specific cell's age interval, i.e., between the $G_1$ to $S$ phase of cell cycle. Hence, the spread of the viruses in cervical tissue not only depends on the time, but also the cell age. By this fact, we introduce a new model that shows the spread of HPV infections on the cervical tissue by considering the age of cells and the time. The model is a four dimensional system of the first order partial differential equations with time and age independent variables, where the cells population is divided into four sub-populations, i.e., susceptible cells, infected cells by HPV, precancerous cells, and cancer cells. There are two types of the steady state solution of the system, i.e., disease-free and cancerous steady state solutions, where the stability is determined by using Fatou's lemma and solving some integral equations. In this case, we use a non-standard method to calculate the basic reproduction number of the system. Lastly, we use numerical simulations to show the dynamics of the age-structured system.

Keywords:

Citation: Eminugroho Ratna Sari, Fajar Adi-Kusumo, Lina Aryati. Mathematical analysis of a SIPC age-structured model of cervical cancer[J]. Mathematical Biosciences and Engineering, 2022, 19(6): 6013-6039. doi: 10.3934/mbe.2022281

Related Papers:

[1]	Vitalii V. Akimenko, Fajar Adi-Kusumo . Stability analysis of an age-structured model of cervical cancer cells and HPV dynamics. Mathematical Biosciences and Engineering, 2021, 18(5): 6155-6177. doi: 10.3934/mbe.2021308
[2]	Shaoli Wang, Jianhong Wu, Libin Rong . A note on the global properties of an age-structured viral dynamic model with multiple target cell populations. Mathematical Biosciences and Engineering, 2017, 14(3): 805-820. doi: 10.3934/mbe.2017044
[3]	Yicang Zhou, Zhien Ma . Global stability of a class of discrete age-structured SIS models with immigration. Mathematical Biosciences and Engineering, 2009, 6(2): 409-425. doi: 10.3934/mbe.2009.6.409
[4]	Jinliang Wang, Xiu Dong . Analysis of an HIV infection model incorporating latency age and infection age. Mathematical Biosciences and Engineering, 2018, 15(3): 569-594. doi: 10.3934/mbe.2018026
[5]	Lu Gao, Yuanshun Tan, Jin Yang, Changcheng Xiang . Dynamic analysis of an age structure model for oncolytic virus therapy. Mathematical Biosciences and Engineering, 2023, 20(2): 3301-3323. doi: 10.3934/mbe.2023155
[6]	Xue-Zhi Li, Ji-Xuan Liu, Maia Martcheva . An age-structured two-strain epidemic model with super-infection. Mathematical Biosciences and Engineering, 2010, 7(1): 123-147. doi: 10.3934/mbe.2010.7.123
[7]	Xia Wang, Yuming Chen . An age-structured vector-borne disease model with horizontal transmission in the host. Mathematical Biosciences and Engineering, 2018, 15(5): 1099-1116. doi: 10.3934/mbe.2018049
[8]	Jinhu Xu . Dynamic analysis of a cytokine-enhanced viral infection model with infection age. Mathematical Biosciences and Engineering, 2023, 20(5): 8666-8684. doi: 10.3934/mbe.2023380
[9]	Tsuyoshi Kajiwara, Toru Sasaki, Yoji Otani . Global stability of an age-structured infection model in vivo with two compartments and two routes. Mathematical Biosciences and Engineering, 2022, 19(11): 11047-11070. doi: 10.3934/mbe.2022515
[10]	Zhiping Liu, Zhen Jin, Junyuan Yang, Juan Zhang . The backward bifurcation of an age-structured cholera transmission model with saturation incidence. Mathematical Biosciences and Engineering, 2022, 19(12): 12427-12447. doi: 10.3934/mbe.2022580

Abstract

1. Introduction

Cervical cancer is one of the malignant cancer that mostly caused by the Human Papillomavirus (HPV) infection. In 2020, there were more than 604,000 new cases of cervical cancer with 342,000 deaths worldwide, see ^[1,2]. HPV16 and HPV18 are the sub-types of HPV that have higher risk as the causal factor of cervical cancer ^[3]. The cancer occurs in cervical epithelial tissue that consists of three zones, i.e., ectocervix, endocervix, and the transition zone between ectocervix and endocervix, and three layers, i.e., basal, middle, and suprabasal ^[4]. In early stages of HPV infection, the non-invasive lesions of abnormal cervical epithelial cells are found. If not treated properly, these abnormal cells have the ability to develop into cervical cancer. This process is known as Cervical Intraepithelial Neoplasia (CIN) ^[5].

There are three steps that CIN must pass regarding to the development of cervical cancer ^[4]. The first step is when the abnormal lesion occurs in one-third of the epithelial tissue. The average time taken from the first infection is three years. In this step, $90\%$ of the infected cells can be cleared while the $10\%$ is still abnormal. The second step is when the abnormal lesion occurs in two-thirds of the epithelial tissue. In the last step, the abnormalities occur in almost all epithelial tissue. The time required from the first to the third step of CIN is between three to ten years. Without appropriate treatment, HPV infection becomes persistent between five to ten years and then becomes invasive.

The virus starts to enter and to infect the basal epithelial cells when the cells are at a certain maturity. The maturity of the cell is closely related to the cell's age, which is divided into four phases, i.e., ${G} _ {1}$ , $S$ phase, ${G} _ {2}$ and $M$ phase (mitosis). The genomes' replication of HPV highly depend on the stage when the cell age is in ${G} _ {1}$ phase towards the $S$ phase. Proteins from HPV, i.e., E6 and E7, play a role in inhibiting the cell from entering the ${G}_{1}$ phase ^[6].

The degeneration of the cancer malignancy are mainly caused by E6 and E7 oncoproteins. The E6 oncoprotein interacts and inactivates the p53 protein, which plays a role as a tumor gene suppressor, that works in the ${G}_{1}$ phase. The inactivity of p53 will stop the cell cycle through its resistance to the complex works that stimulates the cell cycle to enter the next phase. As a result, the cell will continue to enter the $S$ phase without DNA correction. This abnormal cell has ability to proliferate uncontrollably. Therefore, we assume that the age of abnormality will not start from zero. Persistent infections by high-risk HPV strains trigger continued release of HPV virions. This condition leads to the growth of precancerous lesions, which are at risk of becoming cancer cells.

The studies of the spread of cervical cancer at cellular level have been done in ^{[7,8,9,10,11]}, where the authors consider the progression of cervical cancer that involves the interaction between susceptible cells, infected cells, precancerous cells, and cancer cells sub-populations, and free viruses population. The population of free viruses represents the presence of the HPV in the cervical epithelial tissue. The HPV-infected cells are still active in the cell cycle, but their behaviour is changed. The essential parameters that stimulate the number of cancer cells and their metastasis behaviour of cervical cancer were studied in ^[7]. The existence conditions of the unstable equilibrium point and the boundary in the parameter space that the unstable equilibrium point exists have been done in ^[7] and ^[8]. Moreover, the analytical study of the global stability of the disease-free equilibrium point has been done in ^[9]. In this case, the authors found that the therapy given to the precancerous cells can prolong the patient's life expectancy.

The age structured models of cervical cancer have been studied for the HPV transmission between individuals in human populations, see ^[12,13,14]. In cellular level, one of the important characteristic of human cancer is that the small disturbances on the interactions between quiescent cells, proliferating cells, death cells, and the other types of cell in the cell life spans can trigger to uncontrolled cell division. Therefore, the authors in ^[15,16,17] developed two age-structured model for the proliferating and quiescent cell compartments. The model is important to construct a precisely targeted therapy that includes the differentiated cells population by their position within the cell division cycle ^[18]. The other mathematical model that studied the different transition rules which alter the phase distribution of the cells where the focus phases were $G_1$ and $S$ , was done in ^[19].

An age-structured model of HPV infection at the cellular level that includes the interaction between cells population and free virus population has been introduced by ^[10]. The model was age and time-dependent for susceptible cells, infected cells, precancerous cells, cancer cells, and time-dependent of the free virus population. The author in ^[10] was assumed that the virus interacts with susceptible cells at a constant infection rate.

Due to the fact that since the first infection, the virus will be reproduced and spread in the cervical tissue, see ^[20], the appearance of the free virus population as an independent compartment that depends on the time should be reevaluated. The differences with the system in ^[10] is as follows. We remove the free virus compartment and introduce a new parameter that represents the infection rate of the free viruses called force of infection. Since HPV infects normal cervical cells on the specific cell's age interval, i.e., between the $G_1$ to $S$ phase of the cell cycle ^[6,21], the forco of infection parameter is assumed to be age-dependent. Based on the fact that the duration of a cell infection from susceptible, infected, precancerous, and then become cancerous is comparable with the duration of the cell cycle, we add the natural growth rate of the susceptible cells population to our model, which has not been studied in ^[10]. Since the genomes' replication of HPV highly depends on a certain stage of the cell cycle, we assume that the age of the abnormality of the cell does not start from zero.

Various methods were used to determine the basic reproduction number in the recent theoretical age-structured model. The authors in ^[22,23] derived the reproduction number from the biological meanings of the model parameters. The threshold number also can be defined when calculating the endemic steady state (like the work ^[10,24]). In ^[25], determination of basic reproduction number was obtained by transforming the system into Cauchy problem. In ^[12,13,14], the reproduction number was acquired by analyzing characteristic equation at the disease-free steady state. In this paper, we prove that the roots of the characteristic equation at the disease-free steady state will be positive or negative under certain conditions that fit the properties of a basic reproduction number.

The content of this work is organized as follows. In Section 2, a new age-structured mathematical model of cervical cancer involving the force of infection will be introduced. Determination of the steady state solution and detailed steps to obtain basic reproduction number are carried out in Section 3. In Section 4, we focus our study on the local and global stability of the disease-free steady state solution. We also present the existing conditions and the local stability of the cancerous steady state solution. The study is important to understand some patterns to provide a successful targeted therapy from the perspective of mathematics. In Section 5, we will use numerical simulations to show the behavior of the system. Finally, concluding remarks are given in the last section.

2. Model formulation

A cervical cancer model in tissue level with the interactions between the cells and the virus populations was introduced in ^[7]. However, the model is an ODE system that depend only on the time and did not not involve the age dependency on the HPV transmission.

By the fact that the infection of HPV depends on the cell age, see ^[6], we add a new independent variable that represents the age of the cells ( $a$ ), and a new age-dependent parameter called force of infection that replaces the role of the virus compartments in ^[7] in our system. The term "cell age" (in hours) is defined as the time spent by a cell to complete a cell cycle. The cell cycle starts from a new cell as a result of cell division until the next cell division occurs. We denote the newborn cell as $a = 0$ and the maximum age of the cell as $a_\sigma$ .

Let $a_{1}$ is the age of the cell that has the ability to be infected by the virus or become abnormal. Since the first infection, the virus can be reproduced and spread in the tissue, so that the role of the free virus in ^[7] and ^[10] can be integrated in the transmission rate parameter. In our model, the density of the cells population in a tissue is divided into four sub-populations that depend on age $a$ and time $t$ , i.e., susceptible cells, infected cells, precancerous cells, and cancer cells that can be denoted by $S\left(a, t\right)$ , $I\left(a, t\right)$ , $P\left(a, t\right)$ , and $C\left(a, t\right)$ , respectively. We assume that the death rate of cells with age $a$ , i.e., $\mu\left(a\right)$ , of each kind of cell is the same. The total density of the cells population is $N\left(a, t\right) = S\left(a, t\right) + I\left(a, t\right) + P\left(a, t\right) + C\left(a, t\right)$ and the total number of susceptible cells, infected cells, precancerous cells and cancer cells, respectively, from age $a_{1}$ to $a_{2}$ is denoted by

$\begin{equation*} \begin{matrix} {{N}_{s}}\left( t \right) = \int\limits_{{{a}_{1}}}^{{{a}_{2}}}{S\left( a,t \right)}da, & {{N}_{i}}\left( t \right) = \int\limits_{{{a}_{1}}}^{{{a}_{2}}}{I\left( a,t \right)}da, \\ {{N}_{p}}\left( t \right) = \int\limits_{{{a}_{1}}}^{{{a}_{2}}}{P\left( a,t \right)}da, & {{N}_{c}}\left( t \right) = \int\limits_{{{a}_{1}}}^{{{a}_{2}}}{C\left( a,t \right)}da. \end{matrix} \end{equation*}$

It is well-known that most cases of cervical cancer occurs due to HPV infection, and the virus that have entered the tissue will continue to be reproduced by the infected and cancer cells, and then spread to all cells in the cervical tissue. Furthermore, the virus has ability to infect the susceptible cells and cause the abnormality. The value that shows the rate of infection of a susceptible cells is known as force of infection denoted by $\beta\left(a, t\right)$ . If the infection do not handled properly, the cancer cells can grow faster, uncontrollably, and damaging the cells.

The infection rate of susceptible cells of age $a$ which interacts with the infected and cancer cells of age $b$ is denoted by $\lambda\left(a\right)h\left(b\right)$ . Therefore, the force of infection is defined as

$\begin{equation} \beta \left( a,t \right) = \lambda \left( a \right)\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\frac{\left[ I\left( b,t \right)+C\left( b,t \right) \right]}{N\left( b,t \right)}}\,db, \end{equation}$

(2.1)

where $b \in \left[a_{1}, a_{\sigma}\right]$ . The author in ^[7] showed that precancerous cells are controllable when the lesions develop into precancerous. This condition occurs since precancerous cells can regress ^[1,26] with appropriate treatment, and turning back into infected cells with dormant HPV. By that facts, we exclude the precancerous cells from the force of infection.

Based on the fact that the duration of the cell cycle is more or less the same as the duration of cells infection, we add the the natural growth of the susceptible cells population in our model. In this case, we generalize the model in ^[10] where the natural growth of the susceptible cells population has not been studied.

The age-structured model of cervical cancer cells is formulated by system of first order partial differential equation as follows,

$\begin{eqnarray} \frac{\partial S\left( a,t \right)}{\partial a}+\frac{\partial S\left( a,t \right)}{\partial t} & = & \Lambda N\left(a,t\right)-\mu \left( a \right)S\left( a,t \right)-\beta \left( a,t \right)S\left( a,t \right) \\ \frac{\partial I\left( a,t \right)}{\partial a}+\frac{\partial I\left( a,t \right)}{\partial t} & = & \beta \left( a,t \right)S\left( a,t \right)-\left( \delta \left( a \right)+\mu \left( a \right) \right)I\left( a,t \right) \\ \frac{\partial P\left( a,t \right)}{\partial a}+\frac{\partial P\left( a,t \right)}{\partial t} & = & \delta \left( a \right)I\left( a,t \right)-\left( \gamma \left( a \right)+\mu \left( a \right) \right)P\left( a,t \right) \\ \frac{\partial C\left( a,t \right)}{\partial a}+\frac{\partial C\left( a,t \right)}{\partial t} & = & \gamma \left( a \right)P\left( a,t \right)-\mu \left( a \right)C\left( a,t \right) \end{eqnarray}$

(2.2)

where the boundary and the initial conditions are

$\begin{eqnarray} S\left( {{a}_{1}},t \right) = {{\mu }^{*}}\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{N\left( a,t \right)}da, I\left(a_{1},t \right) = P\left(a_{1},t \right) = C\left(a_{1},t \right) = 0, \\ S\left(a,0 \right) = S_{0}(a), I\left(a,0 \right) = I_{0}(a), P\left(a,0 \right) = P_{0}(a), C\left(a,0 \right) = C_{0}(a). \end{eqnarray}$

(2.3)

The term $\Lambda N\left(a, t\right)$ denotes the density of new susceptible cells' because of the cell division, $\mu^{*}$ is the growth rate of susceptible cells that enter age $a_{1}$ , $\delta \left(a \right)$ is progression rate of the infected to precancerous cells, and $\gamma\left(a\right)$ is progression rate to cancer cell. All variables and parameters are shown in Table 1.

Table 1. Variables and parameters of the model.

Variable/Parameter	Meaning
$S(a, t)$	The density of susceptible cells which depend on age $a$ and time $t$
$I(a, t)$	The density of infected cells which depend on age $a$ and time $t$
$P(a, t)$	The density of precancerous cells which depend on age $a$ and time $t$
$C(a, t)$	The density of cancer cells which depend on age $a$ and time $t$
$\Lambda N(a, t)$	The density of new susceptible cells' because of the cell division
$\lambda(a)$	The age-dependent susceptibility of susceptible cells
$h(b)$	The age-dependent infectiousness caused by infected and cancer cells
$\mu\left(a\right)$	The age-dependent death rate of cells
$\mu^{*}$	The growth rate of susceptible cells that enter age $a_{1}$
$\delta \left(a \right)$	The age-dependent progression rate to precancer
$\gamma\left(a\right)$	The age-dependent progression rate to cancer cell

| Show Table

DownLoad: CSV

Moreover, based on Eqs (2.2) and (2.3), then we obtain

$\begin{align} \frac{\partial N\left( a,t \right)}{\partial a}+\frac{\partial N\left( a,t \right)}{\partial t}& = \frac{\partial S\left( a,t \right)}{\partial a}+\frac{\partial S\left( a,t \right)}{\partial t}+\frac{\partial I\left( a,t \right)}{\partial a}+\frac{\partial I\left( a,t \right)}{\partial t} \\ &\quad +\frac{\partial P\left( a,t \right)}{\partial a}+\frac{\partial P\left( a,t \right)}{\partial t}+\frac{\partial C\left( a,t \right)}{\partial a}+\frac{\partial C\left( a,t \right)}{\partial t} \\ & = \Lambda N\left( a,t \right)-\mu \left( a \right)\left[ S\left( a,t \right)+I\left( a,t \right)+P\left( a,t \right)+C\left( a,t \right) \right] \\ & = \Lambda N\left( a,t \right)-\mu \left( a \right)N\left( a,t \right) \end{align}$

(2.4)

with the boundary conditions

$\begin{align} N\left( {{a}_{1}},t \right)& = S\left( {{a}_{1}},t \right)+I\left( {{a}_{1}},t \right)+P\left( {{a}_{1}},t \right)+C\left( {{a}_{1}},t \right) \\ & = \mu^*\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{N\left( a,t \right)da} \end{align}$

(2.5)

and $N\left(a, 0 \right) = {{N}_{0}}\left(a \right)$ .

3. Steady state conditions

Suppose that $\hat{N}^*\left(a\right)$ is the steady state solution of Eq (2.4) with boundary condition (2.5), then $\frac{dN}{dt} = 0$ , and

$\begin{equation} \frac{d\hat{N}^*\left(a\right)}{da} = \Lambda \hat{N}^*\left(a\right)-\mu \left( a \right)\hat{N}^*\left(a\right) \end{equation}$

(3.1)

where the initial condition is

$\begin{equation} \hat{N}^*\left( {{a}_{1}} \right) = \mu^*\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{\hat{N}^*\left(a\right)da}. \end{equation}$

(3.2)

The solution of (3.1) with initial condition (3.2) is

$\begin{equation} \hat{N}^*\left(a\right) = \hat{N}^*\left( {{a}_{1}} \right){{e}^{\int\limits_{{{a}_{1}}}^{a}{\left( \Lambda -\mu \left( \tau \right) \right)d\tau }}}. \end{equation}$

(3.3)

Since the total population (2.4) reaches to its steady state solution (3.3), see ^[27,28], then we have

$\begin{equation} N\left( a,t \right) = \hat{N}^*\left(a\right) = \hat{N}^*\left( {{a}_{1}} \right){{e}^{\int\limits_{{{a}_{1}}}^{a}{\left( \Lambda -\mu \left( \tau \right) \right)d\tau }}}. \end{equation}$

(3.4)

In this paper, we normalize the System (2.2) to determine the dynamics of the system by using a coordinate transformations as follows,

$\begin{equation} s\left( a,t \right) = \frac{S\left( a,t \right)}{\hat{N}^*\left(a\right)},\,\,i\left( a,t \right) = \frac{I\left( a,t \right)}{\hat{N}^*\left(a\right)},\,\,p\left( a,t \right) = \frac{P\left( a,t \right)}{\hat{N}^*\left(a\right)},\,\,c\left( a,t \right) = \frac{C\left( a,t \right)}{\hat{N}^*\left(a\right)}, \end{equation}$

(3.5)

where $s\left(a, t\right)+i\left(a, t\right)+p\left(a, t\right)+c\left(a, t\right) = 1$ . If the transformation (3.5) is substituted to (2.2) and (2.3), then we obtain

$\begin{align} \frac{\partial s\left( a,t \right)}{\partial a}+\frac{\partial s\left( a,t \right)}{\partial t}& = \Lambda -\Lambda s\left( a,t \right)-\beta^* \left( a,t \right)s\left( a,t \right) \end{align}$

(3.6)

$\begin{align} \frac{\partial i\left( a,t \right)}{\partial a}+\frac{\partial i\left( a,t \right)}{\partial t}& = \beta^* \left( a,t \right)s\left( a,t \right)-\left( \delta \left( a \right)+\Lambda \right)i\left( a,t \right) \end{align}$

(3.7)

$\begin{align} \frac{\partial p\left( a,t \right)}{\partial a}+\frac{\partial p\left( a,t \right)}{\partial t}& = \delta \left( a \right)i\left( a,t \right)-\left( \gamma \left( a \right)+\Lambda \right)p\left( a,t \right) \end{align}$

(3.8)

$\begin{align} \frac{\partial c\left( a,t \right)}{\partial a}+\frac{\partial c\left( a,t \right)}{\partial t}& = \gamma \left( a \right)p\left( a,t \right)-\Lambda c\left( a,t \right) \end{align}$

(3.9)

where

$\begin{equation} \beta^* \left( a,t \right) = \lambda\left(a\right)\int\limits_{a_{1}}^{{a_{\sigma}}}{h\left(b \right)\left[i\left( b,t \right)+c\left(b,t\right)\right]}db. \end{equation}$

(3.10)

The initial and boundary conditions of System (3.6)–(3.9) are

$\begin{align} & s\left( {{a}_{1}},t \right) = \frac{S\left( {{a}_{1}},t \right)}{{{N}}\left( {{a}_{1}},t \right)} = \frac{{{\mu }^{*}}\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{N\left( a,t \right)}da}{{{\mu }^{*}}\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{N\left( a,t \right)da}} = 1, \\ & i\left( {{a}_{1}},t \right) = 0,p\left( {{a}_{1}},t \right) = 0,c\left( {{a}_{1}},t \right) = 0, \end{align}$

and $s\left(a, 0 \right) = {{s}_{0}}\left(a \right), i\left(a, 0 \right) = {{i}_{0}}\left(a \right), p\left(a, 0 \right) = {{p}_{0}}\left(a \right), c\left(a, 0 \right) = {{c}_{0}}\left(a \right)$ .

Suppose that $\hat{s}^*\left(a\right)$ , $\hat{i}^*\left(a\right)$ , $\hat{p}^*\left(a\right)$ and $\hat{c}^*\left(a\right)$ are steady state solution of System (3.6)–(3.9), then the steady state system is

$\begin{eqnarray} \frac{d\hat{s}^*\left(a\right)}{da}& = &\Lambda-\lambda \left( a \right)J\hat{s}^*\left(a\right)-\Lambda \hat{s}^*\left(a\right) \end{eqnarray}$

(3.11)

$\begin{eqnarray} \frac{d\hat{i}^*\left(a\right)}{da}& = &\lambda \left( a \right)J\hat{s}^*\left(a\right)- \delta \left( a \right)\hat{i}^*\left(a\right)-\Lambda \hat{i}^*\left(a\right) \end{eqnarray}$

(3.12)

$\begin{eqnarray} \frac{d\hat{p}^*\left(a\right)}{da}& = &\delta \left( a \right)\hat{i}^*\left(a\right)- \gamma \left( a \right)\hat{p}^*\left(a\right)-\Lambda \hat{p}^*\left(a\right) \end{eqnarray}$

(3.13)

$\begin{eqnarray} \frac{d\hat{c}^*\left(a\right)}{da}& = &\gamma \left( a \right)\hat{p}^*\left(a\right)-\Lambda \hat{c}^*\left(a\right) \end{eqnarray}$

(3.14)

where

$\begin{equation} J = \int\limits_{a_{1}}^{{a_{\sigma}}}{h\left( b \right)\left[\hat{i}^*\left(b\right)+\hat{c}^*\left(b\right)\right]}db. \end{equation}$

(3.15)

We describe the solutions of Eqs (3.11)–(3.14) in the following Lemma.

Lemma 3.1. Let $\left({{\hat{s}}^{*}}\left(a \right), {{\hat{i}}^{*}}\left(a \right), {{\hat{p}}^{*}}\left(a \right), {{\hat{c}}^{*}}\left(a \right) \right)$ is the steady state solution of System (3.6)–(3.9). Then $\left({{\hat{s}}^{*}}\left(a \right), {{\hat{i}}^{*}}\left(a \right), {{\hat{p}}^{*}}\left(a \right), {{\hat{c}}^{*}}\left(a \right) \right)$ satisfies (3.11)–(3.14), where

$\begin{align*} {{{\hat{s}}}^{*}}\left( a \right)& = {{e}^{-\int\limits_{{{a}_{1}}}^{a}{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}+\int\limits_{{{a}_{1}}}^{a}{{{e}^{-\int\limits_{\tau }^{a}{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}\Lambda d\tau }, \\ {{\hat{i}}^{*}}\left( a \right)& = \int\limits_{{{a}_{1}}}^{a}{{{e}^{-\int\limits_{\nu }^{a}{\left( \delta \left( j \right)+\Lambda \right)dj}}}}\lambda \left( \nu \right)J\left[ {{e}^{-\int\limits_{{{a}_{1}}}^{\nu }{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}+\int\limits_{{{a}_{1}}}^{\nu }{{{e}^{-\int\limits_{\tau }^{\nu }{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}\Lambda d\tau } \right]d\nu , \\ {{\hat{p}}^{*}}\left( a \right)& = \int\limits_{{{a}_{1}}}^{a}{{{e}^{-\Lambda \left( a-\eta \right)}}\lambda \left( \eta \right)J\left[ {{e}^{-\int\limits_{{{a}_{1}}}^{\eta }{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}+\int\limits_{{{a}_{1}}}^{\eta }{{{e}^{-\int\limits_{\tau }^{\eta }{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}\Lambda d\tau } \right]}\mathit{\rm{}}\int\limits_{\eta }^{a}{{{e}^{-\int\limits_{\nu }^{\eta }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{v}^{a}{\gamma \left( j \right)dj}}}\delta \left( v \right)d\nu }d\eta , \\ {{{\hat{c}}}^{*}}\left( a \right)& = \int\limits_{{{a}_{1}}}^{a}{{{e}^{-\Lambda \left( a-w \right)}}}\lambda \left( w \right)J\left[ {{e}^{-\int\limits_{{{a}_{1}}}^{w}{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}+\int\limits_{{{a}_{1}}}^{w}{{{e}^{-\int\limits_{\tau }^{w}{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}\Lambda d\tau } \right] \\ & \int\limits_{w}^{a}{\gamma \left( \eta \right)}\int\limits_{\eta }^{w}{{{e}^{-\int\limits_{\nu }^{\eta }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{v}^{w}{\gamma \left( j \right)dj}}}\delta \left( v \right)d\nu }d\eta dw. \end{align*}$

Proof. The steady state condition of System (3.6)–(3.9) is determined by taking $s\left(a, t\right)$ , $i\left(a, t\right)$ , $p\left(a, t\right), c\left(a, t\right)$ each be a constant with respect to time $t$ . Thus, the derivatives with respect to $t$ is zero. If $\frac{ds}{dt} = \frac{di}{dt} = \frac{dp}{dt} = \frac{dc}{dt} = 0$ , then we have the steady state solution by solving each equation in (3.11)–(3.14).

If the value of $J$ in Eq (3.15) is equal to zero, then there is no virus transmission in the cell population. By using Lemma 3.1, for $J = 0$ , we have a disease-free steady state solution, i.e., $E_{0} = \left(\hat{s}\left(a \right), \hat{i}\left(a \right), \hat{p}\left(a \right), \hat{c}\left(a \right) \right) = \left(1, 0, 0, 0 \right)$ .

The local stability of the disease-free steady state solution $E_{0}$ is studied by linerizing System (3.6)–(3.9) near the steady state solution $E_{0}$ by using perturbation technique. Suppose, we consider exponential solutions $\left(\tilde{s}(a)e^{\xi t}, \tilde{i}(a)e^{\xi t}, \tilde{p}(a)e^{\xi t}, \tilde{c}(a)e^{\xi t}\right)$ near the steady state solution ^[29], i.e.,

$\begin{eqnarray} s\left( a,t \right)& = &\hat{s}\left( a \right)+\tilde{s}\left( a \right){{e}^{\xi t}} = 1+\tilde{s}\left( a \right){{e}^{\xi t}} \\ i\left( a,t \right)& = &\hat{i}\left( a \right)+\tilde{i}\left( a \right){{e}^{\xi t}} = 0+\tilde{i}\left( a \right){{e}^{\xi t}} = \tilde{i}\left( a \right){{e}^{\xi t}}\\ p\left( a,t \right)& = &\hat{p}\left( a \right)+\tilde{p}\left( a \right){{e}^{\xi t}} = 0+\tilde{p}\left( a \right){{e}^{\xi t}} = \tilde{p}\left( a \right){{e}^{\xi t}}\\ c\left( a,t \right)& = &\hat{c}\left( a \right)+\tilde{c}\left( a \right){{e}^{\xi t}} = 0+\tilde{c}\left( a \right){{e}^{\xi t}} = \tilde{c}\left( a \right){{e}^{\xi t}} \end{eqnarray}$

(3.16)

where $\xi$ can be a real or complex number. If we substitute Eq (3.16) to (3.10), then we obtain

$\begin{equation} \beta^{*} \left( a,t \right) = \lambda \left( a \right)U{{e}^{\xi t}} \end{equation}$

(3.17)

where

$\begin{equation} U = \int\limits_{a_{1}}^{a_{\sigma} }{h\left( b \right)\left[ \tilde{i}\left( b \right)+\tilde{c}\left( b \right) \right]}db. \end{equation}$

(3.18)

Furthermore, if Eqs (3.16) and (3.17) are substituted to System (3.6)–(3.9), then we have the linear part of the system, i.e.,

$\begin{eqnarray} \frac{d\tilde{s}\left( a \right)}{da}& = & -\tilde{s}\left( a \right)\left( \xi +\Lambda \right)-\lambda \left( a \right)U \end{eqnarray}$

(3.19)

$\begin{eqnarray} \frac{d \tilde{i}\left( a \right)}{d a}& = &\lambda \left( a \right)U-\left( \delta \left( a \right)+\Lambda +\xi \right)\tilde{i}\left( a \right) \end{eqnarray}$

(3.20)

$\begin{eqnarray} \frac{d\tilde{p}\left( a \right)}{da}& = &\delta \left( a \right)\tilde{i}\left( a \right)-\left( \xi +\gamma \left( a \right)+\Lambda \right)\tilde{p}\left( a \right) \end{eqnarray}$

(3.21)

$\begin{eqnarray} \frac{d\tilde{c}\left( a \right)}{da}& = &\gamma \left( a \right)\tilde{p}\left( a \right)-\left( \xi +\Lambda \right)\tilde{c}\left( a \right). \end{eqnarray}$

(3.22)

Furthermore, we should solve Eqs (3.20) and (3.22) to analyze the role of $U$ in Eq (3.18), and we obtain

$\begin{eqnarray} \tilde{i}\left( a \right)& = &\int\limits_{{{a}_{1}}}^{a}{{{e}^{-\int\limits_{\tau }^{a}{\delta \left( j \right)dj}}}{{e}^{-\left( \xi +\Lambda \right)\left( a-\tau \right)}}\lambda \left( \tau \right)Ud\tau } \end{eqnarray}$

(3.23)

$\begin{eqnarray} \tilde{c}\left( a \right)& = &\int\limits_{{{a}_{1}}}^{a}{U}{{e}^{-\left( \xi +\Lambda \right)\left( a-\nu \right)}}\lambda \left( \nu \right)\int\limits_{\nu }^{a}{\gamma \left( \tau \right)} \int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu , \end{eqnarray}$

(3.24)

where $\tau < a$ and $\nu < a$ . The characteristic equation is obtained by substituting (3.23) and (3.24) into (3.18) and gives the following results

$\begin{align} U& = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}{{e}^{-\left( \xi +\Lambda \right)\left( b-\nu \right)}}\lambda \left( \nu \right)Ud\nu } \right.} \\ & \left. +\int\limits_{{{a}_{1}}}^{b}{U}{{e}^{-\left( \xi +\Lambda \right)\left( b-\nu \right)}}\lambda \left( \nu \right)\int\limits_{\nu }^{b}{\gamma \left( \tau \right)}\int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu \right]db. \end{align}$

(3.25)

The solution of Eq (3.25) is $U = 0$ or $U > 0$ that satisfies

$\begin{equation} \hat{Q}\left( \xi \right) = 1, \end{equation}$

(3.26)

where

$\begin{align*} \hat{Q}\left( \xi \right)& = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}{{e}^{-\left( \xi +\Lambda \right)\left( b-\nu \right)}}\lambda \left( \nu \right)d\nu } \right.} \\ & \left. +\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\left( \xi +\Lambda \right)\left( b-\nu \right)}}}\lambda \left( \nu \right)\int\limits_{\nu }^{b}{\gamma \left( \tau \right)}\int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu \right]db. \end{align*}$

Since the threshold value is equal to one, Eq (3.26) has solution $\xi > 0$ if it satisfies $\hat{Q}\left(0 \right) > 1$ , where

$\begin{align} \hat{Q}\left( 0 \right)& = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}{{e}^{-\Lambda \left( b-\nu \right)}}\lambda \left( \nu \right)d\nu } \right.} \\ & \left. +\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\Lambda \left( b-\nu \right)}}}\lambda \left( \nu \right)\int\limits_{\nu }^{b}{\gamma \left( \tau \right)}\int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu \right]db. \end{align}$

(3.27)

Otherwise, Eq (3.26) has solution $\xi \le 0$ . It is explained in Lemma 3.2.

Lemma 3.2. Let $\hat{Q}\left(\xi \right) = 1$ .

1) If $\hat{Q}\left(0 \right) > 1$ , then $\xi > 0$ .

2) If $\hat{Q}\left(0 \right)\le 1$ , then $\xi \le 0$ .

Proof. We will prove the first statement. Suppose $\xi \le 0$ , then for $\hat{Q}\left(\xi \right) = 1$ , we have

$\begin{align*} 1& = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}{{e}^{-\left( \xi +\Lambda \right)\left( b-\nu \right)}}\lambda \left( \nu \right)d\nu } \right.} \\ & \left. +\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\left( \xi +\Lambda \right)\left( b-\nu \right)}}}\lambda \left( \nu \right)\int\limits_{\nu }^{b}{\gamma \left( \tau \right)}\int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu \right]db \\ & \ge \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}{{e}^{-\Lambda \left( b-\nu \right)}}\lambda \left( \nu \right)d\nu } \right.} \\ & \left. +\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\Lambda \left( b-\nu \right)}}}\lambda \left( \nu \right)\int\limits_{\nu }^{b}{\gamma \left( \tau \right)}\int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu \right]db. \end{align*}$

The expression can also be written as

$\begin{align*} &\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}{{e}^{-\Lambda \left( b-\nu \right)}}\lambda \left( \nu \right)d\nu } \right.} \\ & \left. +\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\Lambda \left( b-\nu \right)}}}\lambda \left( \nu \right)\int\limits_{\nu }^{b}{\gamma \left( \tau \right)}\int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu \right]db \le 1. \end{align*}$

It means that $\hat{Q}\left(0 \right) \le 1$ , which is contradiction to our sufficient condition that $\hat{Q}\left(0 \right) > 1$ . Hence, our supposition should be $\xi > 0$ .

Furthermore, we will prove the second statement. Suppose that $\xi > 0$ , then for $\hat{Q}\left(\xi \right) = 1$ , we have

$\begin{align*} 1& = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}{{e}^{-\left( \xi +\Lambda \right)\left( b-\nu \right)}}\lambda \left( \nu \right)d\nu } \right.} \\ & \left. +\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\left( \xi +\Lambda \right)\left( b-\nu \right)}}}\lambda \left( \nu \right)\int\limits_{\nu }^{b}{\gamma \left( \tau \right)}\int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu \right]db \\ & < \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}{{e}^{-\Lambda \left( b-\nu \right)}}\lambda \left( \nu \right)d\nu } \right.} \\ & \left. +\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\Lambda \left( b-\nu \right)}}}\lambda \left( \nu \right)\int\limits_{\nu }^{b}{\gamma \left( \tau \right)}\int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu \right]db. \end{align*}$

It implies that

The inequality can also be written as $\hat{Q}\left(0\right) > 1$ , which is contradiction to our sufficient condition that $\hat{Q}\left(0 \right) \le 1$ . Hence our supposition should be $\xi \le 0$ .

Based on Lemma 3.2, then we define ${{R}_{0}} = \hat{Q}\left(0 \right)$ , where the expression of $\hat{Q}\left(0\right)$ has been shown in (3.27). The value of ${{R}_{0}}$ is called basic reproduction number. It estimates the average number of new infections produced by a particular infected cell or cancer cell in a population. In the next section, we will discuss about the stability of the model near the disease-free and cancerous steady state solution.

4. Disease-free and cancerous steady state solutions

4.1. Local and global stability analysis of disease-free steady state solutions

Before carrying out a stability analysis near the disease-free steady state solution, it is necessary to show that the solution of Eq (3.26) is real and unique.

Proposition 4.1. The equation $\hat{Q}(\xi) = 1$ in (3.26) has a unique real solution, if $\hat{Q}'\left(\xi \right) < 0$ , $\underset{\xi \to \infty }{\mathop{\lim }}\, \hat{Q}\left(\xi \right) = 0$ , and $\underset{\xi \to -\infty } {\mathop{\lim }}\, \hat{Q}\left(\xi \right) = \infty$ .

Proof. Let $\xi$ is a real number. By using Leibniz's differentiation rule, we obtain that $\hat{Q}'\left(\xi \right) < 0$ ,

$\begin{align*} \frac{d\hat{Q}}{d\xi }& = \frac{d}{d\xi }\left( \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}{{e}^{-\left( \xi +\Lambda \right)\left( b-\nu \right)}}\lambda \left( \nu \right)d\nu } \right.} \right. \\ & +\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\left( \xi +\Lambda \right)\left( b-\nu \right)}}}\lambda \left( \nu \right)\int\limits_{\nu }^{b}{\gamma \left( \tau \right)}\left. \left. \int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu \right]db \right) \\ & = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}\left( -\left( b-\nu \right) \right)\left( {{e}^{-\left( \xi +\Lambda \right)\left( b-\nu \right)}} \right)\lambda \left( \nu \right)d\nu }db} \\ & +\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)}\int\limits_{{{a}_{1}}}^{b}{\left( -\left( b-\nu \right) \right)\left( {{e}^{-\left( \xi +\Lambda \right)\left( b-\nu \right)}} \right)}\lambda \left( \nu \right)\int\limits_{\nu }^{b}{\gamma \left( \tau \right)} \\ & \int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu db \\ & = -\left[ \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}\left( b-\nu \right)\left( {{e}^{-\left( \xi +\Lambda \right)\left( b-\nu \right)}} \right)\lambda \left( \nu \right)d\nu }db} \right. \\ & +\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)}\int\limits_{{{a}_{1}}}^{b}{\left( b-\nu \right)\left( {{e}^{-\left( \xi +\Lambda \right)\left( b-\nu \right)}} \right)}\lambda \left( \nu \right)\int\limits_{\nu }^{b}{\gamma \left( \tau \right)} \\ & \left. \int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu db \right] \\ & < 0 , \end{align*}$

where $\nu < b$ . Furthermore, $\underset{\xi \to \infty }{\mathop{\lim\limits}}\, \hat{Q}\left(\xi \right) = 0$ and $\underset{\xi \to -\infty }{\mathop{\lim\limits}}\, \hat{Q}\left(\xi \right) = \infty$ . It means that $\hat{Q}\left(\xi\right)$ is monotonic decreasing function. It implies that the characteristic equation (3.26) has a unique real solution.

Next, we will show that Eq (3.26) has the dominant root, namely $\xi^*$ , see Lemma 4.2. We prove the lemma by following the condition that if there is complex roots obtained from Eq (3.26), then the value of the real part is less than $\xi^*$ .

Lemma 4.2. If $\xi^*$ is the real root of Eq (3.26), then $\xi^*$ is the dominant root.

Proof. Suppose that ${{\xi }^{**}} = x+iy$ is a complex root of Eq (3.26). Then

$\begin{align*} &\hat{Q}\left( {{\xi }^{**}} \right) = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}{{e}^{-\left( x+iy \right)\left( b-\nu \right)}}{{e}^{-\Lambda \left( b-\nu \right)}}\lambda \left( \nu \right)d\nu } \right.} \\ & \left. +\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\left( x+iy \right)\left( b-\nu \right)}}{{e}^{-\Lambda \left( b-\nu \right)}}}\lambda \left( \nu \right)\int\limits_{\nu }^{b}{\gamma \left( \tau \right)}\int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu \right]db \nonumber\\ & = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}{{e}^{-x\left( b-\nu \right)}}\left( \cos \left( \left( b-\nu \right)y \right)-i\sin \left( \left( b-\nu \right)y \right) \right){{e}^{-\Lambda \left( b-\nu \right)}}\lambda \left( \nu \right)d\nu } \right.} \\ & +\int\limits_{{{a}_{1}}}^{b}{{{e}^{-x\left( b-\nu \right)}}\left( \cos \left( \left( b-\nu \right)y \right)-i\sin \left( \left( b-\nu \right)y \right) \right){{e}^{-\Lambda \left( b-\nu \right)}}}\lambda \left( \nu \right) \\ & \left. \int\limits_{\nu }^{b}{\gamma \left( \tau \right)}\int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu \right]db. \end{align*}$

If the real part of $\hat{Q}\left({{\xi }^{**}} \right)$ is equal to $\hat{Q}\left({\xi}^*\right)$ , then

$\begin{align*} \hat{Q}\left( {{\xi }^{*}} \right)& = \operatorname{Re}\hat{Q}\left( {{\xi }^{**}} \right) \\ & = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}{{e}^{-x\left( b-\nu \right)}}\cos \left( \left( b-\nu \right)y \right){{e}^{-\Lambda \left( b-\nu \right)}}\lambda \left( \nu \right)d\nu } \right.} \\ & +\int\limits_{{{a}_{1}}}^{b}{{{e}^{-x\left( b-\nu \right)}}\cos \left( \left( b-\nu \right)y \right){{e}^{-\Lambda \left( b-\nu \right)}}}\lambda \left( \nu \right) \\ & \left. \int\limits_{\nu }^{b}{\gamma \left( \tau \right)}\int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu \right]db \\ & \le \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{\nu }^{b}{\delta \left( j \right)dj}}}{{e}^{-\left( x+\Lambda \right)\left( b-\nu \right)}}\lambda \left( \nu \right)d\nu } \right.} +\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\left( x+\Lambda \right)\left( b-\nu \right)}}}\lambda \left( \nu \right) \\ & \left. \int\limits_{\nu }^{b}{\gamma \left( \tau \right)}\int\limits_{\tau }^{\nu }{\delta \left( w \right){{e}^{-\int\limits_{w}^{\tau }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{w}^{\nu }{\gamma \left( j \right)dj}}}}dwd\tau d\nu \right]db \\ & = \hat{Q}\left( x \right) = \hat{Q}\left( \operatorname{Re}{{\xi }^{**}} \right). \end{align*}$

By using the result of Proposition 4.1, for ${\hat{Q}}'\left(\xi \right) < 0$ , we obtain that $\hat{Q}\left({{\xi }^{*}} \right) \le \hat{Q }\left(\operatorname{Re}{{\xi }^{**}} \right)$ . Thus, we prove that $\operatorname{Re}{{\xi }^{**}} \le {{\xi }^{*}}$ .

Corollary 4.3. The disease-free steady-state solution $E_{0} = \left(1, 0, 0, 0 \right)$ is locally asymptotically stable if $R_{0} < 1$ and it is unstable if $R_{0} > 1$ .

Proof. Based on the Proposition 4.1, we suppose that the unique real value root of Eq (3.26) is $\xi^*$ . If $R_{0} > 1$ , based on Lemma 3.2, for $\hat{Q}\left(0\right) > 1$ , then we have ${\xi}^* > 0$ . Since the characteristic equation has a positive root, then we conclude that the disease-free steady state condition $E_{0}$ is unstable.

Meanwhile, if $R_{0} < 1$ , it means $\hat{Q}\left(0\right) < 1$ , then based on Lemma 3.2, the root of (3.26) is real and negative. Suppose that the root of the characteristic equation is $\xi^*$ , then we have ${\xi}^* < 0$ . Furthermore, based on Lemma 4.2, ${\xi}^*$ is the dominant root of Eq (3.26) and the steady state condition of disease-free $E_{0}$ is asymptotically stable.

The stability of the disease-free steady state solution can be interpreted as follows. If the disease-free steady state solution is stable, i.e., for $R_{0} < 1$ , HPV can be removed from the body, and all cells become normal. This condition can be reached if the initial size of each sub-population is in the basin of attraction of $E_{0}$ . It is in line with the work by ^[14,22,30]. Next, we will analyze global stability around $E_0$ using Fatou's Lemma.

The global stability of the disease-free steady state solution $E_{0}$ will be provided by showing that the value of $s\left(a, t\right)\to 1$ , and the values of $i\left(a, t\right)$ , $p\left(a, t\right)$ , and $c\left(a, t\right)$ tend to zero, as $t \to \infty$ . Equation (3.10) can be rewritten by ${{\beta }^{*}}\left(a, t \right) = \lambda \left(a \right)B\left(t \right)$ , where

$\begin{equation} B\left( t \right) = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ i\left( b,t \right)+c\left( b,t \right) \right]}db. \end{equation}$

(4.1)

If Eq (4.1) is substituted into Eqs (3.6)–(3.9), then we obtain the solution along the characteristic curve ${{c}_{1}} = t-a$ as follows,

$\begin{align} s\left( a,t \right)& = {{e}^{-\int\limits_{{{a}_{1}}}^{a}{\Lambda +\lambda \left( j \right)B\left( j+t-a \right)}dj}}\int\limits_{{{a}_{1}}}^{a}{\Lambda {{e}^{\int\limits_{{{a}_{1}}}^{\tau }{\Lambda +\lambda \left( j \right)B\left( j+t-a \right)}dj}}d\tau } +{{e}^{-\int\limits_{{{a}_{1}}}^{a}{\Lambda +\lambda \left( j \right)B\left( j+t-a \right)}dj}} \end{align}$

(4.2)

$\begin{align} i\left( a,t \right)& = {{e}^{-\int\limits_{{{a}_{1}}}^{a}{\left( \delta \left( j \right)+\Lambda \right)}dj}}\int\limits_{{{a}_{1}}}^{a}{{{e}^{\int\limits_{{{a}_{1}}}^{\tau }{\left( \delta \left( j \right)+\Lambda \right)}dj}}\lambda \left( \tau \right)B\left( \tau +t-a \right)s\left( \tau ,\tau +t-a \right)}d\tau \end{align}$

(4.3)

$\begin{align} p\left( a,t \right)& = {{e}^{-\int\limits_{{{a}_{1}}}^{a}{\left( \gamma \left( j \right)+\Lambda \right)dj}}}\int\limits_{{{a}_{1}}}^{a}{{{e}^{\int\limits_{{{a}_{1}}}^{\eta }{\left( \gamma \left( j \right)+\Lambda \right)dj}}}}\delta \left( \eta \right){{e}^{-\int\limits_{{{a}_{1}}}^{\eta }{\left( \delta \left( j \right)+\Lambda \right)}dj}} \\ & \int\limits_{{{a}_{1}}}^{\eta }{{{e}^{\int\limits_{{{a}_{1}}}^{\tau }{\left( \delta \left( j \right)+\Lambda \right)}dj}}\lambda \left( \tau \right)B\left( \tau +t-a \right)s\left( \tau ,\tau +t-a \right)d\tau }d\eta \end{align}$

(4.4)

$\begin{align} c\left( a,t \right)& = {{e}^{-\int\limits_{{{a}_{1}}}^{a}{\Lambda dj}}}\int\limits_{{{a}_{1}}}^{a}{{{e}^{\int\limits_{{{a}_{1}}}^{w}{\Lambda dj}}}}\gamma \left( w \right){{e}^{-\int\limits_{{{a}_{1}}}^{w}{\left( \gamma \left( j \right)+\Lambda \right)dj}}} \\ & \int\limits_{{{a}_{1}}}^{w}{{{e}^{\int\limits_{{{a}_{1}}}^{\eta }{\left( \gamma \left( j \right)+\Lambda \right)dj}}}\delta \left( \eta \right){{e}^{-\int\limits_{{{a}_{1}}}^{\eta }{\left( \delta \left( j \right)+\Lambda \right)}dj}}\int\limits_{{{a}_{1}}}^{\eta }{{{e}^{\int\limits_{{{a}_{1}}}^{\tau }{\left( \delta \left( j \right)+\Lambda \right)}dj}}\lambda \left( \tau \right)B\left( \tau +t-a \right)s\left( \tau ,\tau +t-a \right)d\tau }d\eta }dw. \end{align}$

(4.5)

Now, Eqs (4.3) and (4.5) in $E_0$ are substituted into Eq (4.1), then we have

$\begin{align*} B\left( t \right)& = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ {{e}^{-\int\limits_{{{a}_{1}}}^{a}{\left( \delta \left( j \right)+\Lambda \right)}dj}}\int\limits_{{{a}_{1}}}^{a}{{{e}^{\int\limits_{{{a}_{1}}}^{\tau }{\left( \delta \left( j \right)+\Lambda \right)}dj}}\lambda \left( \tau \right)B\left( \tau +t-a \right)}d\tau \right.} \\ & +{{e}^{-\int\limits_{{{a}_{1}}}^{a}{\Lambda dj}}}\int\limits_{{{a}_{1}}}^{a}{{{e}^{\int\limits_{{{a}_{1}}}^{w}{\Lambda dj}}}\gamma \left( w \right){{e}^{-\int\limits_{{{a}_{1}}}^{w}{\left( \gamma \left( j \right)+\Lambda \right)dj}}}\int\limits_{{{a}_{1}}}^{w}{{{e}^{\int\limits_{{{a}_{1}}}^{\eta }{\left( \gamma \left( j \right)+\Lambda \right)dj}}}\delta \left( \eta \right)}} \\ & \left. {{e}^{-\int\limits_{{{a}_{1}}}^{\eta }{\left( \delta \left( j \right)+\Lambda \right)}dj}}\int\limits_{{{a}_{1}}}^{\eta }{{{e}^{\int\limits_{{{a}_{1}}}^{\tau }{\left( \delta \left( j \right)+\Lambda \right)}dj}}\lambda \left( \tau \right)B\left( \tau +t-a \right)d\tau d\eta dw} \right]db \end{align*}$

If $V\left(a \right) = \lambda \left(a \right)\underset{t\to \infty }{\mathop{\lim \sup }}\, B\left(t \right)$ , then by using Fatou's Lemma and by applying the supremum limit of $t \to \infty$ in both sides, we obtain

$\begin{equation} V\left( a \right)\le \lambda \left( a \right){{C}_{1}} \end{equation}$

(4.6)

where

$\begin{align} {{C}_{1}}& = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ {{e}^{-\int\limits_{{{a}_{1}}}^{a}{\left( \delta \left( j \right)+\Lambda \right)}dj}}\int\limits_{{{a}_{1}}}^{a}{{{e}^{\int\limits_{{{a}_{1}}}^{\tau }{\left( \delta \left( j \right)+\Lambda \right)}dj}}V\left( \tau \right)}d\tau \right.} \\ & +{{e}^{-\int\limits_{{{a}_{1}}}^{a}{\Lambda dj}}}\int\limits_{{{a}_{1}}}^{a}{{{e}^{\int\limits_{{{a}_{1}}}^{w}{\Lambda dj}}}\gamma \left( w \right){{e}^{-\int\limits_{{{a}_{1}}}^{w}{\left( \gamma \left( j \right)+\Lambda \right)dj}}}} \\ & \left. \int\limits_{{{a}_{1}}}^{w}{{{e}^{\int\limits_{{{a}_{1}}}^{\eta }{\left( \gamma \left( j \right)+\Lambda \right)dj}}}\delta \left( \eta \right){{e}^{-\int\limits_{{{a}_{1}}}^{\eta }{\left( \delta \left( j \right)+\Lambda \right)}dj}}\int\limits_{{{a}_{1}}}^{\eta }{{{e}^{\int\limits_{{{a}_{1}}}^{\tau }{\left( \delta \left( j \right)+\Lambda \right)}dj}}V\left( \tau \right)d\tau }d\eta }dw \right]. \end{align}$

(4.7)

By substitution of Eq (4.6) to (4.7), then we have

$\begin{align*} {{C}_{1}} &\le \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ {{e}^{-\int\limits_{{{a}_{1}}}^{a}{\left( \delta \left( j \right)+\Lambda \right)}dj}}\int\limits_{{{a}_{1}}}^{a}{{{e}^{\int\limits_{{{a}_{1}}}^{\tau }{\left( \delta \left( j \right)+\Lambda \right)}dj}}\lambda \left( \tau \right){{C}_{1}}}d\tau \right.} \\ & +{{e}^{-\int\limits_{{{a}_{1}}}^{a}{\Lambda dj}}}\int\limits_{{{a}_{1}}}^{a}{{{e}^{\int\limits_{{{a}_{1}}}^{w}{\Lambda dj}}}\gamma \left( w \right)}{{e}^{-\int\limits_{{{a}_{1}}}^{w}{\left( \gamma \left( j \right)+\Lambda \right)dj}}} \\ & \left. \int\limits_{{{a}_{1}}}^{w}{{{e}^{\int\limits_{{{a}_{1}}}^{\eta }{\left( \gamma \left( j \right)+\Lambda \right)dj}}}\delta \left( \eta \right){{e}^{-\int\limits_{{{a}_{1}}}^{\eta }{\left( \delta \left( j \right)+\Lambda \right)}dj}}\int\limits_{{{a}_{1}}}^{\eta }{{{e}^{\int\limits_{{{a}_{1}}}^{\tau }{\left( \delta \left( j \right)+\Lambda \right)}dj}}\lambda \left( \tau \right){{C}_{1}}d\tau }d\eta }dw \right] \end{align*}$

or ${{C}_{1}}\le {{C}_{1}}{{R}_{0}}$ . It means that ${{C}_{1}}\left(1-{{R}_{0}} \right)\le 0$ . Since ${{R}_{0}} < 1$ and ${{C}_{1}}$ is nonnegative, then the only possibility of ${{C}_{1}}$ is zero. By applying this result, we have $V\left(a \right) = 0$ and $\underset{t\to \infty }{\mathop{\lim \sup }}\, B\left(t \right) = 0$ , and then we have $\underset{t\to \infty }{\mathop{\lim\limits }}\, i\left(t \right) = \underset{t\to \infty }{\mathop{\lim\limits }}\, p\left(t \right) = \underset{t\to \infty }{\mathop{\lim\limits }}\, c\left(t \right) = 0$ and $\underset{t\to \infty }{\mathop{\lim\limits }}\, s\left(t \right) = 1$ . Finally, we conclude the global stability conditions of the disease-free steady state solution in Theorem 4.4.

Theorem 4.4. If $R_{0} < 1$ , then the disease-free steady state solution $E_{0} = \left(1, 0, 0, 0 \right)$ is globally asymptotically stable.

4.2. Existence of cancerous steady state solution

In Corollary 4.3, we have shown that for ${{R}_{0}} > 1$ , the disease-free steady state solution is unstable and there exists a cancerous steady state solution. The existence conditions of the cancerous steady state solution is explained in Theorem 4.5.

Theorem 4.5. If ${{R}_{0}} > 1$ , then there exists a cancerous steady state solution.

Proof. Recall the steady state solutions of System (3.6)–(3.9) in Lemma 3.1. Suppose that $J$ in Eq (3.15) is not equal to zero, then we have $E_{1} = \left({{\hat{s}}^{*}}\left(a \right), {{\hat{i}}^{*}}\left(a \right), {{\hat{p}}^{*}}\left(a \right), {{\hat{c}}^{*}}\left(a \right) \right)$ and $J = J\hat{Q}\left(J \right)$ where

$\begin{align} \hat{Q}\left( J \right)& = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)}\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\int\limits_{w}^{b}{\left( \delta \left( j \right)+\Lambda \right)dj}}}\lambda \left( w \right)} \left[ {{e}^{-\int\limits_{{{a}_{1}}}^{w}{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}+\int\limits_{{{a}_{1}}}^{w}{{{e}^{-\int\limits_{\tau }^{w}{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}\Lambda d\tau } \right]dwdb \\ & +\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)}\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\Lambda \left( b-w \right)}}}\lambda \left( w \right)\left[ {{e}^{-\int\limits_{{{a}_{1}}}^{w}{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}+\int\limits_{{{a}_{1}}}^{w}{{{e}^{-\int\limits_{\tau }^{w}{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}\Lambda d\tau } \right] \\ & \int\limits_{w}^{b}{\gamma \left( \eta \right)}\int\limits_{\eta }^{w}{{{e}^{-\int\limits_{\nu }^{\eta }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{v}^{w}{\gamma \left( j \right)dj}}}\delta \left( v \right)d\nu }d\eta dwdb. \end{align}$

(4.8)

The solution of $J = J\hat{Q}\left(J \right)$ is $J = 0$ or $J > 0$ that satisfies

$\begin{equation} \hat{Q}\left( J \right) = 1. \end{equation}$

(4.9)

By changing the order of integration, for $J = 0$ , Eq (4.8) can be written as

We see that System (3.6)–(3.9) has a cancerous steady state solution $E_{1}$ , if Eq (4.9) holds for $J$ . We denote that $\hat{Q}\left(0 \right) = R_{0}$ . Since $i\left(b \right)+c\left(b \right) < 1$ , then we obtain

$\begin{equation*} \hat{Q}\left( J \right) = \frac{1}{J}\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ i\left( b \right)+c\left( b \right) \right]}db < \frac{1}{J}\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)}db \end{equation*}$

where $\hat{Q}\left(J \right)\to 0$ as $J \to \infty$ . Thus, for $R_{0} > 1$ , Eq (4.9) holds for a unique positive value of $J$ .

4.3. Local stability of cancerous steady state solution

We use a perturbation technique to study the local stability of the cancerous steady state solution. Based on the steady state solution of System (3.6)–(3.9) that has been attained in Lemma 3.1, we obtain the perturbation transformations as follows,

$\begin{align} s\left( a,t \right)& = {{{\hat{s}}}^{*}}\left( a \right)+{{{\tilde{s}}}^{*}}\left( a \right){{e}^{\xi t}} \\ & = {{e}^{-\int\limits_{{{a}_{1}}}^{a}{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}+\int\limits_{{{a}_{1}}}^{a}{{{e}^{-\int\limits_{\tau }^{a}{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}\Lambda d\tau }+{{{\tilde{s}}}^{*}}\left( a \right){{e}^{\xi t}} \end{align}$

(4.10)

$\begin{align} i\left( a,t \right)& = {{{\hat{i}}}^{*}}\left( a \right)+{{{\tilde{i}}}^{*}}\left( a \right){{e}^{\xi t}} \\ & = \int\limits_{{{a}_{1}}}^{a}{{{e}^{-\int\limits_{\nu }^{a}{\left( \delta \left( j \right)+\Lambda \right)dj}}}\lambda \left( \nu \right)J\left[ {{e}^{-\int\limits_{{{a}_{1}}}^{\nu }{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}+\int\limits_{{{a}_{1}}}^{\nu }{{{e}^{-\int\limits_{\tau }^{\nu }{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}\Lambda d\tau } \right]d\nu } +{{{\tilde{i}}}^{*}}\left( a \right){{e}^{\xi t}} \end{align}$

(4.11)

$\begin{align} p\left( a,t \right)& = {{{\hat{p}}}^{*}}\left( a \right)+{{{\tilde{p}}}^{*}}\left( a \right){{e}^{\xi t}} \\ & = \int\limits_{{{a}_{1}}}^{a}{{{e}^{-\Lambda \left( a-\eta \right)}}\lambda \left( \eta \right)J\left[ {{e}^{-\int\limits_{{{a}_{1}}}^{\eta }{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}+\int\limits_{{{a}_{1}}}^{\eta }{{{e}^{-\int\limits_{\tau }^{\eta }{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}\Lambda d\tau } \right]} \\ & \int\limits_{\eta }^{a}{{{e}^{-\int\limits_{\nu }^{\eta }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{v}^{a}{\gamma \left( j \right)dj}}}\delta \left( v \right)d\nu }d\eta +{{{\tilde{p}}}^{*}}\left( a \right){{e}^{\xi t}} \end{align}$

(4.12)

$\begin{align} c\left( a,t \right)& = {{{\hat{c}}}^{*}}\left( a \right)+{{{\tilde{c}}}^{*}}\left( a \right){{e}^{\xi t}} \\ & = \int\limits_{{{a}_{1}}}^{a}{{{e}^{-\Lambda \left( a-w \right)}}}\lambda \left( w \right)J\left[ {{e}^{-\int\limits_{{{a}_{1}}}^{w}{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}+\int\limits_{{{a}_{1}}}^{w}{{{e}^{-\int\limits_{\tau }^{w}{\left( \lambda \left( j \right)J+\Lambda \right)dj}}}\Lambda d\tau } \right] \\ & \int\limits_{w}^{a}{\gamma \left( \eta \right)}\int\limits_{\eta }^{w}{{{e}^{-\int\limits_{\nu }^{\eta }{\delta \left( j \right)dj}}}{{e}^{-\int\limits_{v}^{w}{\gamma \left( j \right)dj}}}\delta \left( v \right)d\nu }d\eta dw +{{{\tilde{c}}}^{*}}\left( a \right){{e}^{\xi t}}. \end{align}$

(4.13)

If Eqs (4.11) and (4.13) are substituted into Eq (3.10) and based on Eq (3.15), we have

$\begin{align} {{\beta }^{*}}\left( a,t \right)& = \lambda \left( a \right)\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ {{{\hat{i}}}^{*}}\left( b \right)+{{{\tilde{i}}}^{*}}\left( b \right){{e}^{\xi t}}+{{{\hat{c}}}^{*}}\left( b \right)+{{{\tilde{c}}}^{*}}\left( b \right){{e}^{\xi t}} \right]}db = \lambda \left( a \right)J+\lambda \left( a \right){{U}^{*}}{{e}^{\xi t}}, \end{align}$

(4.14)

where

$\begin{equation} {{U}^{*}}\, = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ {{{\tilde{i}}}^{*}}\left( b \right)+{{{\tilde{c}}}^{*}}\left( b \right) \right]}db. \end{equation}$

(4.15)

Then we have the linearized perturbed system in ${{\tilde{s}}^{*}}\left(a \right), {{\tilde{i}}^{*}}\left(a \right), {{\tilde{p}}^{*}}\left(a \right)$ , and ${{\tilde{c}}^{*}}\left(a \right)$ , as follows,

$\begin{align*} \frac{d{{\tilde{s}}^{*}}\left( a \right)}{da}& = -{{\tilde{s}}^{*}}\left( a \right)\left[ \Lambda +\xi +\lambda \left( a \right)J \right]-\lambda \left( a \right){{U}^{*}}{{\hat{s}}^{*}}\left( a \right) \\ \frac{d{{\tilde{i}}^{*}}\left( a \right)}{da}& = -\left( \xi +\delta \left( a \right)+\Lambda \right){{\tilde{i}}^{*}}\left( a \right)+\lambda \left( a \right)\left( {{U}^{*}}{{{\hat{s}}}^{*}}\left( a \right)+J{{{\tilde{s}}}^{*}}\left( a \right) \right) \\ \frac{d{{\tilde{p}}^{*}}\left( a \right)}{da}& = \delta \left( a \right){{\tilde{i}}^{*}}\left( a \right)-\left( \xi +\gamma \left( a \right)+\Lambda \right){{\tilde{p}}^{*}}\left( a \right) \\ \frac{d {{{\tilde{c}}}^{*}}\left( a \right)}{d a}& = \gamma \left( a \right){{\tilde{p}}^{*}}\left( a \right)-\left( \Lambda +\xi \right){{\tilde{c}}^{*}}\left( a \right). \end{align*}$

Since $U^* \ne 0$ , then we can use the transformation

$\begin{equation*} \bar{s} = \frac{{{{\tilde{s}}}^{*}}}{{{U}^{*}}},\quad\bar{i} = \frac{{{{\tilde{i}}}^{*}}}{{{U}^{*}}},\quad \bar{p} = \frac{{{{\tilde{p}}}^{*}}}{{{U}^{*}}},\quad \rm{and}\quad \bar{c} = \frac{{{{\tilde{c}}}^{*}}}{{{U}^{*}}}, \end{equation*}$

to remove $U^*$ from the system. Thus, we have

$\begin{align} \frac{d\bar{s}\left( a \right)}{da}& = -\bar{s}\left( a \right)\left[ \Lambda +\xi +\lambda \left( a \right)J \right]-\lambda \left( a \right){{{\hat{s}}}^{*}}\left( a \right) \end{align}$

(4.16)

$\begin{align} \frac{d\bar{i}\left( a \right)}{da}& = -\left( \xi +\delta \left( a \right)+\Lambda \right)\bar{i}\left( a \right)+\lambda \left( a \right)\left( {{{\hat{s}}}^{*}}\left( a \right)+J\bar{s}\left( a \right) \right) \end{align}$

(4.17)

$\begin{align} \frac{d\bar{p}\left( a \right)}{da}& = \delta \left( a \right)\bar{i}\left( a \right)-\left( \xi +\gamma \left( a \right)+\Lambda \right)\bar{p}\left( a \right) \end{align}$

(4.18)

$\begin{align} \frac{d\bar{c}\left( a \right)}{da}& = \gamma \left( a \right)\bar{p}\left( a \right)-\left( \Lambda +\xi \right)\bar{c}\left( a \right). \end{align}$

(4.19)

By using the transformation above, Eq (4.15) becomes

$\begin{equation} 1\, = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)\left[ \bar{i}\left( b \right)+\bar{c}\left( b \right) \right]}db, \end{equation}$

(4.20)

and the solution of (4.16)–(4.19) is

$\begin{align} \bar{s}\left( a \right)& = \int\limits_{{{a}_{1}}}^{a}{-\lambda \left( \tau \right){{{\hat{s}}}^{*}}\left( \tau \right)}{{e}^{-\int\limits_{\tau }^{a}{\left( \Lambda +\xi +\lambda \left( j \right)J \right)dj}}}d\tau \end{align}$

(4.21)

$\begin{align} \bar{i}(a)& = \int\limits_{{{a}_{1}}}^{a}{{{e}^{-\left( \xi +\Lambda \right)\left( a-w \right)}}}\lambda \left( w \right){{{\hat{s}}}^{*}}\left( w \right) \left[ {{e}^{-\int\limits_{w}^{a}{\delta \left( j \right)}dj}}-J\int\limits_{w}^{a}{\lambda \left( \tau \right)}{{e}^{-\int\limits_{\tau }^{w}{\lambda \left( j \right)J}dj}}{{e}^{-\int\limits_{\tau }^{a}{\delta \left( j \right)}dj}}d\tau \right]dw \end{align}$

(4.22)

$\begin{align} \bar{p}\left( a \right)& = \int\limits_{{{a}_{1}}}^{a}{{{e}^{-\int\limits_{\tau }^{a}{\left( \xi +\gamma \left( j \right)+\Lambda \right)dj}}}}\delta \left( \tau \right)\bar{i}\left( \tau \right)d\tau \end{align}$

(4.23)

$\begin{align} \bar{c}(a) & = \int\limits_{{{a}_{1}}}^{a}{{{e}^{-\left( \xi +\Lambda \right)\left( a-w \right)}}\lambda \left( w \right){{{\hat{s}}}^{*}}\left( w \right)}\int\limits_{w}^{a}{\gamma \left( \eta \right)}\int\limits_{\eta }^{w}{{{e}^{-\int\limits_{\nu }^{w}{\gamma \left( j \right)dj}}}}\delta \left( \nu \right) \\ & \left( {{e}^{-\int\limits_{v}^{\eta }{\delta \left( j \right)}dj}}-J\int\limits_{\nu }^{\eta }{\lambda \left( \tau \right)}{{e}^{-\int\limits_{\tau }^{v}{\lambda \left( j \right)Jdj}}}{{e}^{-\int\limits_{\tau }^{\eta }{\delta \left( j \right)}dj}}d\tau \right)dvd\eta dw. \end{align}$

(4.24)

If Eqs (4.22) and (4.24) are substituted into (4.20), and the right-hand side of (4.20) is denoted by $\bar{Q}\left(\xi\right)$ , then we get

$\begin{equation} \bar{Q}\left( \xi \right) = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)}\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\left( \xi +\Lambda \right)\left( b-w \right)}}}\lambda \left( w \right){{\hat{s}}^{*}}\left( w \right)\Omega \left( b,w \right)dwdb \end{equation}$

(4.25)

where

$\begin{align} \Omega \left( b,w \right)& = {{e}^{-\int\limits_{w}^{b}{\delta \left( j \right)}dj}}-J\int\limits_{w}^{b}{\lambda \left( \tau \right)}{{e}^{-\int\limits_{\tau }^{w}{\lambda \left( j \right)J}dj}}{{e}^{-\int\limits_{\tau }^{b}{\delta \left( j \right)}dj}}d\tau +\int\limits_{w}^{b}{\gamma \left( \eta \right)}\int\limits_{\eta }^{w}{{{e}^{-\int\limits_{\nu }^{w}{\gamma \left( j \right)dj}}}}\delta \left( \nu \right)\left( {{e}^{-\int\limits_{v}^{\eta }{\delta \left( j \right)}dj}} \right. \\ & \left. -J\int\limits_{\nu }^{\eta }{\lambda \left( \tau \right)}{{e}^{-\int\limits_{\tau }^{v}{\lambda \left( j \right)Jdj}}}{{e}^{-\int\limits_{\tau }^{\eta }{\delta \left( j \right)}dj}}d\tau \right)dvd\eta. \end{align}$

(4.26)

Let

$\begin{equation} \Psi \left( b,w \right) = {{e}^{-\int\limits_{w}^{b}{\delta \left( j \right)}dj}}-J\int\limits_{w}^{b}{\lambda \left( \tau \right)}{{e}^{-\int\limits_{\tau }^{w}{\lambda \left( j \right)J}dj}}{{e}^{-\int\limits_{\tau }^{b}{\delta \left( j \right)}dj}}d\tau \end{equation}$

(4.27)

and $\Phi \left(w, \nu \right) = {{e}^{-\int\limits_{v}^{w}{\gamma \left(j \right)dj}}}\delta \left(\nu \right)$ , then Eq (4.26) becomes

$\begin{equation*} \Omega \left( b,w \right) = \Psi \left( b,w \right)+\int\limits_{w}^{b}{\gamma \left( \eta \right)}\int\limits_{\eta }^{w}{\Phi \left( w,\nu \right)\Psi \left( \eta ,\nu \right)d\nu d\eta }. \end{equation*}$

Next, we will show the conditions that the root of $\bar{Q}\left(\xi \right) = 1$ is unique and real.

Theorem 4.6. If $\Psi \left(b, w \right)\ge 0$ where $a_1 \le w\le b\le {{a}_{\sigma }}$ , then

${\bar{Q}}'\left(\xi \right) < 0, \underset{\xi \to \infty }{\mathop{\lim\limits }}\, \bar{Q}\left(\xi \right) = 0, \underset{\xi \to -\infty }{\mathop{\lim\limits }}\, \bar{Q}\left(\xi \right) = \infty$ , and $\bar{Q}\left(0 \right) < 1$ .

Proof. For the first statement. If $\Psi \left(b, w \right)\ge 0$ , then $\Omega \left(b, w \right) \ge 0$ . Based on (4.25), then $\bar{Q}\left(\xi \right) \ge 0$ . Thus

$\begin{align*} {\bar{Q}}'\left( \xi \right)& = \frac{d}{d\xi }\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)}\int\limits_{{{a}_{1}}}^{b}{{{e}^{-\left( \xi +\Lambda \right)\left( b-w \right)}}}\lambda \left( w \right){{{\hat{s}}}^{*}}\left( w \right)\Omega \left( b,w \right)dwdb \\ & = -\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)}\int\limits_{{{a}_{1}}}^{b}{\left( b-w \right){{e}^{-\left( \xi +\Lambda \right)\left( b-w \right)}}}\lambda \left( w \right){{{\hat{s}}}^{*}}\left( w \right)\Omega \left( b,w \right)dwdb \\ & < 0. \end{align*}$

Furthermore, based on (4.25) it is clear that $\underset{\xi \to \infty }{\mathop{\lim\limits }}\, \bar{Q}\left(\xi \right) = 0, \underset{\xi \to -\infty }{\mathop{\lim\limits }}\, \bar{Q}\left(\xi \right) = \infty$ .

Now we will prove for the second statement. Based on (4.9), then we have

$\begin{align*} \bar{Q}\left( 0 \right)& = 1+J\int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)}\int\limits_{w}^{b}{\bar{s}\left( \eta \right)\lambda \left( \eta \right)} \\ & \left( {{e}^{-\int\limits_{\tau }^{b}{\delta \left( j \right)}dj}}+\int\limits_{\eta }^{w}{\gamma \left( \nu \right)}\int\limits_{\nu }^{\eta }{\delta \left( \tau \right)}{{e}^{-\int\limits_{\tau }^{\eta }{\delta \left( j \right)}dj}}{{e}^{-\int\limits_{\tau }^{w}{\gamma \left( j \right)dj}}}d\tau dv \right)d\eta db. \end{align*}$

Since based on (4.21) that $\bar{s}\left(a\right) < 0$ , then $\bar{Q}\left(0\right) < 1$ . It completes the proof.

Based on Theorem 4.6, the function of $\bar{Q}$ satisfies ${\bar{Q}}'\left(\xi \right) < 0, \underset{\xi \to \infty }{\mathop{\lim\limits }}\, \bar{Q}\left(\xi \right) = 0, \underset{\xi \to -\infty }{\mathop{\lim\limits }}\, \bar{Q}\left(\xi \right) = \infty$ . It follows that the function of $\bar{Q}$ is a monotone decreasing. Furthermore, derived from the properties of $\bar{Q}$ , Eq (4.20) has a unique root. We also have proved in Theorem 4.6 that $\bar{Q}\left(0\right) < 1$ , then the characteristic equation (4.20) has a negative root. Next it is proved that the negative root obtained is the dominant root.

Lemma 4.7. If $\xi^*$ is real root of $\bar{Q}\left(\xi \right) = 1$ , then $\xi^*$ is the dominant root.

Proof. Proving this lemma means if there is a complex root obtained from $\bar{Q}\left(\xi \right) = 1$ , then this root has a real part whose value is less than $\xi^*$ . Suppose $\bar{\xi } = \bar{x}+i\bar{y}$ is the root $\bar{Q}\left(\xi \right) = 1$ in the form of a complex number, then we have

$\begin{align*} \bar{Q}\left( {\bar{\xi }} \right)& = \bar{Q}\left( \bar{x}+i\bar{y} \right) \\ & = \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)}\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\bar{x}\left( b-w \right)}}\left( \cos \left( \left( b-w \right)\bar{y} \right)-i\sin \left( \left( b-w \right)\bar{y} \right) \right)}{{e}^{-\Lambda \left( b-w \right)}} \right. \\ & \lambda \left( w \right){{{\hat{s}}}^{*}}\left( w \right)\left( {{e}^{-\int\limits_{w}^{b}{\delta \left( j \right)}dj}}-J\int\limits_{w}^{b}{\lambda \left( \tau \right)}{{e}^{-\int\limits_{w}^{\tau }{\lambda \left( j \right)J}dj}}{{e}^{-\int\limits_{\tau }^{b}{\delta \left( j \right)}dj}}d\tau \right. \\ & +\int\limits_{w}^{b}{\gamma \left( \eta \right)}\int\limits_{w}^{b}{{{e}^{-\int\limits_{v}^{b}{\gamma \left( j \right)dj}}}\delta \left( v \right)}{{e}^{-\int\limits_{w}^{\eta }{\delta \left( j \right)}dj}}dvd\eta -J\int\limits_{w}^{b}{\gamma \left( \eta \right)}\int\limits_{w}^{b}{{{e}^{-\int\limits_{v}^{b}{\gamma \left( j \right)dj}}}\delta \left( v \right)} \\ & \left. \left. \int\limits_{w}^{v}{\lambda \left( \tau \right)}{{e}^{-\int\limits_{w}^{\tau }{\lambda \left( j \right)Jdj}}}{{e}^{-\int\limits_{\tau }^{v}{\delta \left( j \right)}dj}}d\tau dvd\eta \right)dw \right]db. \end{align*}$

As a result, the real part of $\bar{Q}\left({\bar{\xi }} \right)$ is equal to $\bar{Q}\left({{{\bar{\xi }}}^{ *}} \right)$ , then we acquire

$\begin{align*} \bar{Q}\left( {{{\bar{\xi }}}^{*}} \right)& = \operatorname{Re}\bar{Q}\left( {\bar{\xi }} \right) \\ & \le \int\limits_{{{a}_{1}}}^{{{a}_{\sigma }}}{h\left( b \right)}\left[ \int\limits_{{{a}_{1}}}^{b}{{{e}^{-\bar{x}\left( b-w \right)}}}{{e}^{-\Lambda \left( b-w \right)}}\lambda \left( w \right){{{\hat{s}}}^{*}}\left( w \right) \right. \\ & \left( {{e}^{-\int\limits_{w}^{b}{\delta \left( j \right)}dj}}-J\int\limits_{w}^{b}{\lambda \left( \tau \right)}{{e}^{-\int\limits_{w}^{\tau }{\lambda \left( j \right)J}dj}}{{e}^{-\int\limits_{\tau }^{b}{\delta \left( j \right)}dj}}d\tau \right. \\ & +\int\limits_{w}^{b}{\gamma \left( \eta \right)}\int\limits_{w}^{b}{{{e}^{-\int\limits_{v}^{b}{\gamma \left( j \right)dj}}}\delta \left( v \right)}{{e}^{-\int\limits_{w}^{\eta }{\delta \left( j \right)}dj}}dvd\eta -J\int\limits_{w}^{b}{\gamma \left( \eta \right)}\int\limits_{w}^{b}{{{e}^{-\int\limits_{v}^{b}{\gamma \left( j \right)dj}}}\delta \left( v \right)} \\ & \left. \left. \int\limits_{w}^{v}{\lambda \left( \tau \right)}{{e}^{-\int\limits_{w}^{\tau }{\lambda \left( j \right)Jdj}}}{{e}^{-\int\limits_{\tau }^{v}{\delta \left( j \right)}dj}}d\tau dvd\eta \right)dw \right]db \\ & = \bar{Q}\left( {\bar{x}} \right) = \bar{Q}\left( \operatorname{Re}\bar{\xi } \right). \end{align*}$

Based on the Theorem 4.6, we derive ${\bar{Q}}'\left(\xi \right) < 0$ and $\bar{Q}\left({{{\bar{\xi }}}^{*}} \right)\le \bar{Q}\left(\operatorname{Re}\bar{\xi } \right)$ , so it proves $\operatorname{Re}\bar{\xi }\le {{\bar{\xi }}^{*}}$ .

Furthermore, the stability of the cancerous steady state solution can be seen in Corollary 4.8.

Corollary 4.8. If assumption in Theorems 4.5, 4.6, and Lemma 4.7 hold, then the cancerous steady state solution is locally asymptotically stable.

Proof. Let $R_0 > 1$ , then cancerous steady state condition of the model exists as guaranteed by Theorem 4.5. Based on the Theorem 4.6, the solution of $\bar{Q}\left(\xi \right) = 1$ is unique and $\xi < 0$ . Furthermore, based on Lemma 4.7 the root obtained is the dominant root. Because we obtain negative root, then the cancerous steady state condition is locally asymptotically stable.

5. Numerical simulation

In this section, we show some illustrations of the solutions of System (3.6)–(3.9). We study the model in a constant condition. This implies that the value of all parameters is all constants. The parameter values that are used for simulations are shown in the Table 2.

Table 2. Parameters value.

Parameter	Value	Reference
$\Lambda$	0.1	^[31]
$\lambda(a)$	0.05	^[32]
$h(b)$	0.2	Assumed
$\delta (a)$	0.15	^[10]
$\gamma (a)$	0.04	^[7,10]

| Show Table

DownLoad: CSV

The interaction between the cells is usually begin shortly after the proliferation, which is in $G_1$ phase. We assume that the time average, which is needed from the proliferation to $G_{1}$ , is 1 hour. On the other hand, we assume that one cycle life span required by a cell is 100 hours. Based on this facts and assumptions, we define the age interval is $\left[1,100\right]$ . The basic reproduction number is calculated by substituting the parameter value in Table 2 into (3.27), then we obtain

$\begin{align*} R_0& = \int\limits_{1}^{100}{0.2\left[ \int\limits_{1}^{b}{{{e}^{-\int\limits_{\nu }^{b}{0.15dj}}}{{e}^{-0.1\left( b-\nu \right)}}0.05d\nu } \right.} \\ & \left. +\int\limits_{1}^{b}{{{e}^{-0.1\left( b-\nu \right)}}}0.05\int\limits_{\nu }^{b}{0.04}\int\limits_{\tau }^{\nu }{0.15{{e}^{-\int\limits_{w}^{\tau }{0.15dj}}}{{e}^{-\int\limits_{w}^{\nu }{0.04dj}}}}dwd\tau d\nu \right]db \\ & \approx 2.83. \end{align*}$

Since $R_{0} > 1$ , then based on the Theorem 4.5, the virus exists and then spreads in all parts of the tissue, see Figure 1. The density of susceptible cells (see Figure 1a) decreases significantly but then increases at a lower value than before, and goes to the constant value. The density of the infected cells also decreases. The decreasing of the infected cells can be caused by the death of infected cells or infected cells progressing to precancerous cells. Simultaneously, the development of the cancerous cells population is initially at the low speed even had a decrease moment, but then approaching near the 10th hours, it approaches the steady state solution. The behaviour of precancerous and cancer cells population (see Figure 1a) indicates that the cancer cells are uncontrollable. We can see in Figure 1b that the age distribution of precancerous cells population converges to a positive distribution as time evolves.

Figure 1. The behaviour of cells profiles.

DownLoad: Full-Size Img PowerPoint

We also consider the dynamics near the disease-free steady state solution. By using $h = 0.05$ , then we have $R_0 \approx 0.707 < 1$ . In Figure 1c, we show that the density of the susceptible cells increases while the density of infected, precancerous, and cancer cells decreases. We focused on showing the profile of the precancerous cell, exhibits that the age distribution of precancerous cells population converges to zero as time evolves. Moreover, illustrates the phase portraits of $s$ , $p$ , and $c$ , where the solutions converge to $s = 1$ when the initial value of $(s, i, p, c)$ is $(1, 0.1, 0, 0), (1, 0.2, 0.01, 0), (1, 0.3, 0.02, 0)$ . This example shows the solution of the infected, precancerous and cancer cells populations are decreasing and then vanish. All cells will be susceptible or normal (see Theorem 4.4). It is also showed in , where the initial value of $(s, i, p, c)$ is at $(0.01, 0.1, 0.01, 0.1), (0.011, 0.2, 0.02, 0.2), (0.013, 0.3, 0.03, 0.3)$ , the solutions tends to cancerous steady state solution. Thus, the disease-free steady state solution is unstable. Simultaneously, the cancerous steady state solution becomes stable, see Corollary 4.8.

Figure 2. Phaseportrait of susceptible, precancer and cancer cells with different initial value.

DownLoad: Full-Size Img PowerPoint

Furthermore, if we change the value of $\delta$ from $0.15$ to $0.165$ and substitute it to the system, the density of infected cells decreases. However, the density of precancerous and cancer cells will increase. This behaviour can be seen from , the higher the value of $\delta$ , the more the infected cells to progress to precancerous cells. Since $\delta$ means the progression rate from infected to precancerous cells, it is attempted to reduce the value of $\delta$ , so that infected cells do not have the possibility to become precancerous cells. In this case, we need a control therapy by reducing the value of $\delta$ . On the other hand, based on Theorem 4.4, if the sufficient condition hold, then by the appropriate therapy the infected, precancerous, and cancer cells can be eliminated, see ^[9].

Figure 3. The behaviour changes of cells due to changes in the values of

$\delta$ when

$R_0 > 1$ .

DownLoad: Full-Size Img PowerPoint

6. Concluding remarks

This paper proposes a new mathematical model of cervical cancer based on the age-structured of cells at the tissue level. The model is an extension of the one in ^[7], which is, by adding the age of cells to the system. We consider the patterns of the steady state solutions and their stabilities to determine the dynamics of HPV infections based on the age of cells. Our model is a sequel work of the one in ^[10] with a different scenario, that is, by adding the natural growth of the cells and the force of infections parameter. The force of infections represents the HPV infection rate to the normal cells.

The steady state condition of the system has been discussed in Lemma 3.1. Furthermore, the stability analysis of the model was presented through basic reproduction number, $R_{0}$ . If the value of $R_{0} < 1$ , then the disease-free steady state solution is globally asymptotically stable. As our analytical result, a cancerous steady state solution will exist and locally asymptotically stable based on certain conditions; see Theorem 4.5 and Corollary 4.8. By the numerical simulations, we have shown that different changes of $\delta$ will affect the cells behaviour. Parameter $\delta$ indicates the progression rate from infected to precancerous cells. If this value is lowered, then it is in line with the decreasing of precancerous cells density (see Figure 3b). This result related to ^[7] by minimizing the precancerous cells population, then it will minimize the risk of cancer cell growth.

Mathematical research constantly develops a model that addresses modern and basic cancer treatments, such as using chemotherapy ^[33,34,35], immunotherapy ^[36,37], TNF- $\alpha$ inhibitor therapy ^[38]. Moreover, it is beneficial to analyze the optimal control problem ^[39,40] and impulse control ^[41]. The possibility of metastasis conditions for cervical cancer also needs to be discussed in a mathematical point of view. However, some results that should be important are not discussed in our manuscript, such as cancer treatment, optimal control and impulse control, the metastasis condition for cervical cancer. These studies are still an open problem in this paper.

Acknowledgments

The first author would like to thank LPDP Indonesia for the Doctoral scholarship. This research is partially funded by Universitas Gadjah Mada through research scheme "Rekognisi Tugas Akhir" 2021 with the letter of assignment number 3143/UN1.P.III/D1T.L1T/PT/2021.

Conflict of interest

All authors declare that they have no conflict of interest.

References

[1]	S. V. Graham, The human papillomavirus replication cycle, and its links to cancer progression: a comprehensive review, Clin. Sci., 131 (2017), 2201–2221. https://doi.org/10.1042/CS20160786 doi: 10.1042/CS20160786
[2]	H. Sung, J. Ferlay, R. L. Siegel, M. Laversanne, I. Soerjomataram, A. Jemal, et al., Global cancer statistics 2020: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries, CA: Cancer J. Clin., 71 (2021), 209–249. https://doi.org/10.3322/caac.21660 doi: 10.3322/caac.21660
[3]	G. M. Clifford, J. S. Smith, M. Plumme, N. Muñoz, S. Franceschi, Human papillomavirus types in invasive cervical cancer worldwide: a meta-analysis, Brit. J. Cancer, 88 (2003), 63–69. https://doi.org/10.1038/sj.bjc.6600688 doi: 10.1038/sj.bjc.6600688
[4]	T. Sasagawa, H. Takagi, S. Makinoda, Immune responses against human papillomavirus (HPV) infection and evasion of host defense in cervical cancer, J. Infect. Chemother., 18 (2012), 807–815. https://doi.org/10.1007/s10156-012-0485-5 doi: 10.1007/s10156-012-0485-5
[5]	E. M. Burd, Human papillomavirus and cervical cancer, Clin. Microbiol. Rev., 16 (2003), 1–17. https://doi.org/10.1128/CMR.16.1.1-17.2003 doi: 10.1128/CMR.16.1.1-17.2003
[6]	C. A. Moody, L. A. Laimins, Human papillomavirus oncoproteins: pathways to transformation, Nat. Rev. Cancer, 10 (2010), 550–560. https://doi.org/10.1038/nrc2886 doi: 10.1038/nrc2886
[7]	T. S. N. Asih, S. Lenhart, S. Wise, L. Aryati, F. Adi-Kusumo, M. S. Hardianti, et al., The dynamics of HPV infection and cervical cancer cells, Bull. Math. Biol., 78 (2016), 4–20. https://doi.org/10.1007/s11538-015-0124-2 doi: 10.1007/s11538-015-0124-2
[8]	T. S. N. Asih, M. Masrukan, The analysis and interpretation of the all exist unstable equilibrium points of cervical cancer mathematical modeling, Proc. ICMSE, 4 (2017), 127–129.
[9]	L. Aryati, T. S. Noor-Asih, F. Adi-Kusumo, M. S. Hardianti, Global stability of the disease free equilibrium in a cervical cancer model: a chance to recover, Far East J. Math. Sci., 103 (2018), 1535–1546. https://doi.org/10.17654/MS103101535 doi: 10.17654/MS103101535
[10]	V. V. Akimenko, F. Adi-Kusumo, Stability analysis of an age-structured model of cervical cancer cells and HPV dynamics, Math. Biosci. Eng., 18 (2021), 6155–6177. https://doi.org/10.3934/mbe.2021308 doi: 10.3934/mbe.2021308
[11]	K. Allali, Stability analysis and optimal control of HPV infection model with early-stage cervical cancer, Biosystems, 199 (2021), 104321. https://doi.org/10.1016/j.biosystems.2020.104321 doi: 10.1016/j.biosystems.2020.104321
[12]	T. Malik, A. Gumel, E. Elbasha, Qualitative analysis of an age and sex structured vaccination model for human papillomavirus, Discrete Contin. Dynam. Syst. Ser. B, 18 (2013), 2151–2174. https://doi.org/10.3934/dcdsb.2013.18.2151 doi: 10.3934/dcdsb.2013.18.2151
[13]	M. Al-Arydah, R. Smith, An age-structured model of human papillomavirus vaccination, Math. Comput. Simul., 82 (2011), 629–652. https://doi.org/10.1016/j.matcom.2011.10.006 doi: 10.1016/j.matcom.2011.10.006
[14]	M. Al-Arydah, T. Malik, An age-structured model of the human papillomavirus dynamics and optimal vaccine control, Int. J. Biomath., 10 (2017), 1750083. https://doi.org/10.1142/S1793524517500838 doi: 10.1142/S1793524517500838
[15]	L. Spinelli, A. Torricelli, P. Ubezio, B. Basse, Modelling the balance between quiescence and cell death in normal and tumour cell populations, Math. Biosci., 202 (2006), 349–370. https://doi.org/10.1016/j.mbs.2006.03.016 doi: 10.1016/j.mbs.2006.03.016
[16]	Z. Liu, J. Chen, J. Pang, P. Bi, S. Ruan, Modeling and analysis of a nonlinear age-structured model for tumor cell populations with quiescence, J. Nonlinear Sci., 28 (2018), 1763–1791. https://doi.org/10.1007/s00332-018-9463-0 doi: 10.1007/s00332-018-9463-0
[17]	M. Gyllenberg, G. F. Webb, A nonlinear structured population model of tumor growth with quiescence, J. Math. Biol., 28 (1990), 671–694. https://doi.org/10.1007/BF00160231 doi: 10.1007/BF00160231
[18]	B. Basse, P. Ubezio, A generalised age-and phase-structured model of human tumour cell populations both unperturbed and exposed to a range of cancer therapies, Bull. Math. Biol., 69 (2007), 1673–1690. https://doi.org/10.1007/s11538-006-9185-6 doi: 10.1007/s11538-006-9185-6
[19]	G. S. Chaffey, D. J. Lloyd, A. C. Skeldon, N. F. Kirkby, The effect of the $G_1$ -S transition checkpoint on an age structured cell cycle model, PloS One, 9 (2014), e83477. https://doi.org/10.1371/journal.pone.0083477 doi: 10.1371/journal.pone.0083477
[20]	M. Nowak, R. M. May, Virus Dynamics: Mathematical Principles of Immunology and Virology, Oxford University Press, UK, 2000.
[21]	S. Patil, R. S. Rao, N. Amrutha, D. S. Sanketh, Analysis of human papilloma virus in oral squamous cell carcinoma using p16: An immunohistochemical study, J. Int. Soc. Prev. Community Dent., 4 (2014), 61–66. https://doi.org/10.4103/2231-0762.131269 doi: 10.4103/2231-0762.131269
[22]	J. Yang, Z. Qiu, X. Li, Global stability of an age-structured cholera model, Math. Biosci. Eng., 11 (2014), 641–665. https://doi.org/10.3934/mbe.2014.11.641 doi: 10.3934/mbe.2014.11.641
[23]	Q. Richard, Global stability in a competitive infection-age structured model, Math. Model. Nat. Phenom., 15 (2020), 54. https://doi.org/10.1051/mmnp/2020007 doi: 10.1051/mmnp/2020007
[24]	X. Rui, X. Tian, F. Zhang, Global dynamics of a tuberculosis transmission model with age of infection and incomplete treatment, Adv. Differ. Equations, 242 (2017), 1–34. https://doi.org/10.1186/s13662-017-1294-z doi: 10.1186/s13662-017-1294-z
[25]	X. Tian, R. Xu, N. Bai, J. Lin, Bifurcation analysis of an age-structured SIRI epidemic model, Math. Biosci. Eng., 17 (2020), 7130–7150. https://doi.org/10.3934/mbe.2020366 doi: 10.3934/mbe.2020366
[26]	C. M. Martin, J. J. O'Leary, Histology of cervical intraepithelial neoplasia and the role of biomarkers, Best Pract. Res. Clin. Obstet. Gynaecol., 25 (2011), 605–615. https://doi.org/10.1016/j.bpobgyn.2011.04.005 doi: 10.1016/j.bpobgyn.2011.04.005
[27]	X. Li, J. Liu, M. Martcheva, An age-structured two-strain epidemic model with super-infection, Math. Biosci. Eng., 7 (2010), 123. https://doi.org/10.3934/mbe.2010.7.123 doi: 10.3934/mbe.2010.7.123
[28]	A. Khan, G. Zaman, Global analysis of an age-structured SEIR endemic model, Chaos Solitons Fract., 108 (2018), 154–165. https://doi.org/10.1016/j.chaos.2018.01.037 doi: 10.1016/j.chaos.2018.01.037
[29]	X. Li, J. Yang, M. Martcheva, Age Structured Epidemic Modeling, Springer Nature, Switzerland, 2020.
[30]	H. Inaba, Threshold and stability results for an age-structured epidemic model, J. Math. Biol., 28 (1990), 411–34. https://doi.org/10.1007/BF00178326 doi: 10.1007/BF00178326
[31]	A. K. Miller, K. Munger, F. R. Adler, A mathematical model of cell cycle dysregulation due to Human Papillomavirus infection, Bull. Math. Biol., 79 (2017), 1564–1585. https://doi.org/10.1007/s11538-017-0299-9 doi: 10.1007/s11538-017-0299-9
[32]	S. Park, S. Chung, K. M. Kim, K. C. Jung, C. Park, E. R. Hahm, et al., Determination of binding constant of transcription factor myc–max/max–max and E-box DNA: the effect of inhibitors on the binding, Biochim. Biophys. Acta, Gen. Subj., 1670 (2004), 217–228. https://doi.org/ 10.1016/j.bbagen.2003.12.007 doi: 10.1016/j.bbagen.2003.12.007
[33]	F. Ansarizadeh, M. Singh, D. Richards, Modelling of tumor cells regression in response to chemotherapeutic treatment, Appl. Math. Modell., 48 (2017), 96–112. https://doi.org/10.1016/j.apm.2017.03.045 doi: 10.1016/j.apm.2017.03.045
[34]	F. J. Solis, S. E. Delgadillo, Evolution of a mathematical model of an aggressive–invasive cancer under chemotherapy, Comput. Math. Appl., 69 (2015), 545–558. https://doi.org/10.1016/j.camwa.2015.01.013 doi: 10.1016/j.camwa.2015.01.013
[35]	E. R. Sari, D. Lestari, E. Yulianti, R. Subekti, Stability analysis of a mathematical model of tumor with chemotherapy, J. Phys. Conf. Ser., 1321 (2019), 022072. https://doi.org/10.1088/1742-6596/1321/2/022072 doi: 10.1088/1742-6596/1321/2/022072
[36]	R. Eskander, K. S. Tewari, Immunotherapy: an evolving paradigm in the treatment of advanced cervical cancer, Clin. Ther., 37 (2015), 20–38. https://doi.org/10.1016/j.clinthera.2014.11.010 doi: 10.1016/j.clinthera.2014.11.010
[37]	P. K. Roy, A. K. Roy, E. N. Khailov, F. Al Basir, E. V. Grigorieva, Model of the optimal immunotherapy of psoriasis by introducing IL-10 and IL-22 inhibitor, J. Biol. Syst., 28 (2020), 609–639. https://doi.org/10.1142/S0218339020500084 doi: 10.1142/S0218339020500084
[38]	A. K. Roy, F. Al Basir, P. K. Roy, A vivid cytokines interaction model on psoriasis with the effect of impulse biologic (TNF- $\alpha$ inhibitor) therapy, J. Theor. Biol., 474 (2019), 63–77. https://doi.org/10.1016/j.jtbi.2019.04.007 doi: 10.1016/j.jtbi.2019.04.007
[39]	A. Khan, G. Zaman, R. Ullah, N. Naveed, Optimal control strategies for a heroin epidemic model with age-dependent susceptibility and recovery-age, AIMS Math., 6 (2021), 1377–1394. https://doi.org/10.3934/math.2021086 doi: 10.3934/math.2021086
[40]	A. K. Roy, M. Nelson, P. K. Roy, A control-based mathematical study on psoriasis dynamics with special emphasis on IL-21 and IFN- $\gamma$ interaction network, Math. Methods Appl. Sci., 44 (2021), 13403–13420. https://doi.org/10.1002/mma.7635 doi: 10.1002/mma.7635
[41]	A. K. Roy, P. K. Roy, E. Grigorieva, Mathematical insights on psoriasis regulation: Role of $\text {Th}_1$ and $\text {Th}_2$ cells, Math. Biosci. Eng., 15 (2018), 717–738. https://doi.org/10.3934/mbe.2018032 doi: 10.3934/mbe.2018032

This article has been cited by:

1.	Juan Carlos Sierra-Rojas, Ramón Reyes-Carreto, Cruz Vargas-De-León, Jorge Fernando Camacho, Khalid Hattaf, Modeling and Mathematical Analysis of the Dynamics of HPV in Cervical Epithelial Cells: Transient, Acute, Latency, and Chronic Infections, 2022, 2022, 1748-6718, 1, 10.1155/2022/8650071
2.	Eminugroho Ratna Sari, Lina Aryati, Fajar Adi-Kusumo, An age-structured SIPC model of cervical cancer with immunotherapy, 2024, 9, 2473-6988, 14075, 10.3934/math.2024685
3.	Bushra Bajjah, Mahmut Modanli, Francisco R. Villatoro, Finite Difference Method for Infection Model of HPV with Cervical Cancer under Caputo Operator, 2024, 2024, 1607-887X, 1, 10.1155/2024/2580745
4.	Shayidan Abuduwaili, Lei Wang, Zhidong Teng, Abidan Ailawaer, Ramziya Rifhat, Estimating real-time reproduction number for HPV infection in Xinjiang, China, 2025, 2025, 2731-4235, 10.1186/s13662-025-03896-x

Reader Comments

Your name:*

Email:*
© 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

Mathematical Biosciences and Engineering

3.9

Metrics

Article views(2750) PDF downloads(189) Cited by(4)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(3) / Tables(2)

Mathematical Biosciences and Engineering

Mathematical analysis of a SIPC age-structured model of cervical cancer

Related Papers:

Abstract

1. Introduction

2. Model formulation

3. Steady state conditions

4. Disease-free and cancerous steady state solutions

4.1. Local and global stability analysis of disease-free steady state solutions

4.2. Existence of cancerous steady state solution

4.3. Local stability of cancerous steady state solution

5. Numerical simulation

6. Concluding remarks

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Catalog

Mathematical Biosciences and Engineering

Mathematical analysis of a SIPC age-structured model of cervical cancer

Related Papers:

Abstract

1. Introduction

2. Model formulation

3. Steady state conditions

4. Disease-free and cancerous steady state solutions

4.1. Local and global stability analysis of disease-free steady state solutions

4.2. Existence of cancerous steady state solution

4.3. Local stability of cancerous steady state solution

5. Numerical simulation

6. Concluding remarks

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog