Age-structured delayed SIPCV epidemic model of HPV and cervical cancer cells dynamics I. Numerical method

— 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 coefﬁcients 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 ﬁrst 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.


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 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 motivates 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 proliferation 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 tolerant with respect to its own cells [27], [30], [36], [45]; (ii) death rate of HPV is densitydependent 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 densitydependent 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 simulation 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 transport 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 nonlinear 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-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 implementations 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 instance 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 exact 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 relative 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 oscillating dynamics of population when the quantity 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 d s (a, t), infectious d q (a, t), precancerous d r (a, t) and cancer cells d c (a, t) with maximum lifespan a d , an age reproductive window of non-cancer cells [a r , a m ] and cancer cells [a c , a g ], a c < a r , a g < a m . 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 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 d v (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) < d q (a, t), 0 < d * r (a, t) < d r (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 abnormal rapid proliferation (a c < a r , a g < a m ) when cells do not have time for maturation ("nondifferentiable" cells). Cancer cells differ from the normal ones in their lack of response to normal fertility control mechanism [27], [30], [36], [45]. That is why precancerous cells turn to the cancer in our model only through the proliferation when a precancerous cell may divide into two new cancer ones. Fertility rate of precancerous cells turned to the cancer is proportional to the density-dependent saturated rate µ(N r (t)) = ρNr(t) 1+wNr(t) [48], where the quantity of precancerous cells of re-productive age and older is N r (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 ≤ µ(N r ) < ρw −1 for N r ≥ 0. These assump-tions lead to the following age-structured epidemic model in domain 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: I(0, t) = am ar β q (a, t)I(a, t)da 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) where the birth (fertility) rates of the susceptible, infectious, precancerous and cancer cells are β s (a, t), β q (a, t), β p (a, t) and β c (a, t), respectively; ϕ(a) is an initial density of susceptible cells, V 0 (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: We impose the following restrictions on the density-dependent HPV death rate and cells' birth and death rates [27], [30], [31], [36], [45]: The positiveness of derivative of densitydependent death rate of HPV in (15) means that increasing of HPV quantity changes the characteristics of intracellular space that result in the organism immune response through the activation of cell immunity (T-killers) and humoral immunity (B-lymphocytes) that leads to the elimination of viruses (i.e. monotone increasing of their death rate). We assume that all coefficients and initial values of system (1)-(11) are twice continuously differentiable functions and have the private derivatives of the second order by all their arguments.

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 , and we define the following sets ( Fig. 1): . We define also the auxiliary set of age intervals: 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:  Renewal boundary condition (6) completes the problem (19), (20). In original variables, solution of problem (21), (22), has a form (k = 1, ..., K): 1 (a−t) are defined from the initial (10) and boundary (6) conditions: 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 − t k−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: 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 ∈ [a r + u, a m + u] 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: is taken from the previous section. Solution of problem (31), (32) has a form (k = 1, ..., K): Functions F (k) 2 (v) are defined from the initial (10), and boundary (7) conditions: where the second term in the right side of Equation (39) should be omitted for l = 2, function F (k) 2 (v) in Eq. (40) has been already defined on the previous steps be-cause its argument v ∈ [a r + u, a m + u] and v > u (see example 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).
where I(v+t, t) is taken from the previous section. Function µ(v +t−σ, N P (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): Functions F (k) 3 (v) and N P (ξ) are defined from the initial (10) and boundary (8) conditions: u ∈ −a 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 defined on the previous steps because its argument v ∈ [a r + u, a m + u] and v > u (see example in Fig.2). Because the initial value of P (a, t) is trivial, solution (43) satisfies the compatibility (continuity) condition of solution, P (a, t) ∈ C(Q).

VI. CANCER CELL POPULATION DYNAMICS
Since cancer cells have a specific age reproductive window [a c , a k ] differed from the one of noncancer cells, we introduce a new auxiliary set of age intervals for the problem (4), (9), (10): 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: where function µ(N P (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): Functions F (k) 4 (v) are defined from the initial (10) and boundary (9) conditions: where the auxiliary functions Φ (k) 4l (v) are defined as: where the second term in the right side of Equation (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 ∈ [a c + u, a g + u] and v > u (see example in Fig.2). Because the initial value of C(a, t) is trivial, solution (57) satis-fies the compatibility (continuity) condition of solution, C(a, t) ∈ C(Q).

VII. HPV POPULATION DYNAMICS
Nonlinear death rate d v (V ) of Cauchy problem for nonlinear delayed ODE (5), (11), satisfies 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]: 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.

VIII. NUMERICAL IMPLEMENTATIONS
The numerical age-specific densities S j i , 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 The mesh spacing h, the same for the variables a and t, is chosen so that 0 < h ≤ min(θ, a c ), (54)) coincides always with some point t j . The numerical HPV quantity V j is defined at the points t j . We introduce the special iand Mindexes: i c = a c /h, i r = a r /h, i g = a g /h, i m = a m /h, M θ = T /θ. Numerical density of susceptible cells S j i is defined from Eq. (23): for k = 1: for k > 1: where F k ≥ 1, l = 0, ..., N, Page 11 of 23 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ...
where we use the trapezoidal rule [22] for approximating the integrals in Eqs. (27)-(29) with the order of approximation O(h 2 ): The smooth functionṼ (t j − θ) is obtained from the initial value V 0 (t) and numerical virus density V j by parabolic interpolating procedure [22]. We imply here that sums like  72)). The numerical density of infectious cells I 1 i is given: for k = 1: for k > 1: where F k ≥ 1, l = 0, ..., N, where in Eqs.(82), (86), (88), we use the trapezoidal rule [22] for approximating the inte-grals with the order of approximation O(h 2 ). The numerical density of precancerous cells P j i is given: for k = 1: where F are given: Page 13 of 23 V V Akimenko, F Adi-Kusumo, Age-structured delayed SIPCV epidemic model of HPV and cervical ...
where we use in Eqs.(95), (96), (101) -(104) the trapezoidal rule [22] for approximating the integrals with accuracy O(h 2 ). The numerical density of cancer cells C j i is given: for k = 1: for k > 1: where F where we use in Eqs. (111), (114) the trapezoidal rule [22] for approximating the integrals with accuracy O(h 2 ). The exact solution (66) of initial problem (5), (11) for HPV density can be written in a more convenient form for numerical implementation: Using the trapezoidal rule [22] for approximating the integrals in Eq. (115) yields the numerical density of virus population: The above results lead to the following Theorem. Along with the study of the order of approximation 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) are exact solution of the system (1)-(11) with the following particular coefficients and initial values: β c (t)= (ε−µ(t)β r (t) (exp(−εa r )−exp(−εa m ))) 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 implementation in simulation of epidemic dynamics in various applications.
is measured by the two dimensionless norms of functional spaces H and C [8], [51]:  Table 2 and are shown in Fig. 3 ( -the values of δ  Table III. C (V ) → 0 with x → 0 in all experiments that illustrates and confirms the second 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 δ 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. Nevertheless, we observe the significand decreasing of δ (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 T s given in Table III. Since the refine of difference grid leads necessarily to increasing of T s , numerical 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 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 corresponds 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/a d = 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: nonoscillating and oscillating convergence of all subpopulation quantities to the nontrivial equilibrium states. The dynamical regime of the first typenon-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, respectively.
It should be noted that the non-oscillating dynamics of dysplasia (precancerous cells) and cancer growth shown in Figs. 10 and 11 by dotted lines is a most realistic type of tumor tissue Fig. 9.
Graphs of NI (t) for non-oscillating dynamics for 3 different initial values ϕ(a).

growth.
The second type of dynamical regimes -oscillating 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). 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 numerical 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 susceptible, infectious, precancerous and cancer cells and unstructured population of hu-man papilloma virus (HPV) (SIPCV epidemic model). The model of NI (t) for oscillating dynamics for 2 different initial values ϕ(a).  is based on the competi-tive system of initialboundary value problems for delayed semi-linear transport equations with integral boundary conditions and initial problem for delayed nonlinear ODE. For this system we obtained the exact solution in the form of recurrent formulae in which the den-sities of all subpopulations depend from the integrals from solution taken at previous instance 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 approximation 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 numerical method in applications, we estimated the difference between computed solution and benchmark 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 a d ) 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 a d ) 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 experiments obtained in [8] for nonlinear age-structured model of population dynamics.
The numerical experiments with model parameters revealed two types of asymptotically stable dynamical regimes of SPICV population non-oscillating and oscillating convergence of solution to the positive steady states. The nonoscillating type of SPICV population dynamics corresponds to the observed in practice dynamics of tumor growth-localization of dysplasia (precancerous cells) and cancer in biological tissue without metastases. Overall, development of agestructured SIPCV epidemic model, derivation of its ex-act solution and design of corresponding numerical methods provide the theoretical instrument for study the dynamics of susceptible, infected, precancerous, cancerous cells and viruses populations that help us better understand the features of human papilloma virus infectious disease.