www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Age-structured delayed SIPCV epidemic model of HPV and cervical cancer cells dynamics I. Numerical method Vitalii V. Akimenko∗, Fajar Adi-Kusumo† ∗Institute of Biology and Medicine,, Taras Shevchenko National University of Kyiv,Ukraine vitaliiakm@gmail.com †Faculty of Mathematics and Natural Sciences Universitas Gadjah Mada, Yogyakarta, Indonesia. f adikusumo@ugm.ac.id Received: 29 April 2021, accepted: 2 October 2021, published: 6 December 2021 Abstract— The numerical method for simulation dynamics of nonlinear epidemic model of age- structured sub-populations of susceptible, infectious, precancerous and cancer cells and unstructured population of human papilloma virus (HPV) is developed (SIPCV model). Cell population dynamics is described by the initial-boundary value problem for the delayed semi-linear hyperbolic equations with age- and time-dependent coefficients and HPV dynamics is described by the initial problem for nonlinear delayed ODE. The model considers two time-delay parameters: the time between viral entry into a target susceptible cell and the production of new virus particles, and duration of the first stage of delayed immune response to HPV population growing. Using the method of characteristics and method of steps we obtain the exact solution of the SIPCV epidemic model in the form of explicit recurrent formulae. The numerical method designed for this solution and used the trapezoidal rule for integrals in recurrent formulae has a second order of accuracy. Numerical experiments with vanished mesh spacing illustrate the second order of accuracy of numerical solution with respect to the benchmark solution and show the dynamical regimes of cell- HPV population with the different phase portraits. Keywords- SIPCV epidemic model; Age- structured model; HPV. Numerical epidemiology; Method of characteristics; I. INTRODUCTION Human papilloma virus (HPV) is the most common sexually transmitted infection that can cause dysplasia (precancerous tissue) and cervical cancer [27], [36], [45]. The importance of this problem in medicine motivated extensive studied of HPV-epidemic models over the last decades. The problem is considered at two levels: (i) social level (trans-mission of HPV between people and Copyright: © 2021 Akimenko et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: Vitalii V. Akimenko, Fajar Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical cancer cells dynamics I. Numerical method, Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 1 of 23 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... effectiveness of HPV vaccination of population) [13], [43], and (ii) molecular or tissue level (cell- HPV population dynamics) [17], [48], [54]. Early epidemic models of the second level (ii) are based on the systems of nonlinear ODE and consider the dynamics of several compartments-subclasses of cells and virus. Since viruses are non-living things which do not replicate their population dynamics can be efficiently described by the unstructured models. But such models are not efficient enough for the cell populations because they neglect the life history of biological cell as living organ-ism and provide only restricted description of system in many applications. The need for integration of cell life-histories with population dynamics moti- vates the development of age-structured, and more generally, physiologically structured, population models in cell population dynamics and tumor growth modelling [18], [21], [29], [32], [35], [38], [39], [40], [41], [42], [44], [47], [49], [50]. The model considered in this paper is based on the epidemic models studied in [3], [48], and includes the age-structured models of susceptible, infective, precancerous, cancer cells populations and unstructured model of human papilloma virus population (SIPCV epidemic model). We use the L. Hayflick limit theory [20] for modelling prolif- eration in all cell subpopulations. Biological cell assumed to divide into mother and daughter cells with different developmental potential [23]. The model considers life history of mother cells: birth, maturing up to the age when they can proliferate, limited number of divisions at the reproductive age when it gives birth to several daughter cells, aging up to the final reproductive age and death. The features of new model are: (i) death rates of infected, precancerous and cancerous cells do not depend from the HPV abundance since the immune response of biological organism is tol- erant with respect to its own cells [27], [30], [36], [45]; (ii) death rate of HPV is density- dependent function due to the immune response on the virus population growing [27], [36], [45]; (iii) interaction strength between susceptible and HPV is a product of the Lotka-Voltera incidence rate and result in the growth of infective cells [3]; (iv) infective cells partially move to the precancerous subclass and partially apoptose when viruses leave infectious cells and ready to infect new susceptible cells [31], [48]; (v) precancerous cells move to the cancer subclass with the non-linear density- dependent saturated rate [48]; (vi) two time-delay parameters describe the time between viral entry into a target susceptible cell and the production of new virus particles [26] and a duration of the first stage of delayed immune response to HPV population growing [27], [36], [45]. Development of population dynamics models of mathematical epidemiology necessarily leads to the development of methods of numerical epidemiology in simula- tion of cell-virus interaction dynamics. There are two main approaches to numerical simulation of physiologically structured models of population dynamics. The first one uses the difference-finite or element-finite approximation of the boundary-initial value problem for the trans- port equation [1], [7], [10], [19], [34], [44], [51], [52]. This approach is preferred when solving the problems with complex equations, complex domain of definition or/and non-local boundary conditions for which it is difficult or impossible to obtain an exact solution. On the other hand, only biologically correct monotone and conservative difference schemes are applicable for the transport and reaction-diffusion equations in practice. Since in work [33] it was proved that the physically correct monotone linear difference scheme of the second order of approximation and higher for the transport or reaction-diffusion equations does not exist, one has to design of the new special non- linear monotone difference schemes of the second order of approximation and higher for these types of equations [7], [10], [19]. The second approach is based on method of characteristics of the theory of hyperbolic equations of the first order which reduces the PDE - transport equation to the ODE of the first order and allows for applying the corresponding well-known numerical methods to this equation. This approach includes the escalator boxcar train (EBT) method [25], method of char- Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 2 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... acteristics [1], [3], [5], and different combinations and variations of them [2], [14], [15], [16], [37], [52], [56]. Numerical methods of this type provide the biologically correct numerical solution of the second and higher order of approximation in most cases. But despite the importance of study and development of numerical epidemiology nowadays there are no works with numerical implemen- tations of such methods for the physiologically structured epidemic models so far. Using the method of characteristics [3], [4], [5], [6], [8], [9], [11], [24], [53], [55], [58] and method of steps from the theory of delayed differential equations [12], [28], [46], [57], [59], we obtain an exact solution of the SIPCV epidemic model. This solution is given in form of the recurrent formulae (like in works [3], [4], [8], [55]) in which the densities of all subpopulations are defined through the integrals from solution taken at previous in- stance of time. For this exact solution we create the numerical method of the second order of accuracy using the trapezoidal rule [22] for approximation of all integrals in the recur-rent formulae. The ex- act solution of the age-structured delayed SIPCV epidemic system provides us the advantages in developing of the more fast and accurate numerical method in comparison with the other methods used for the age- and size- structured models [1], [5], [25]. Evaluation of the convergence rate of the approximate solution to the exact one is considered in the next paper, the second part of our work. For illustration of accuracy of the designed numerical method we evaluate in experiments the residual of deviation of numerical solution from the benchmark solution of SIPCV model obtained for the special particular coefficients of equations and initial values. Simulations show that the rela- tive numerical error converges pointwise to zero with h → 0 (where h is a mesh spacing) by the quadratic low that illustrates and confirms the second order of accuracy of numerical solution. These results are in good agreement with the evaluations of numerical errors obtained earlier in experiments [8] with the numerical method designed for the nonlinear age-structured models of population dynamics by the same approach. Numerical experiments with model parameters of the system reveal two types of the asymptotically stable dynamical regimes-non-oscillating and os- cillating dynamics of population when the quan- tity of cells and HPV converge to the non-trivial asymptotically stable states of the system. These dynamical regimes correspond to the localization of dysplasia (precancer cells) and cancer tumor in biological tissue without metastases. Overall, the numerical method obtained in this paper provides the reliable and accurate theoretical instrument for simulation and study of age-structured SIPCV epidemic models. II. MODEL We consider a SIPCV epidemic model that consists of susceptible (noninfected), in-fectious (without significant changing of morphology, CIN I and CIN II stages [30], [36]), precancerous (with changed by virus morphology - dysplasia, but is differentiable yet, CIN III stage [30], [36]), cancer (nondifferentiable) cells and human papilloma virus (HPV) that moves freely between cells. The age-specific densities of susceptible, infectious, pre-cancerous and cancer cells are denoted as S(a,t), I(a,t), P(a,t) and C(a,t). The dynamics of cell subclasses (subpopulations) is described by the nonlinear age-structured model with age- and time-dependent death rates of susceptible ds(a,t), infectious dq(a,t), precancerous dr(a,t) and cancer cells dc(a,t) with maximum lifespan ad, an age reproductive window of non-cancer cells [ar,am] and cancer cells [ac,ag], ac < ar, ag < am. We assume that adaptive behavior of the HPV makes the immune response of biological organism (both T-killers cells and humoral immunity) tolerant with respect to infectious and precancerous cells and noneffective with respect to cancer cells that is their death rates does not depend from the HPV abundance [27], [30], [36]. Organism recognizes cancer cells as its own and does not at-tempt to destroy them. Since viruses are not living things and cannot reproduce (multiply) until they enter a living cell, we use the ODE model Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 3 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... instead of age-structured model for describing the dynamics of HPV quantity V (t). The model considers the time-dependent recruitment rate of viruses Λ(t), their density-dependent death rate dv(V (t − σ)) with delay parameter σ - the duration of the first stage of delayed immune response to HPV population growing [36]. The interaction strength be-tween susceptible and HPV is a product of the Lotka-Voltera incidence rate α(a − θ,t − θ)V (t − θ)S(a,t) (α(a,t) is a time and age dependent infection rate, θ is a delay parameter (θ < σ) - time between viral entry into a target cell and the production of new virus particles [26]), and result in the growth of infectious cells, with partial move of them to the precancerous subclass with rate δ(a,t)I(a,t) (δ(a,t) is a time and age dependent progression rate from infectious to precancerous cells - dysplasia). Model considers also the partial apoptosis of infectious and precancerous cells with rate n(t) ad∫ 0 (d∗q(a,t)I(a,t) + d ∗ r(a,t)P(a,t))da, (where n(t) is a time dependent mean number of virions produced by one cell, d∗q(a,t) and d∗r(a,t) are time and age dependent death rates of infectious and precancerous cells as a result of virus replication), when viruses leave destroyed cells and ready to infect new susceptible cells. We assume that 0 < d∗q(a,t) < dq(a,t), 0 < d∗r(a,t) < dr(a,t), that is some of infectious and precancerous cells may heal themselves or, at least, can slow down the replication of the HPV within themselves and die when they reach the maximum age or as a result of exposure to some external factors. Cancer cells cannot be defined in absolute terms and usually they are recognized only by the ab- normal rapid proliferation (ac < ar, ag < am) when cells do not have time for maturation (”non- differentiable” cells). Cancer cells differ from the normal ones in their lack of response to nor- mal fertility control mechanism [27], [30], [36], [45]. That is why precancerous cells turn to the cancer in our model only through the prolifer- ation when a precancerous cell may divide into two new cancer ones. Fertility rate of precancer- ous cells turned to the cancer is proportional to the density-dependent saturated rate µ(Nr(t)) = ρNr(t) 1+wNr(t) [48], where the quantity of precancerous cells of re-productive age and older is Nr(t) = ad∫ ar P(a,t)da, ρ ∈ (0, 1] is a progression rate from precancerous to cancerous cells, w ≥ 1 is a coefficient of satura-tion, 0 ≤ µ(Nr) < ρw−1 for Nr ≥ 0. These assump-tions lead to the fol- lowing age-structured epidemic model in domain Q = {(a,t) |a ∈ (0,ad), t ∈ (0,T)}: ∂S(a,t) ∂t + ∂S(a,t) ∂a = −ds(a,t)S(a,t) −α(a−θ,t−θ)V (t−θ)S(a,t), (1) ∂I(a,t) ∂t + ∂I(a,t) ∂a =−(dq(a,t)+δ(a,t))I(a,t) + α(a−θ,t−θ)V (t−θ)S(a,t), (2) ∂P(a,t) ∂t + ∂P(a,t) ∂a = −dr(a,t)P(a,t) + δ(a,t)I(a,t), (3) ∂C(a,t) ∂t + ∂C(a,t) ∂a = −dc(a,t)C(a,t), (4) ∂V (t) ∂t = Λ(t) −dv(V (t−σ))V (t) + n(t) ad∫ 0 (d∗q(a,t)I(a,t)+d ∗ r(a,t)P(a,t))da (5) Equations (1) - (5) are completed by the boundary conditions and initial values: S(0, t) = am∫ ar βs(a,t)S(a,t)da (6) I(0, t) = am∫ ar βq(a,t)I(a,t)da (7) P(0, t) = (1−µ(Nr(t)) am∫ ar βr(a,t)P(a,t)da (8) Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 4 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... C(0, t) = ag∫ ac βc(a,t)C(a,t)da + µ(Nr(t)) am∫ ar βr(a,t)P(a,t)da (9) S(a, 0) =ϕ(a),I(a,0) =P(a,0) =C(a,0) = 0 (10) V (t) = V0(t), t ∈ [−σ, 0], (11) where the birth (fertility) rates of the suscepti- ble, infectious, precancerous and cancer cells are βs(a,t), βq(a,t), βp(a,t) and βc(a,t), respec- tively; ϕ(a) is an initial density of susceptible cells, V0(t) is an initial value of HPV quantity. Function α(a−θ,t−θ) is a rate of infection which takes into account infection of susceptible cells over the age θ (survived after incubation period) at the unit of time: α(a−θ,t−θ) = { 0, if a < θ α0(a,t−θ), if a ≥ θ, (12) where α0(a,t) is an auxiliary rate of infection defined for a ∈ [θ,ad], t ∈ [−θ,T−θ]. We impose the following restrictions on the density-dependent HPV death rate and cells’ birth and death rates [27], [30], [31], [36], [45]: δ(a,t),dc(a,t),dr(a,t),dq(a,t),ds(a,t)>0, (13) βs(a,t),βq(a,t),βr(a,t),βc(a,t) > 0, (14) α0(a,t) ≥ 0, Λ(t) ≥ 0,n(t) > 0, (15) dv(V ) > 0, ∂dv(V ) dV ≥ 0 for V > 0, (16) ϕ(a) ≥ 0, ad∫ 0 ϕ(a)da > 0,V0(t) ≥ 0 (17) The positiveness of derivative of density- dependent death rate of HPV in (15) means that increasing of HPV quantity changes the charac- teristics of intracellular space that result in the organism immune response through the activation of cell immunity (T-killers) and humoral immunity (B-lymphocytes) that leads to the elimination of viruses (i.e. monotone increasing of their death rate). We assume that all coefficients and ini- tial values of system (1)-(11) are twice continu- ously differentiable functions and have the private derivatives of the second order by all their argu- ments. III. SUSCEPTIBLE CELL POPULATION DYNAMICS Since the time delay θ in equation (1) provides the formal linearization, we can apply the method of steps [3], [28], from the theory of linear delayed ODE (1), (6), (10) and the method of characteristic [3], [4], [8], [21], [53], [55] from the theory of line-ar hyperbolic equations of the first order to the initial-boundary problem for the quasi-linear transport equation. Without loss of generality time cut is covered by the set of consequent time periods [tk−1, tk] (tk = kad, k = 1, ...,K, t0 = 0, tK = T), and we define the following sets (Fig. 1): Q (1) k = {(a,t)|t ∈ [(k − 1)ad,a + (k − 1)ad), a ∈ [0,ad]} (18) Q (2 ) k = {(a,t)|t ∈ [a + (k − 1)ad,kad], a ∈ [0,ad]} where Q = K⋃ k=1 ( Q (1) k ⋃ Q (2) k ) . We define also the auxiliary set of age intervals: Ω(k) = {[−a(k)l ,−a (k) l+1]|a (k) l = lar + (k − 1)ad, l = 0, ...,L− 1,a(k)L = am + (k − 1)ad, a (k) L+1 = kad} (19) L = { [am/ar] + 1, if am/ar − [am/ar] > 0, am/ar, if am/ar − [am/ar] = 0, (20) where k = 1, ...,K, and [a] - is an integer part of real number a. Using new characteristic variable v = a − t, and time t we reduce the problem (1), (6), (10) to Cauchy problem for the linear homogeneous delayed ODE: ∂S ∂t = −ds(v + t,t)S(v + t,t) −α(v + t−θ,t−θ)V (t−θ)S(v + t,t) (21) S(v, 0) = ϕ(v). (22) Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 5 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... Fig. 1. Splitting of domain Q̄ and the age intervals[ a (1) l ,a (1) l+1 ] , for l = 0, ...,L. Renewal boundary condition (6) completes the problem (19), (20). In original variables, solution of problem (21), (22), has a form (k = 1, ...,K): S(a,t) =   F (k−1) 1 (a−t)W1(a−t,tk−1, t), if (a,t) ∈ Q(1 )k F (k) 1 (a−t)W1(a−t,t−a,t), if (a,t) ∈ Q(2 )k , (23) W1(v,tk,t) =e − t∫ tk (ds(v+ξ,ξ)+α(v+ξ−θ,ξ−θ)V (ξ−θ))dξ (24) Functions F(k)1 (a−t) are defined from the initial (10) and boundary (6) conditions: F (0) 1 (v) = ϕ(v),v ∈ [0,ad] (25) F (k) 1 (v) = Φ (k) 1l (v),v ∈ [ −a(k)l ,−a (k) l−1 ] , (26) l = 1, ...,L + 1, where the set of auxiliary functions Φ(k)1l (v) is defined as: Φ (k) 11 (u) = am+u∫ ar+u βs(v −u,−u) F (k−1) 1 (v) W1(v,tk−1,−u)dv (27) u ∈ [ −a(k)1 ,−a (k) 0 ] , Φ (k) 1l (u) = −a(k)l−2∫ ar+u βs(v−u,−u)Φ (k) 1(l−1)(v)W1(v,−v,−u)dv + l−3∑ j=0 −a(k)j∫ −a(k)j+1 βs(v−u,−u)Φ (k) 1(j+1) (v)W1(v,−v,−u)dv + am+u∫ −a(k)0 βs(v−u,−u)F (k−1) 1 (v)W1(v,tk−1,−u)dv, (28) u ∈ [ −a(k)l ,−a (k) l−1 ] , l = 2, ...,L, Φ (k) 1(L+1) (u) = am+u∫ ar+u βs(v −u,−u) F (k) 1 (v) W1(v,−v,−u)dv (29) u ∈ [ −a(k)L+1,−a (k) L ] . Two parts of the solution (23) have to be linked continuously, that is S(a,t) ∈ C(Q) at the points of characteristics a = t − tk−1 in directions a = const, . By analogy with [3], [4], [8], we can write the compatibility (continuity) condition of solution S(a,t) ∈ C(Q) in the form: ϕ(0) = am∫ ar βs(a, 0)ϕ(a)da (30) where the second term in the right side of Eq. (28) should be omitted for l = 2, function F(k)1 (v) in Eq. (29) has been already defined on the previous steps because its argument v ∈ [ar + u,am + u] Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 6 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... and v > u. The schematic diagram of definition of the functional sequence Φ(1)1l (v) for the example L = 4 is given in Fig. 2. IV. INFECTIOUS CELL POPULATION DYNAMICS Using the same approach as in the previous section we reduce the problem (2), (7), (10), to the Cauchy problem for linear nonhomogeneous delayed ODE: ∂I ∂t =−d̃q|(v+t,t)I(v+t,t) (31) + α(v+t−θ,t−θ)V (t−θ))S(v+t,t) I(v, 0) = 0, (32) where d̃q(v,t) = dq(v,t) + δ(v,t), function V (t) is taken from the previous instant t − θ, func- tion S(v + t,t) is taken from the previous sec- tion. Solution of problem (31), (32) has a form (k = 1, ...,K): I(a,t) =   F (k−1) 2 (a− t)W2(a− t,tk−1, t) +Z (k−1) 2 (a− t,tk−1, t), if (a,t) ∈ Q(1)k , F (k) 2 (a− t)W2(a− t,t−a,t) +Z (k) 2 (a− t,t−a,t), if (a,t) ∈ Q(2)k , (33) W2(v,tk, t) = exp  − t∫ tk d̃q(v + ξ,ξ)dξ   , (34) Z (k) 2 (v,tk, t) = t∫ tk W2(v,ξ,t)α(v + ξ −θ,ξ −θ) V (ξ −θ)S(v + ξ,ξ)dξ. (35) Functions F(k)2 (v) are defined from the initial (10), and boundary (7) conditions: F (0) 2 (v) = 0, v ∈ [0,ad], (36) F (k) 2 (v) = Φ (k) 2l (v), v ∈ [ −a(k)l ,−a (k) l−1 ] , (37) l = 1, ...,L + 1, where the set of auxiliary functions Φ(k)2l (v) is defined as: Φ (k) 21 (u)= am+u∫ ar+u βq(v−u,−u)(F (k−1) 2 (v)W2(v,tk−1,−u) + Z (k−1) 2 (v,tk−1,−u))dv, (38) u ∈ [ −a(k)1 ,−a (k) 0 ] , Φ (k) 2l (u)= −a(k)l−2∫ ar+u βq(v−u,−u)(Φ (k) 2(l−1)(v)W2(v,−v,−u) + Z (k) 2 (v,−v,−u))dv + l−3∑ j=0 −a(k)j∫ −a(k)j+1 βq(v −u,−u)(Φ (k) 2(j+1) (v) W2(v,−v,−u) + Z (k) 2 (v,−v,−u))dv + am+u∫ −a(k)0 βq(v −u,−u)(F (k−1) 2 (v) W2(v,tk−1,−u)+Z (k−1) 2 (v,tk−1,−u))dv, (39) u ∈ [ −a(k)l ,−a (k) l−1 ] , l = 2, ...,L, Φ (k) 2(L+1) (u)= am+u∫ ar+u βq(v−u,−u)(F (k) 2 (v)W2(v,−v,−u) + Z (k) 2 (v,−v,−u))dv (40) u ∈ [ −a(k)L+1,−a (k) L ] , where the second term in the right side of Equa- tion (39) should be omitted for l = 2, func- tion F(k)2 (v) in Eq. (40) has been already de- fined on the previous steps be-cause its argument v ∈ [ar + u,am + u] and v > u (see exam- ple in Fig.2). Since I(a,t) has the trivial initial value, two parts of the solution (33) satisfy the compatibility (continuity) condition of solution, I(a,t) ∈ C(Q). Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 7 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... Fig. 2. Schematic diagram of definition of the sequence Φ(1)1l (v), l = 1, ...,L + 1, L = 4. V. PRECANCEROUS CELL POPULATION DYNAMICS The problem (3), (8), (10) is reduced to the Cauchy problem for nonlinear non-homogeneous ODE ∂P ∂t = −dr(v + t,t)P(v + t,t) + δ(v + t,t)I(v + t,t) (41) P(v, 0) = 0, (42) where I(v+t,t) is taken from the previous section. Function µ(v+t−σ,NP (t)) will be defined below in the recurrent formulae for the travelling wave solution. Solution of problem (41), (42) has a form (k = 1, ...,K): P(a,t) = P(v,t) =   F (k−1) 3 (a− t)W3(a− t,tk−1, t) +Z (k−1) 3 (a− t,tk−1, t), if (a,t) ∈ Q(1 )k , F (k) 3 (a− t)W3(a− t,t−a,t) +Z (k) 3 (a− t,t−a,t), if (a,t) ∈ Q(2 )k , , (43) W3(v,tk, t) = exp  − t∫ tk dr(v + ξ,ξ)dξ   , (44) Z (k) 3 (v,tk−1, t) = t∫ tk−1 W3(v,ξ,t)δ(v + ξ,ξ) I(v + ξ,ξ)dξ (45) Functions F(k)3 (v) and NP (ξ) are defined from the initial (10) and boundary (8) conditions: F (0) 3 (v) = 0,v ∈ [0,ad] (46) F (k) 3 (v) = Φ (k) 3l (v),v ∈ [ −a(k)l ,−a (k) l−1 ] , (47) l = 1, ...,L + 1, where the auxiliary functions Φ(k)3l (v) are defined as: Φ (k) 31 (u) = (1 −µ(Nr(−u))) am+u∫ ar+u βr(v −u,−u) ( F (k−1) 3 (v)W3(v,tk−1,−u) +Z (k−1) 3 (v,tk−1,−u) ) dv, (48) u ∈ [ −a(k)1 ,−a (k) 0 ] Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 8 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... Np(−u) = ad+u∫ ar+u ( F (k−1) 3 (v)W3(v,tk−1,−u) +Z (k−1) 3 (v,tk−1,−u) ) dv,u∈[−a(k)1 ,−a (k) 0 ], (49) Φ (k) 3l (u) = (1 −µ(Nr(−u)))( −a(k)l−2∫ ar+u βr(v−u,−u) ( Φ (k) 3(l−1)(v)W3(v,−v,−u) ) +Z (k) 3 (v,−v,−u) ) dv l−3∑ j=0 −a(k)j∫ −a(k)j+1 βr(v −u,−u) ( Φ (k) 3(j+1) (v) W3(v,−v,−u) +Z (k) 3 (v,−v,−u) ) dv + am+u∫ −a(k)0 βr(v −u,−u) ( F (k−1) 3 (v) W3(v,tk−1,−u) +Z (k−1) 3 (v,tk−1,−u) ) dv ) u ∈ [ −a(k)l ,−a (k) l−1 ] , l = 2, ...,L (50) Np(−u) = −a(k)l−2∫ ar+u ( Φ (k) 3(l−1)(v)W3(v,−v,−u) +Z (k) 3 (v,−v,−u) ) dv + l−3∑ j=0 −a(k)j∫ −a(k)j+1 ( Φ (k) 3(j+1) (v)W3(v,−v,−u) +Z (k) 3 (v,−v,−u) ) dv + ad+u∫ −a(k)0 ( F (k−1) 3 (v)W3(v,tk−1,−u) +Z (k−1) 3 (v,tk−1,−u) ) dv, (51) u ∈ [ −a(k)l ,−a (k) l−1 ] , l = 2, ...,L. Φ (k) 3(L+1) (u) = 2(1 −µ(Np(−u))) am+u∫ ar+u βr(v−u,−u)F (k) 3 (v)W3(v,−v,−u) +Z (k) 3 (v,−v,−u) ) dv, (52) u ∈ [ −a(k)L+1,−a (k) L ] , Np(−u) = ad+u∫ ar+u ( F (k) 3 (v)W3(v,−v,−u) +Z (k) 3 (v,−v,−u) ) dv (53) u ∈ [ −a(k)L+1,−a (k) L ] , where the second term in the right side of Eqs. (50), (51) should be omitted for l = 2, function F (k) 3 (v) in Eqs. (52), (53) has been already de- fined on the previous steps because its argument v ∈ [ar + u,am + u] and v > u (see example in Fig.2). Because the initial value of P(a,t) is triv- ial, solution (43) satisfies the compatibility (conti- nuity) condition of solution, P(a,t) ∈ C(Q). VI. CANCER CELL POPULATION DYNAMICS Since cancer cells have a specific age reproduc- tive window [ac,ak] differed from the one of non- cancer cells, we introduce a new auxiliary set of age intervals for the problem (4), (9), (10): Ω̃(k) = {[ −ã(k)l ,−ã (k) l+1 ]∣∣∣ã(k)l = lac + (k − 1)ad, l = 0, ...,L− 1, ã(k)L = ag + (k − 1)ad, ã (k) L+1 = kad} (54) L̃= { [ag/ac]+1, if ag/ac−[ag/ac]>0, ag/ac, if ag/ac − [ag/ac] = 0, (55) Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 9 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... where k = 1, ...,K, and [a] - is an integer part of real number a. In new characteristic variable v = a−t = const problem (4), (9), (10) is reduced to the Cauchy problem for linear homogeneous ODE: ∂C ∂t = −dc(v + t,t) C(v + t,t) (56) C(v, 0) = 0, (57) where function µ(NP (t − σ)) is taken from the previous instant of time, P(a,t) is taken from the previous section. Solution of problem (56), (57) has a form (k = 1, ...,K): (a,t) = (v,t) (58) =   F (k−1) 4 (a− t)W4(a−t,tk−1, t), if (a,t) ∈ Q(1 )k , F (k) 4 (a−t)W4(a−t,t−a,t), if (a,t) ∈ Q(2 )k , (59) W4(v,tk, t) = exp  − t∫ tk dc(v + ξ,ξ)dξ   . (60) Functions F(k)4 (v) are defined from the initial (10) and boundary (9) conditions: F (0) 4 (v) = 0, v ∈ [0,ad] (61) F (k) 4 (v) = Φ (k) 4l (v), v ∈ [ −ã(k)l ,−ã (k) l−1 ] , l = 1, ..., L̃ + 1, (62) where the auxiliary functions Φ(k)4l (v) are defined as: Φ (k) 41 (u) = ag+u∫ a+u βc(v −u,−u)F (k−1) 4 (v) W4(v,tk−1,−u)dv + µ(Np(−u)) × am+u∫ ar+u βp(v −u,−u)P(v −u,−u)dv, (63) u ∈ [ −a(k)1 ,−a (k) 0 ] , Φ (k) 4l (u) = −ã(k)l−2∫ a+u β(v −u,−u)Φ(k) 4(l−1)(v) W4(v,−v,−u)dv+ l−3∑ j=0 −ã(k)j∫ −ã(k)j+1 β(v−u,−u) × Φ(k) 4(j+1) (v)W4(v,−v,−u)dv + ag+u∫ −ã(k)0 βc(v −u,−u)F (k−1) 4 (v) W4(v,tk−1,−u)dv + µ(Np(−u)) am+u∫ ar+u βp(v −u,−u)P(v −u,−u)dv, (64) u ∈ [ −ã(k)l ,−ã (k) l−1 ] , l = 2, ..., L̃, Φ (k) 4(L̃+1) (u) = ag+u∫ ac+u βc(v −u,−u) F (k) 4 (v) W4(v,−v,−u)dv + µ(Np(−u)) × am+u∫ ar+u βp(v−u,−u)P(v−u,−u)dv, (65) u ∈ [ −ã(k) L̃+1 ,−ã(k) L̃ ] , where the second term in the right side of Equa- tion (62) should be omitted for l = 2, function F (k) 4 (v) in Eq. (63) has been already defined on the previous steps be-cause its argument v ∈ [ac + u,ag + u] and v > u (see example in Fig.2). Because the initial value of C(a,t) is trivial, so- lution (57) satis-fies the compatibility (continuity) condition of solution, C(a,t) ∈ C(Q). VII. HPV POPULATION DYNAMICS Nonlinear death rate dv(V ) of Cauchy problem for nonlinear delayed ODE (5), (11), satisfies Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 10 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... restrictions (16) and, hence, it is a locally Lipschitz continuous function. Hence, the unique smooth solution of problem (5), (11) exists and can be ob-tained by the method of steps [28]: V (t) = V (0) exp  − t∫ 0 dv(V (ξ −σ))dξ   + t∫ 0 exp  − t∫ η dv(V (ξ −σ))dξ  (Λ(η) + n(η) ad∫ 0 (d∗q(a,η)I(a−η,η) + d ∗ p(a,η) P(a−η,η))da) dη, (66) where functions I(a,t) and P(a,t) are taken from the previous sections. Overall, the results obtained above are summarized in the Theorem 1. Theorem 1. Let the coefficients of SIPCV epi- demic system (1)-(11) satisfy conditions (12)-(17) and initial value ϕ(a) satisfies the compatibility condition (30), then the unique continuous solution of problem (1)-(11), P(a,t) ∈ C(Q̄), C(a,t) ∈ C(Q̄), V (t) ∈ C([0,T]) exists and can be ob- tained by explicit recurrent formulae (23)-(29), (33)-(40), (43)-(53), (58)-(65),, (66), respectively. VIII. NUMERICAL IMPLEMENTATIONS The numerical age-specific densities Sji , I j i , P j i , C j i of susceptible, infectious, precancerous, cancer cells, respectively, are de-fined in the nodes of uniform grid ω̄h = {(ai, tj) |ai = ih,i = 1, ...,N,tj = jh, j = 1, ...,M,N = ad/h,M = T/h} . The mesh spacing h, the same for the variables a and t, is chosen so that 0 < h ≤ min(θ,ac), θ/h − [θ/h] = 0, ac/h − [ac/h] = 0, ar/h − [ar/h] = 0, ag/h− [ag/h] = 0, am/h− [am/h] = 0, ad/h− [ad/h] = 0. In this case each point a (k) l ∈ Ω (k) (Eq. (19)) and ã (k) l ∈ Ω̃ (k) (Eq. (54)) coincides always with some point tj. The numerical HPV quantity Vj is defined at the points tj. We introduce the special i- and M- indexes: ic = ac/h, ir = ar/h, ig = ag/h, im = am/h, Mθ = T/θ. Numerical density of susceptible cells Sji is defined from Eq. (23): S0i = ϕ(ai), i = 0, ...,N,j = 0, (67) for k = 1: S j i =   ( F (0) 1 ) i−j (W01) j i−j, 0 , if 1 ≤ j < i ≤ N,( F (1) 1 ) j−i (W1) j j−i, j−i , if 0 ≤ i ≤ j ≤ N, (68) for k > 1: S j i =   ( F (k−1) 1 ) j−(k−2)N−i (W1) j j−i(k−1)N , if 1 ≤ (j − (k − 1)N) < i,( F (k) 1 ) j−(k−1)N−i (W1) j j−i ,j−i , if i ≤ (j − (k − 1)N) ≤ N, (69) where ( F (k) 1 ) l , (W01) j i,p, and (W1) j i,p are given:( F (0) 1 ) l = ϕ(al), l = 0, ...,N, (70)( F (k) 1 ) l = Tr ( g (k) 1 ,e (k) 1 , ir, im, l ) , (71) k ≥ 1, l = 0, ...,N, ( g (k) 1 ) l,i =   βs(ai, tl) ( F (0) 1 ) i−l (W01) l i−l,0 , if k = 1, βs(ai, tl+(k−1)N ) ( F (k−1) 1 ) l−i+N (W1) l+(k−1)N l+(k−1)N−i, (k−1)N , if k > 1, (72) ( e (k) 1 ) l,i = βs(ai, tl+(k−1)N ) ( F (k) 1 ) l−i (W1) l+(k−1)N l+(k−1)N−i,l+(k−1)N−i , (73) k = 1, ...,K, Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 11 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... (W01) j i,0 = exp ( −0.5h ( (f1) 0 i + ( f (0) 1 )j i+j ) − h j−1∑ l=1 (f1) l i+l ) , (74) (W1) j i,p=   1, if j = p, exp ( −0.5h ( (f1) p p−i+(f1) j j−i ) −h j−1∑ l=p+1 (f1) l l−i ) , if j > p, (75) (f1) j i = ds(ai, tj) + α(ai −θ,tj −θ)Ṽ (tj −θ), (76) where we use the trapezoidal rule [22] for approx- imating the integrals in Eqs. (27)-(29) with the order of approximation O(h2): Tr(x,y,ir,im,l) =   h ( 0.5 ( (x)l,ir+(x)l,im ) + im−1∑ i=ir+1 (x)l,i ) , if 0 ≤ l ≤ ir, h ( 0.5 ( (x)l,ir+(x)l,l ) + l−1∑ i=ir+1 (x)l,i ) +h ( 0.5 ( (y)l,l+(y)l,im ) + im−1∑ i=l+1 (y)l,i ) , if ir < l < im, h ( 0.5 ( (y)l,ir +(y)l,im ) + im−1∑ i=ir+1 (y)l,i ) , if im ≤ l ≤ N, (77) The smooth function Ṽ (tj − θ) is obtained from the initial value V0(t) and numerical virus density V j by parabolic interpolating procedure [22]. We imply here that sums like j∑ i=l fi in all equations are omitted if j < l. The numerical func- tion ( F (k) 1 ) l depends from the numerical function ( e (k) 1 ) li (Eqs. (71), (73)) which uses the value of( F (k) 1 ) i−l obtained at the previous steps (index l > i − l in Eq.(72)). The numerical density of infectious cells I1i is given: I0i = 0, i = 0, ...,N, j = 0, (78) for k = 1: I j i =   ( Z (0) 2 )j i−j ,0 if 1 ≤ j < i ≤ N,( F (1) 2 ) j−i (W2) j j−i ,j−i+ ( Z (1) 2 )j j−i,j−i , if 0 ≤ i ≤ j ≤ N, (79) for k > 1: I j i =   ( F (k−1) 2 ) j−(k−2)N−i (W2) j j−i,(k−1)N + ( Z (k−1) 2 )j j−i,(k−1)N , if 1 ≤ (j − (k − 1)N) < i,( F (k) 2 ) j−(k−1)N−i (W2) j j−i,j−i + ( Z (k) 2 )j j−i,j−i , if i ≤ (j − (k − 1)N) ≤ N, (80) where ( F (k) 2 ) i , (W2) j i,p and ( Z (k) 2 )j i,p are given: ( F (0) 2 ) l = 0, l = 0, ...,N, (81)( F (k) 2 ) l = Tr ( g (k) 2 ,e (k) 2 , ir, im, l ) , (82) k ≥ 1, l = 0, ...,N, ( g (k) 2 ) l,i =   βq(ai, tl) ( Z (0) 2 )l i−l,0 , if k = 1, βq(ai, tl+(k−1)N ) (( F (k−1) 2 ) l−i+N (W2) l+(k−1)N l+(k−1)N−i,(k−1)N + ( Z (k−1) 2 )l+(k−1)N l+(k−1)N−i,(k−1)N ) , if k > 1, (83) Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 12 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ...( e (k) 2 ) l,i = βq(ai, tl+(k−1)N ) (( F (k) 2 ) l−i (W2) l+(k−1)N l+(k−1)N−i,l+(k−1)N−i + ( Z (k) 2 )l+(k−1)N l+(k−1)N−i,l+(k−1)N−i ) , (84) k = 1, ...,K, (W02) j i,0 = exp ( −0.5h ( d̃q(ai, t0) + d̃q(ai+j, tj) ) −h j−1∑ l=1 d̃q(ai+l, tl) ) (85) (W2) j i,p=   1, if j = p, exp ( −0.5h ( d̃q(ap−i, tp) +d̃q(aj−i, tj) ) −h j−1∑ l=p+1 d̃q(al−i, tl) ) , if j > p, (86) ( Z (0) 2 )j i,0 = 0.5h ( (f02) j i,0 + (f02) j i,j ) + h j−1∑ l=1 (f02) j i,l, i = j, ...,N, (87) ( Z (k) 2 )j i,p =   0, if j = p, 0.5h ( (f2) j i,p + (f2) j i,j ) +h j−1∑ l=p+1 (f2) j i,l, i = 1, ...,N, if j > p, (88) (f02) j i,p= (W02) j i,pα(ai+p−θ,tp−θ)Ṽ (tp−θ)S p i+p, (89) (f2) j i,p= (W2) j i,pα(ap−i−θ,tp−θ)Ṽ (tp−θ)S p p−i, (90) where in Eqs.(82), (86), (88), we use the trape- zoidal rule [22] for approximating the inte-grals with the order of approximation O(h2). The nu- merical density of precancerous cells Pji is given: P0i = 0, i = 0, ...,N, j = 0, (91) for k = 1: P j i =   ( Z (0) 3 )j i−j,0 , if 1 ≤ j < i ≤ N,( F (1) 3 ) j−i (W3) j j−i,j−i+ ( Z (1) 3 )j j−i,j−i , if 0 ≤ i ≤ j ≤ N, (92) P j i =   ( F (k−1) 3 ) j−(k−2)N−i (W3) j j−i, (k−1)N + ( Z (k−1) 3 )j j−i,(k−1)N , if 1 ≤ (j − (k − 1)N) < i,( F (k) 3 ) j−(k−1)N−i (W3) j j−i,j−i + ( Z (k) 3 )j j−i,j−i , if i ≤ (j − (k − 1)N) ≤ N, (93) where ( F (k) 3 ) i , (W3) j i,p and ( Z (k) 3 )j i,p are given: ( F (0) 3 ) l = 0, l = 0, ...,N. (94)( F (k) 3 ) l = (1 −µ(k)l )Tr ( g (k) 3 ,e (k) 3 , ir, im, l ) , k ≥ 1, l = 1, ...,N, (95) µ (k) l = Tr ( g̃ (k) 3 , ẽ (k) 3 , ir, id, l ) ,k ≥ 1 (96) ( g (k) 3 ) l,i =   βr(ai, tl) ( Z (0) 3 )l i−l,0 , if k = 1, βr(ai, tl+(k−1)N ) (( F (k−1) 3 ) l−i+N (W3) l+(k−1)N l+(k−1)N−i,(k−1)N + ( Z (k−1) 3 )l+(k−1)N l+(k−1)N−i,(k−1)N ) if k > 1, (97) ( e (k) 3 ) l,i = βr(ai, tl+(k−1)N ) (( F (k) 3 ) l−i (W3) l+(k−1)N l+(k−1)N−i,l+(k−1)N−i + ( Z (k) 3 )l+(k−1)N l+(k−1)N−i,l+(k−1)N−i ) , (98) Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 13 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... ( g (k) 3 ) l,i =   ( Z (0) 3 )l i−l,0 , if k = 1,( F (k−1) 3 ) l−i+N (W3) l+(k−1)N l+(k−1)N−i,(k−1)N + ( Z (k−1) 3 )l+(k−1)N l+(k−1)N−i,(k−1)N , if k > 1, (99) ( ẽ (k) 3 ) l,i = ( F (k) 3 ) l−i (W3) l+(k−1)N l+(k−1)N−i,l+(k−1)N−i + ( Z (k) 3 )l+(k−1)N l+(k−1)N−i,l+(k−1)N−i , (100) (W03) j i,0 = exp (−0.5h (dr(ai, t0) + dr(ai+j, tj)) −h j−1∑ l=1 dr(ai+l, tl) ) , (101) (W3) j i,p =   1, if j = p, exp (−0.5h (dr(ap−i, tp) +dr(aj−i, tj)) −h j−1∑ l=p+1 dr(al−i, tl) ) , if j > p, (102) ( Z (0) 3 )j i,0 = 0.5h ( (f03) j i,0 + (f03) j i,j ) + h j−1∑ l=1 (f03) j i,l, i = j, ...,N, (103) ( Z (k) 3 )j i,p =   0, if j = p, 0.5h ( (f3) j i,p + (f3) j i,j ) +h j−1∑ l=p+1 (f3) j i,l, i = 1, ...,N, if j > p, (104) (f03) j i,p = (W03) j i,p δ(ai+p, tp)I p i+p, (105) (f3) j i,p = (W3) j i,p δ(ap−i, tp)I p p−i, (106) where we use in Eqs.(95), (96), (101) - (104) the trapezoidal rule [22] for approximating the inte- grals with accuracy O(h2). The numerical density of cancer cells Cji is given: C0i = 0, i = 0, ...,N,j = 0, (107) for k = 1: C j i =   0, if 1 ≤ j < i ≤ N,( F (1) 4 ) j−i (W4) j j−i,j−i , if 0 ≤ i ≤ j ≤ N, (108) for k > 1: C j i =   ( F (k−1) 4 ) j−(k−2)N−i (W4) j j−i, (k−1)N , if 1 ≤ (j − (k − 1)N) < i,( F (k) 4 ) j−(k−1)N−i (W4) j j−i,j−i , if i ≤ (j − (k − 1)N) ≤ N, (109) where ( F (k) 4 ) i and (W4) j i,p are given:( F (0) 4 ) l = 0, l = 0, ...,N, (110) ( F (k) 4 ) l = Tr ( g (k) 4 ,e (k) 4 , ic, in, l ) + µ (k) l h( P l+(k−1)N ir βr(ar, tl+(k−1)N ) + P l+(k−1)N im ×βr(am, tl+(k−1)N ) +2 im−1∑ p=ir+1 P l+(k−1)Np βr(ap, tl+(k−1)N )   , k ≥ 1, l = 1, ...,N, (111) ( g (k) 4 ) l,i =   0, if k = 1, βc(ai, tl+(k−1)N ) ( F (k−1) 4 ) l−i+N (W4) l+(k−1)N l+(k−1)N−i,(k−1)N , if k > 1, (112) ( e (k) 4 ) l,i = βc(ai, tl+(k−1)N ) ( F (k) 4 ) l−i (W4) l+(k−1)N l+(k−1)N−i,l+(k−1)N−i (113) Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 14 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... (W4) j i,p=   1, if j = p, exp (−0.5h (dc(ap−i, tp) +dc(aj−i, tj)) −h j−1∑ l=p+1 dc(al−i, tl) ) , if j > p, (114) where we use in Eqs. (111), (114) the trape- zoidal rule [22] for approximating the integrals with accuracy O(h2). The exact solution (66) of initial problem (5), (11) for HPV density can be written in a more convenient form for numerical implementation: V (t) = V (t−h) exp  − t∫ t−h dv(V (ξ −σ))dξ   + t∫ t−h exp  − t∫ η dv(V (ξ −σ))dξ  (Λ(η)+ n(η) ad∫ 0 (d∗q(a,η)I(a,η)+d ∗ r(a,η)P(a,η))da ) dη. (115) Using the trapezoidal rule [22] for approximat- ing the integrals in Eq. (115) yields the numerical density of virus population: V j = V j−1 exp ( −0.5h ( dv(Ṽ (tj −σ)) +dv(Ṽ (tj−1 −σ)) )) + 0.5h ( (g5) j + (g5) j−1 exp ( −0.5h ( dv(Ṽ (tj −σ)) +dv(Ṽ (tj−1 −σ)) ))) , (116) (g5) j = Λ(tj) + n(tj) ( 0.5h((e5) j 0 + (e5) j id ) +h id−1∑ i=1 (e5) j i ) , (117) (e5) j i = d ∗ q(ai, tj)I j i + d ∗ r(ai, tj)P j i . (118) The above results lead to the following Theo- rem. Theorem 2. Let the conditions of the The- orem 1 are hold. The numerical method for the recurrent formulae (23)-(30), (33)-(40), (43)- (53), (58)-(66) - exact solution of the Sys- tem (1)-(11), developed for the uniform grid ω̄h = {(ai, tj) |ai = ih,i = 1, ...,N,tj = jh, j = 1, ...,M,N = ad/h,M = T/h} is given by Eqs. (67)-(118). The method is based on the trape- zoid rules and approximates the exact solution of the System (1)-(11) with accuracy ψ = O(h2). Along with the study of the order of approxima- tion of numerical method it is im-portant to obtain the estimations of the velocity of conversation of numerical solution to exact one. This question is studied in the next, second part of our paper. In the next section we consider a series of numerical experiments which illustrate the theoretical results ob-tained in Theorems 1 and 2. IX. NUMERICAL EXPERIMENTS The theoretical results obtained in Theorem 2 are illustrated by the numerical experiments. In the first group of experiments we study the residual as a measure of deviation of obtained numerical solution from some benchmark solution of System (1)-(11) taken from [8], [37], [56]. It is easy to verify by direct substitution that functions S(a,t) = I(a,t) = P(a,t) = C(a,t) = exp(−εa)(1 + exp(−t))−1, (119) V (t) = (2σ + t) 0.5 , (120) are exact solution of the system (1)-(11) with the following particular coefficients and initial values: βs(a,t)=βq(a,t) =ε(exp(−εar)−exp(−εam))−1 = const, (121) µ(t) = ρ (ε (1 + exp(−t)) (exp(−εar) −exp(−εad)) −1 + w ) , (122) βr(t)=ε(1−µ(t))−1(exp(−εar)−exp(−εam))−1, (123) βc(t)= (ε−µ(t)βr(t) (exp(−εar)−exp(−εam))) ×(exp(−εac)−exp(−εag))−1, (124) Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 15 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... ds(t) =ε−α(1+exp(−t))−0.5−(1 + exp(t))−1, (125) dq(t) =ε + α(1 + exp(−t))−0.5 −δ − (1 + exp(t))−1, (126) dr(t) = ε + δ − (1 + exp(t))−1, (127) d∗q(t) = 0.2dq(t), (128) d∗r(t) = 0.2dr(t), (129) dc(t) = ε− (1 + exp(t))−1, (130) dv(V (t−σ)) = 0.05(1 + V (t−σ))0.5, (131) Λ(t) = 0.5(2σ + t) −0.5 + 0.05 ( 1 + (σ + t) 0.5 )0.5 (2σ + t) 0.5 −nε−1 ( d∗q(t) + d ∗ r(t) ) ×(1−exp(−εad)) (1+exp(−t)) −1 , (132) S(a, 0) = I(a, 0) = P(a, 0) = C(a, 0) = 0.5 exp(−εa), (133) V0(t) = (2σ + t) 0.5 , t ∈ [−σ, 0]. (134) Exact solutions (119) satisfy the continuity condi- tion (30) and conditions: I(0, 0) = am∫ ar βs(a, 0)I(a, 0)da, (135) P(0, 0) = (1 −µ(Nr(0)) am∫ ar βr(a, 0)P(a, 0)da, (136) C(0, 0) = ag∫ ac βc(a, 0)C(a, 0)da + µ(Nr(0)) am∫ ar βr(a, 0)P(a, 0)da. (137) The set of constants used in experiments is given in Table I. Our goal here is to pro-vide a detailed description of the numerical experiments and obtained numerical results which will help the potential users to evaluate the numerical method of characteristics and increase the chances of its im- plementation in simulation of epidemic dynamics in various applications. TABLE I SET OF CONSTANTS Constant Value Constant Value ar 0.3 σ 0.2 am 0.9 n 1.5 ac 0.1 α 0.01 ag 0.4 δ 0.01 ad 1 ρ 0.1 T 10 w 0.5 θ 0.1 ε 4 The time of simulation in all experiments equals to 10 cell’s lifespan T = 10ad. Deviation of the numerical solution Y ji = (S j i ,I j i ,P j i ,C j i ) from the benchmark one Y (ai, tj) = (S(ai, tj),I(ai, tj),P(ai, tj),C(ai, tj)) and V j from V (tj) with x = h/ad → 0 (dimensionless parameter) is measured by the two dimensionless norms of functional spaces H and C [8], [51]: δ (1) H (Y ) = M∑ j=0 N∑ i=0 ∣∣∣Y ji −Y (ai, tj)∣∣∣  M∑ j=0 N∑ i=0 ∣∣∣Y ji ∣∣∣  −1, (138) δ (1) C (Y ) = max 0≤i≤N 0≤j≤M ∣∣∣Y ji −Y (ai, tj)∣∣∣   max 0≤i≤N 0≤j≤M ∣∣∣Y ji ∣∣∣  −1, (139) δ (2) H (V ) = M∑ j=0 ∣∣V j −V (tj)∣∣   M∑ j=0 ∣∣V j∣∣  −1, (140) δ (2) C (V ) = max 0≤j≤M ∣∣V j −V (tj)∣∣( max 0≤j≤M ∣∣V j∣∣)−1, (141) All program projects were created in Microsoft Visual C 2019, platform Microsoft .Net Frame- work and were launched on a PC (CPU i5 of 8- th generation 1.6 GHz, RAM 16 Gb). The results Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 16 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... of numerical experiments are presented in Table 2 and are shown in Fig. 3 ( � - the values of δ(1)H (S), δ (1) C (S) from Table II depending from x = h/ad, - - - - the graphs of corresponding regressions y(x) (with denoted equations)), Fig. 4 (δ(1)H (I), δ (1) C (I) and y(x)), Fig. 5 (δ(1)H (P), δ (1) C (P) and y(x)), Fig. 6 (δ(1)H (C), δ (1) C (C) and y(x)) and Fig. 7 (δ (2) H (V ), δ (2) C (V ) and y(x)). The session time of simulation Ts (sec) (diagnostic session time of Microsoft Visual C 2019) depending from parameter x in each numerical experiment is shown in Table III. TABLE II THE VALUES OF δ(1)H (Y ), δ (1) C (Y ), Y j i = (S j i ,I j i ,P j i ,C j i ), AND δ (2) H (V ), δ (2) C (V ) FOR THE DIFFERENT VALUES OF x. x 0.005 0.01 0.025 0.05 δ (1) H (S) 0.00038 0.0015 0.0094 0.0371 δ (1) C (S) 0.00074 0.0023 0.0184 0.0716 δ (1) H (I) 0.00035 0.0014 0.0088 0.0347 δ (1) C (I) 0.00066 0.0026 0.0164 0.0639 δ (1) H (P) 0.00035 0.0014 0.0086 0.0342 δ (1) C (P) 0.00065 0.0026 0.0161 0.0629 δ (1) H (C) 0.00077 0.0031 0.0193 0.0757 δ (1) C (C) 0.0015 0.0054 0.0368 0.1397 δ (1) H (V ) 0.00021 0.0008 0.0051 0.0210 δ (1) C (V ) 0.00036 0.0014 0.0091 0.0382 TABLE III THE SESSION TIME OF SIMULATION Ts IN SECONDS (DIAGNOSTIC SESSION TIME OF MI-CROSOFT VISUAL C 2019), DEPENDING FROM THE x. x 0.005 0.01 0.025 0.05 Ts (sec) 651 52 4 2 The equations of regressions in Figs. 3 - 7 are built automatically in MS Excel at the points of δ(1)H (Y ), δ (1) C (Y ), δ (2) H (V ), δ (2) C (V ) depending from the values of x. They exhibit the quadratic (parabolic) low of convergence δ(1)H (Y ) → 0, δ (1) C (Y ) → 0, (Y j i = (S j i ,I j i ,P j i ,C j i )) and δ (2) H (V ) → 0, δ (2) C (V ) → 0 with x → 0 in all ex- periments that illustrates and confirms the second Fig. 3. � - values of δ(1)H (S), δ (1) C (S), - - - graphs of regressions y(x), x = h/ad. Fig. 4. � - values of δ(1)H (I), δ (1) C (I), - - - graphs of regressions y(x), x = h/ad. Fig. 5. � - values of δ(1)H (P), δ (1) C (P), - - - graphs of regressions y(x), x = h/ad. order of accuracy of numerical solution. The most numerical error among all 5 subclasses is obtained for the cancer cell population (Fig. 6, Table II - δ(1)H (C), δ (1) C (C)). The large deviation (> 3%) of numerical cancer cell population density from the benchmark solution for x ≥ 0.025 can be Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 17 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... Fig. 6. � - values of δ(1)H (C), δ (1) C (C), - - - graphs of regressions y(x), x = h/ad. Fig. 7. � - values of δ(1)H (V ), δ (1) C (V ), - - - graphs of regressions y(x), x = h/ad. explained by the complexity of nonlinear boundary condition (9) which includes the integral from the numerical precancerous cell population density com-puted with some numerical error. Neverthe- less, we observe the significand decreasing of δ (1) H (C), δ (1) C (C) for x ≤ 0, 01 with follow-ing convergence of them to 0. Besides the accuracy of numerical method, we consider also the valuable in practice parameter - session time of simulation Ts given in Table III. Since the refine of difference grid leads necessarily to increasing of Ts, numer- ical experiments with benchmark solution provide the optimal range of mesh spacing h for which the appropriate accuracy of numerical solution can be reachable for an acceptable time of simulation. According to the data of Table III the critical point of mech spacing after which the time of simulation begins to grow rapidly is h = 0.01 (or x = 0.01). On the other hand, the results of Fig. 8. Graphs of NS(t) for non-oscillating dynamics for 3 different initial values ϕ(a). Table II show that for x = 0.01 the deviation of numerical solution from the benchmark one (relative numerical error) is less or equal than 0, 5% for all subclasses of population that corre- sponds to the acceptable accuracy of simulation in most biological systems. From Tables II and III, it follows that the optimal value of the dimensionless parameter is x = h/ad = 0.01 that corresponds to the mesh spacing h = 0.01. This result is in good agreement with the recommended value of time spacing obtained in the numerical experiments for nonlinear age-structured model of population dynamics in [8]. In the second group of numerical experiments we study the asymptotically stable dynamical regimes of population. Simulations reveal two types of asymptotically stable dynamics: non- oscillating and oscillating convergence of all sub- population quantities to the nontrivial equilibrium states. The dynamical regime of the first type - non-oscillating dynamics is shown in Figs. 8 - 12. In particular, in Fig. 8 it is shown the graphs of susceptible cell quantity dynamics for three different initial values ϕ(a). The corresponding graphs of infected cell quantity, precancerous cell quantity, cancer cell quantity and HPV quantity dynamics are shown in Figs. 9, 10, 11, 12, respec- tively. It should be noted that the non-oscillating dy- namics of dysplasia (precancerous cells) and can- cer growth shown in Figs. 10 and 11 by dotted lines is a most realistic type of tumor tissue Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 18 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... Fig. 9. Graphs of NI (t) for non-oscillating dynamics for 3 different initial values ϕ(a). Fig. 10. Graphs of NP (t) for non-oscillating dynamics for 3 different initial values ϕ(a). Fig. 11. Graphs of NC (t) for non-oscillating dynamics for 3 different initial values ϕ(a). growth. The second type of dynamical regimes - os- cillating dynamics is shown in Figs. 13 - 17. In Fig. 13, it is shown the susceptible cell quantity dynamics for two different initial values ϕ(a). The corresponding graphs of infected cell quantity, Fig. 12. Graphs of V (t) for non-oscillating dynamics for 3 different initial values ϕ(a). Fig. 13. of NS(t) for oscillating dynamics for 2 different initial values ϕ(a). precancerous cell quantity, cancer cell quantity and HPV quantity dynamics are shown in Figs. 14, 15, 16, 17, respectively. Asymptotically stable regimes of SIPCV model shown in Figs. 8 - 12 and 13 - 17 demonstrate the localization of dysplasia (precancerous cells) and cancer in biological tissue without metastases. Thus, the results of simulations exhibit that the numerical method designed for the age-structured SIPCV epidemic model can be applied for numer- ical study and modelling of cell-HPV population dynamics with high accuracy. X. CONCLUSION In this study we consider a new epidemic model of age-structured sub-populations of sus- ceptible, infectious, precancerous and cancer cells and unstructured population of hu-man papilloma virus (HPV) (SIPCV epidemic model). The model Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 19 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... Fig. 14. of NI (t) for oscillating dynamics for 2 different initial values ϕ(a). Fig. 15. of NP (t) for oscillating dynamics for 2 different initial values ϕ(a). Fig. 16. of NC (t) for oscillating dynamics for 2 different initial values ϕ(a). is based on the competi-tive system of initial- boundary value problems for delayed semi-linear transport equations with integral boundary con- ditions and initial problem for delayed nonlinear ODE. For this system we obtained the exact so- lution in the form of recurrent formulae in which Fig. 17. of V (t) for oscillating dynamics for 2 different initial values ϕ(a). the den-sities of all subpopulations depend from the integrals from solution taken at previous in- stance of time. The main ideas and method of derivation of such exact solutions are taken from works [3], [8]. The new result obtained in this work is the numerical implementation of recurrent formulae for exact solution and development of the numerical method of the second order of approx- imation for simulation of susceptible, infectious, precancerous, can-cer cells and human papilloma virus population dynamics. Since evaluation of the accuracy of numerical solution and session time of simulation are essential to successful use of nu- merical method in applications, we estimated the difference between computed solution and bench- mark solution of model and session time on the refined difference grid. Numeri-cal experiments showed that the relative numerical error of solution may be reduced up to 0.1% for the very small value of mesh spacing parameter h = 0.005 (that is 0.5% of ad) but for large value of session time of simulation. We recommend to use in applications the optimal value of mesh spacing h = 0.01 (1% of ad) that pro-vides the small value of relative numerical error (less than 0,5%) for acceptable session time. This recommendation is in a good agreement with the results of numerical experi- ments obtained in [8] for nonlinear age-structured model of population dynamics. The numerical experiments with model pa- rameters revealed two types of asymptotically stable dynamical regimes of SPICV population Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 20 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... non-oscillating and oscillating convergence of so- lution to the positive steady states. The non- oscillating type of SPICV population dynamics corresponds to the observed in practice dynamics of tumor growth-localization of dysplasia (pre- cancerous cells) and cancer in biological tissue without metastases. Overall, development of age- structured SIPCV epidemic model, derivation of its ex-act solution and design of corresponding nu- merical methods provide the theoretical instrument for study the dynamics of susceptible, infected, precancerous, cancerous cells and viruses popu- lations that help us better understand the features of human papilloma virus infectious disease. ACKNOWLEDGMENT We would like to acknowledge Department of Mathematics, Universitas Gadjah Mada for the organizational support of these studies. REFERENCES [1] L. M. Abia, O. Angulo, J. C. Lopez-Marcos, Age- structured population models and their numerical solu- tion, Ecol. Modell. 188(1) (2005), 112-136. [2] L. M. Abia, O. Angulo, J. C. Lopez-Marcos, M. A. Lopez-Marcos, Numerical schemes for a size-structured cell population model with equal fission, Math. Comput. Modell. 50 (2009) 653-664. [3] V. V. Akimenko, An age-structured SIR epidemic model with the fixed incubation pe-riod of infection, Comput. Math. Appl. 73 (2017) 1485?150. [4] V. V. Akimenko, Asymptotically stable states of non- linear age-structured monocyclic cell population model I. Travelling wave solution. Math. Comput. Simul. 133 (2017) 2-23. [5] V. V. Akimenko, Asymptotically stable states of non- linear age-structured monocyclic cell population model II. Numerical simulation. Math. Comput. Simul. 133 (2017) 24-38. [6] V. V. Akimenko, Continuous monocyclic and polycyclic age structured models of population dynamics. Commun. Biomath. Sci. 2 (2) (2019) 85-95. [7] V. V. Akimenko, Modelling of two-dimensional transport processes by using of non-linear monotonous second order schemes, Cybern. Syst. Anal. 39 (6)(2003) 839-853. [8] V. V. Akimenko, Nonlinear age-structured models of polycyclic population dynamics with death rates as a power functions with exponent n. Math.Comput.Simul. 133 (2017) 175-205. [9] V. V. Akimenko, Stability analysis of delayed age- structured resource-consumer mod-el of population dy- namics with saturated intake rate. Front. Ecol. Evol. 9:531833 (2021). doi: 10.3389/fevo.2021.531833. [10] V. V. Akimenko, The maximum principle and nonlin- ear monotone schemes for parabolic equations, Comput. Math. Math. Phys. 39 (4) (1999), 590-600. [11] V. V. Akimenko, R. Anguelov, Steady states and out- breaks of two-phase nonlinear age-structured model of population dynamics with discrete time delay, J. Biol. Dynam. 11 (1) (2017) 75-101. [12] V. V. Akimenko, V. K?ivan, Asymptotic stability of delayed predator age-structured population models with an Allee effect. Math. Biosci. 306 (2018) 170-179. [13] M. Al-arydah, R. Smith, An age-structured model of human papillomavirus vaccination. Math. Comput. Simul. 82 (2011) 629?652. [14] O. Angulo, J. C. L?opez-Marcos, M. A. L?opez-Marcos, Analysis of an efficient inte-grator for a size-structured population model with a dynamical resource, Comput. Math. Applic. 68 (9) (2014) 941?961. [15] O. Angulo, J. C. Lopez-Marcos, M. A. Lopez-Marcos, F. A. Milner, A numerical method for nonlinear age- structured population models with finite maximum age, J. Math. Anal. Appl. 361 (2010) 150?160. [16] O. Angulo, J. C. L?opez-Marcos, Numerical integra- tion of fully nonlinear size-structured population models, Appl. Numer. Math. 50 (2004) 291-327. [17] 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 (10) (2018) 1535-1546. [18] F. Billy, J. Clairambault, F. Delaunay, C. Feillet, N. Robert, Age-structured cell population model to study the influence of growth factors on cell cycle dynamics. Math. Biosci. Eng. 10 (1) (2013) 1-17. [19] V. S. Borisov, S. Sorek, On monotonicity of difference schemes for computational physics, SIAM J. Sci. Comput. 25 (5) (2004), 1557 ? 1584. [20] F. M. Burnet, Intrinsic Mutagenesis: A Genetic Ap- proach to Ageing, John Wiley & Sons Inc., New York, 1974. [21] F. Brauer, C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology, Springer, New York, 2012. [22] R. L. Burden, J. D. Faires, A. M. Burden, Numerical Analysis, Cengage Learning, Boston, 2015. [23] F. Chang, D. Drubin, Cell division: Why daughters cannot be like their mothers. Curr. Biol. 6 (1996) 651- 654. [24] A. Coronel, L. Friz, I. Hess, M. Zegarra, On the existence and uniqueness of an inverse problem in epidemiology. Applicable Analysis 100(3) (2019) DOI: 10.1080/00036811.2019.1608964. [25] A. M. de Ross, Numerical methods for structured population models: The escalator boxcar train, Numer. Methods for Partial Differential Equations 4 (1988) 173- 195. [26] N. M. Dixit and A. S. Perelson. Complex patterns of viral load decay under antiretroviral therapy: influence of Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 21 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... pharmacokinetics and intracellular delay. J. Theor. Biol., 226 (2004) 95?109. [27] P. J. Eifel, A. H. Klopp, J. S. Berek, P. A. Kon- stantinopoulos, Cancer of the cervix, vagina and vulva, in: V. T. DeVita Jr, T. S. Lawrence, S. A. Rosenberg (Eds.), Cancer: Principles and Practice of Oncology (11th Edition), Wolters Kluwer, Philadelphia, 2019, pp. 2083 - 2151. [28] L. E. Elsgolts, S. B. Norkin, Introduction to the Theory and Application of Differential Equations with Deviating Arguments, Mathematics in Science and Engineering, V. 105. Academic Press, New York - London, 1973. [29] J. Z. Farkas, G.F. Webb, Mathematical analysis of a clonal evolution model of tumour cell proliferation. J. Evol. Equ. 17 (2017), 275?308. [30] L. M. Franks, M. A. Knowles, What is cancer? in: M. A. Knowles, P. J. Selby (Eds.), Introduction to the Cellular and Molecular Biology of Cancer, Oxford University Press, Oxford, 2005, pp. 1 ? 24. [31] K. Fernald, M. Kurokawa, Evading apoptosis in cancer. Trends Cell Biol. 23(12) (2013) 620?633. [32] A. Gandolfi, M. Iannelli, G. Marinoschi, An age- structured model of epidermis growth, J. Math. Biol. 62 (1) (2011) 111-141. [33] S.K. Godunov, A method for numerical calculation of discontinuous solutions to flu-id-dynamic equations, Mathem. Sb. 47 (3) (1959) 271-306. [34] P. Krzyzanowski, D. Wrzosek, D. Wit, Discontinuous Galerkian method for piece-wise regular solution to the nonlinear age-structured population model, Math. Biosci. 203(2) (2006) 277-300. [35] W. Krzyzanski, Pharmacodynamic models of age- structured cell populations. J. Pharmacokinet. Pharma- codyn. 42 (2015) 573?589. [36] C.A. Kunos, F.W. Abdul-Karim, D.S. Dizon, R. De- bernardo, Cervix uteri, in: D.S. Chi, D.S. Dizon, A. Berchuck, C. Yashar (Eds.), Principles and Practice of Gynecologic Oncology (7th Edition), Wolters Kluwer, Philadelphia, 2017, pp. 946 ? 983. [37] X. Y. Li, Variational iteration method for nonlinear age- structured population models. Comput. Math. Applic. 58 (11-12) (2009) 2177-2181. [38] J. Li, F. Brauer, Continuous-Time Age-Structured Mod- els in Population Dynamics and Epidemiology, in: F. Brauer, P. van den Driessche, J. Wu (Eds.), Mathemat- ical Epidemiology, Lecture Notes in Mathematics 1945, Springer, Berlin, 2008, pp. 205-227. [39] X.-Z. Li J. Yang, M. Martcheva, Age Structured Epi- demic Modelling, Springer, New York, 2020. [40] 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. [41] Z. Liu, C. Guo, H. Li, L. Zhao, Analysis of a nonlinear age-structured tumor cell population model. Nonlinear Dynam. 98 (2019) 283-300. [42] Z. Liu, C. Guo, J. Yang, H. Li, Steady states analysis of a nonlinear age-structured tumor cell population model with quiescence and bidirectional transition. Acta Appl. Math. 169 (2020) 455 ? 474. [43] 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 (8) (2013) 2151-2174. [44] M. Martcheva, An Introduction to Mathematical Epi- demiology, Springer, New York, 2015. [45] V. R. Martin, S. V. Temple, Cervical cancer, in: C. H. Yarbro, D. Wujcik, B. H. Gobel (Eds.) Cancer Nursing: Principles and Practice (7th ed.), Jones and Bartlett Publishers, Sudbury Massachusetts, 2011, pp. 1188 ? 1204. [46] V. Martsenyuk, M. Karpinski, S. Rajba, J. Nikodem, K. Warwas, L. Wieclaw, T. Gancarczyk, Global asymptotic stability and nonlinear analysis of the model of the square immunopixels array based on delay lattice differential equations. Symmetry (2020), 12, 40. [47] J. Muller, Ch. Kuttler, Methods and Models in Math- ematical Biology, Deterministic and Stochastic ap- proaches, Lecture Notes on Mathematical Modelling in Life Sciences, Springer, Berlin, 2015. [48] T. S. Noor-Asih, S. Lenhart, S. Wise, L. Aryati, F. Adi- Kusumo, M.S. Hardianti, J. Forde, The dynamics of HPV infection and cervical cancer cells. Bull. Math. Biol. 78 (2016) 4?20. [49] B. Perthame, Transport Equations in Biology. Frontiers in mathematics, Birkhauser Verlag, Basel, 2007. [50] I. Roeder, M. Herberg, M. Horn, An ?age?-structured model of hematopoietic stem cell organization with ap- plication to chronic myeloid leukemia. Bull. Math. Biol. 71 (2009) 602?626.4. [51] A.A. Samarskii, The Theory of Difference Schemes, v.240. CRC Press, 2001. [52] D. Sulsky, Numerical solution of structured population models. I Age structure, J. Math. Biol. 31 (8) (1993) 817- 839. [53] G.F. Webb, Population Models Structured by Age, Size and Spatial Position, in: P.Magal, S.Ruan (Eds.), Struc- tured Population Models in Biology and Epidemiology, Lecture Notes in Mathematics 1936, Springer, Berlin, pp. 1-49, 2008. [54] C.B.J. Woodman, S.I. Collins, L.S. Young, The natural history of cervical HPV in-fection: unresolved issues. Nat. Rev. Cancer 7 (2007) 11?22. [55] Z. Yin, Y. Yu, Zh. Lu, Stability analysis of an age- structured SEIRS model with time delay, Mathematics 8 (2020) 455. [56] S.A. Yousefi, M. Behroozifar, M. Dehghan, Numerical solution of the nonlinear age-structured population mod- els by using the operational matrices of Bernstein polyno- mials, Appl. Math. Modell. 36 (2012) 945?963. [57] Z. Zhang, S. Kumari, R.K. Upadhyay, A delayed e- epidemic SLBS model for computer virus, Adv. Dif- fer. Equ. 414 (2019). https://doi.org/10.1186/s13662-019- 2341-8 Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 22 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ... [58] R. Zhang, D. Li, S. Liu, Global analysis of an age- structured SEIR model with immigration of population and nonlinear incidence rate. J. Appl. Anal. Comput. 9 (4) (2019) 1470-1492. [59] J. Zhou, L. Song, J. Wei, Mixed types of waves in a dis- crete diffusive epidemic mod-el with nonlinear incidence and time delay, J. Differential Equations 268 (8) (2020) 4491-4524. Biomath 10 (2021), 2110027, http://dx.doi.org/10.11145/j.biomath.2021.10.027 Page 23 of 23 http://dx.doi.org/10.11145/j.biomath.2021.10.027 Introduction Model Susceptible cell population dynamics Infectious cell population dynamics Precancerous cell population dynamics Cancer cell population dynamics HPV population dynamics Numerical implementations Numerical experiments Conclusion References