Age-structured Delayed SIPCV Epidemic Model of HPV and Cervical Cancer Cells Dynamics II. Convergence of Numerical Solution

—The numerical method for simulation of age-structured SIPCV epidemic model with age-structured sub-classes of susceptible, infectious, precancerous and cancer cells and unstructured population of human papilloma virus (HPV) dynamics with incubation period is developed. Convergence of the numerical approximations is studied both theoretically and numerically. We prove the stability and second rate of convergence of the approximate solutions to the exact solution of the SIPCV epidemic nonlinear system. The numerical experiments based on the grid reﬁned method conﬁrm and illustrate the second order of accuracy of the obtained numerical method and show the various dynamical regimes of population dynamics. Simulations for model parameters of the system reveal two unstable dynamical regimes of SIPCV population which correspond to the cancer tumor growth and formation of metastases in organism.


I. INTRODUCTION
Continuous age-structured models have proved to be useful in modelling of a broad set of problems arising in life science because they describe the life history of biological organisms in population through the all stages of their lifecycle: birth, maturation, reproduction and death [3], [4], [6][7][8], [10], [11], [16], [18], [20][21][22][23][24], [26], [30][31][32][33]. Development of the numerical methods for such models plays an important role in modelling of population and evolutionary dynamics and were studied in many works in the different fields of life science [1], [2], [5], [9], [12][13][14]. The qualitative analysis of numerical methods implies the study of consistence (approximation), well-posedness (existence, uniqueness and stability of approximate solution or correctness of difference scheme) and convergence of numerical approximation. Proving of the convergence of numerical approximation to exact solution in direct manner is not easy problem in general case. The solution of this problem is facilitated by a classical Lax-Richtmyer Equivalence Theorem [29] (see also the corresponding theorems in [12], [27], [28]) which says that from the consistency and well-posedness of difference equations follows the convergence of approximate solution to the exact one. This fundamental property of difference schemes motivated extensive development of the stability theory and using it in proving of convergence of numerical solution of linear difference schemes [27][28][29]. But, in the case of nonlinear difference equations the definition of stability is actually equivalent to the definition of convergence of numerical approximation (p.132 in [28]) that makes the qualitative analysis of nonlinear difference schemes more complicated.
The difference scheme is said to be correct or well-posed if it has the unique and stable approximate solution [28]. Since in the first part of our work [9] we focus on the study of existence and uniqueness of exact solution of nonlinear SIPCV model (Theorem 1) and develop the consistent numerical method of the second order of approximation (Theorem 2), in this paper we complete our study and prove the convergence of an approximate solution to the exact solution of the nonlinear system with the second order of accuracy. The convergence of numerical solution follows directly from the stability of numerical method, and the rate of convergence is defined by the evaluation of residual of deviation of the approximate solution from the exact one [28].
The rate of convergence of numerical solution in our paper is studied both theoretically and numerically. In the first part of our work [9] we consider the convergence of approximate solution to the benchmark solution of system in the series of numerical experiments. In this article we study the convergence of approximate solution using the grid re-fined method [19], [28] that correspond to the definition of stability of nonlinear difference schemes [28]. All numerical experiments confirm and illustrate the second rate of convergence of approximate solution of nonlinear SIPCV model. The evaluations of numerical error obtained in this paper are in good agreement with the results of numerical experiments presented in our earlier works [7], [9]. Numerical experiments with model parameters of the system reveal two types of unstable dynamical regimes interesting from biological point of view. The first regime corresponds to the unlimited exponential growth of cancer cells population (when a basic reproduction number greater than one) and asymptotically stable oscillating dynamics of the other cell subpopulations and HPV population. The second regime corresponds to the unlimited exponential growth of dysplasia (precancer cells) which induces the unstable dynamics of cancer cells population and, from the other hand, the asymptotically stable oscillating dynamics of the other cell subpopulations and HPV population. In both cases the unstable cancer cells population dynamics mean the cancer tumor growth and formation of metastases in organism.
Thus, the numerical method of the second order of accuracy designed in our work provides the reliable theoretical instrument for simulation and theoretical study of nonlinear age-structured SIPCV epidemic models.

II. MODEL
We recall that the dynamics of age-specific densities of susceptible, infectious, pre-cancerous and cancer cells S(a, t), I(a, t), P (a, t), C(a, t), respectively, and time-dependent density of HPV V (t) is governed by the nonlinear age-structured model [9]: Equations (1)-(5) are completed by the boundary conditions and initial values: S(a, 0) = ϕ(a), I(a, 0) = 0, P (a, 0) = 0, 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 − θ)∈ C(Q) 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: where α 0 (a, t) is an auxiliary rate of infection defined for a ∈ [θ, a d ], t ∈ [−θ, T −θ]. We impose the following restrictions on the density-dependent HPV death rate and cell's birth and death rates [9]: 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.

SOLUTION ON REFINED GRIDS
The convergence of approximate solution is studied numerically by the method of refined grids [19], [28]. In the first part of our work [9], we considered the convergence of approximate solution to the benchmark solution in numerical experiments for the model with the specific set of coefficients and initial values. Using the method of refined grids, we study the convergence of numerical approximation for the arbitrary coefficients and initial values of model that allows us for considering more realistic, biologically motivated birth and death rates of susceptible, infectious, precancerous and cancer cells. In experiments we use the constants from Table 1 [9] and the following coefficients and initial values: We evaluate the residual of two solutionsỸ j i andỸ j i computed on two different grids with parameters (h, i = 0, .., N , j = 0, .., M ) and (0.5h, i = 0, .., 2N , j = 0, .., 2M ), respectively. The difference between previousỸ j i and current Y j i solutions (Y j i = (S j i , I j i , P j i , C j i )) and betweeñ V j andṼ j , respectively, is estimated by means of the following dimensionless functionalities: The results of numerical experiments for the different pairs of dimensionless parameter x = h/a d are presented in the Table.
the same as shown in Table 2 in work [9]. The biggest among all subclasses values of difference between two numerical solutionsỸ j i andỸ j i are obtained for the cancer cells subpopulation. This is a consequence of the large numerical error of approximate density of cancer cells subpopulation observed in the experiments in [9]. The results of numerical experiments presented in Table 1 are shown in Fig. 1 ( -the values of δ Table 1 depending from x = h/a d , ----the graphs of corresponding regressions y(x) (with denoted equations)), Fig. 2 C (C) and y(x)) and C (V ) and y(x)). We draw the plots of regressions for all data points obtained in experiments which provide the experimental evaluation of the convergence rate of approximate solutions ( Fig. 1 -Fig. 5). The equations of regressions in all experiments illustrate the stability of numerical method and convergence of numerical approximation with the second order of accuracy. This conclusion is in good agreement with the results of numerical experiments obtained in the first part of our study for the bench-mark solution [9].
In the second group of numerical experiments we study the unstable dynamical regimes of SIPCV population. Simulations reveal existence of two types of such dynamics interesting from the biological point of view. The dynamical regime of the first type for two different initial values of ϕ(a) is shown in Figs. 6 -10. In particular, in Figs. 6, 7, 8, 10 it is shown the oscillating, asymptotically stable dynamics of susceptible, infected and precancerous cells quantity, and HPV quantity, while in Fig. 9 it is shown the unstable dynamics of cancer cells quantity with exponential infinite Vitalii V. Akimenko, Fajar Adi-Kusumo, Age-structured Delayed SIPCV Epidemic Model of HPV ...     Graphs of NI (t) for oscillating dynamics for 2 different initial values ϕ(a).
growth. Such behavior of cancer cells is typical for the population with a basic reproduction number greater than one [3], [11], [16], [21], [22]. From the biological point of view unstable dynamics of cancer cells means that cancer tumor is not localized in biological tissue and tends to form the metastases.
The second type of unstable dynamical regimes corresponds to the asymptotically stable oscillating dynamics of susceptible, infected cells quantity and HPV quantity like those shown in Figs. 6, 7, 10, and unstable dynamics of dysplasia (precancerous cells) and cancer cells quantity shown in Figs. 11, 12 for 2 different initial values of ϕ(a). In this case unstable behavior of cancer cell subpopula-    Graphs of NP (t) for unstable dynamics for 2 different initial values ϕ(a).

Fig. 12.
Graphs of NC (t) for unstable dynamics for 2 different initial values ϕ(a). tion can be induced by the un-limited growth of dysplasia but not by the rapid proliferation and/or low mortality of can-cer cells. Thus, the basic reproduction number of cancer cell subpopulation can be less than one, and their unstable dynamics may be induced by the precancerous cells which play a role of donor for cancer cells in this case.
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.

V. CONCLUSIONS AND DISCUSSION
In this paper we study both theoretically and numerically the convergence of approximate solution of a numerical method designed for the epidemic model of age-structured sub-populations of susceptible, infectious, precancerous and cancer cells and un-structured population of human papilloma virus (HPV) (SIPCV epidemic model) and presented in the first part of our work [9]. Theoretical analysis based on the evaluation of residual of approximate solution shows the stability and converges of the numerical solution to the exact solution of system with the second order of accuracy. The numerical experiments confirm and illustrate the second rate of convergence of numerical approximation on the refined grids. Although this result was expected, since we have obtained the second rate of convergence of approximate solution to the benchmark solution of system in the numerical experiments in work [9], the study of convergence of numerical approximation on the refined grids allows us to consider the model with more realistic biologically motivated coefficients. The results of numerical experiments are in good agreement with the results obtained in [9] and confirm conclusion 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 the large value of session time of simulation.
The numerical experiments with model parameters revealed two types of unstable dynamical regimes of SPICV population. In the first case we observed the unlimited exponential growth of cancer cells population with a basic reproduction number greater than one when the other cells subpopulations and HPV population showed the asymptotically stable oscillating dynamics. In the second case the unlimited exponential growth of dysplasia (precancer cells) induced the unstable dynamics of cancer cells population while the other cells subpopulations and HPV showed the asymptotically stable oscillating dynamics. In both cases the unlimited cancer cells population dynamics tend to the tumor growth and formation of metastases in organism. Overall, in our papers (first and second parts of our project) we developed an age-structured SIPCV epidemic model of HPV and cervical cancer cells dynamics, obtained an exact solution of this model, designed and studied the corresponding numerical method for the exact solution for simulation the dynamics of susceptible, infected, precancerous, cancerous cells and HPV populations with high accuracy, important for study the human papilloma cancer diseases in theoretical biology and medicine.

ACKNOWLEDGMENT
We would like to acknowledge Department of Mathematics, Universitas Gadjah Mada for the organizational support of these studies. Using Eqs. (40), (67), (68) from [9] we can estimate the residual of susceptible cells density (ψ S ) j i = S(a i , t j )−S j i . We use here the formula for numerical error of trapezoidal rule of integration from [17], the Taylor remainder of zero order for the approximation of natural exponen-tial function and the norm of functional space |f (x, y)|. The numerical formula of derivation of residual estimation is given in Appendix A. For k = 1, 0 < t j < a i ≤ a d , 0 < t j ≤ θ < a r we have: , ∂α(a, t) ∂a where A 0 is a real number, A S > 0 is a constant, the value of which depends from the norm of coefficients and norm of their partial derivatives. Expression means the maximum of absolute value of sum of ∂f (ξ l ) ∂ξ by index l at the intervals ξ l ∈ [t l , t l+1 ], i.e. we consider the set of points {ξ l }, ξ l ∈ [t l , t l+1 ], for which the absolute value of this sum takes the maximum value. Existence of such set of points {ξ l } at the intervals ξ l ∈ [t l , t l+1 ] follows from the continuous differentiability of function f (ξ) at the interval ξ ∈ [0, T ]. For k = 1, 0 < a i ≤ t j ≤ θ < a r , we have: , ∂β s (a, t) ∂a where A 1 , A 2 are real constants, A S > 0 is a constant, the value of which depends from the norm of coefficients and norm of their partial derivatives.

APPENDIX B. ESTIMATION OF THE RESIDUAL OF INFECTIOUS CELLS SUBCLASS DENSITY.
From Eqs. (33), (78) from [9] we obtain the residual of infectious cells density , ∂d q (a, t) ∂t , ∂α(a, t) ∂a where ξ l = lh; A 3 , A 4 , A 5 are real constants, A (0) I > 0 is a constant, the value of which depends from the norm of coefficients and norm of their partial derivatives. We use here the formula for the numerical error of trapezoidal rule of integration from [17], the Taylor remainder of zero order for the approximation of natural exponential function and the norm of functional space C.
where ξ l = lh, a m = mh; A 6 . . . A 14 are real constants, A I > 0 is a constant, the value of which depends from the norm of coefficients and norm of their partial derivatives.

APPENDIX C. ESTIMATION OF THE RESIDUAL OF PRECANCEROUS CELLS SUBCLASS DENSITY.
By analogy with Eq. (35), from Eqs. (43), (91), (92) from [9] we obtain the residual of precancerous cells , ∂d q (a, t) ∂a , ∂d r (a, t) ∂a , ∂α(a, t) ∂a where where η p = ph, η q = qh, ξ l = lh; A 15 , ..., A 19 are real constants, A P > 0 is a constant, the value of which depends from the norm of coefficients and norm of their partial derivatives.
Continuing the estimation of residuals of susceptible, infectious, precancerous, cancer cells density and virus density at the next time intervals for 1 < k ≤ K, θ ≤ t j ≤ T , we obtain the evaluation of residuals at the final instant of time T . Finally, for the norm of residuals where the positive constants A * S , A * I , A * P , A * C , A * V depend from the θ M θ , (a m − a r ) M , a M d , (a g − a c ) M , T , norm of coefficients of system (1)- (11) and norm of their partial derivatives of the first and second order.

APPENDIX F. AUXILIARY NUMERICAL FORMULAE FOR INTEGRALS.
For evaluation of the residual of approximate solution we use the formula for the numerical error of trapezoidal rule of integration from [17], the Taylor remainder of zero order for the approximation of natural exponential function and the norm of functional space C. For the expressions with 1-d integral we use the approach: Vitalii V. Akimenko, Fajar Adi-Kusumo, Age-structured Delayed SIPCV Epidemic Model of HPV ...
where x i = ih, h = (b(t j ) − a(t j ))/N , x ia = a(t j ), x i b = b(t j ); A 0 and υ are real constants, ξ j is some point of interval [x j , x j+1 ] given in the formula of numerical error of trapezoidal rule [17], T z (i a , i b , f (x i , t j ), x i )is a trapezoidal rule for the integral b(tj ) a(tj ) f (x, t j )dx given by Eq. (37), A f is a positive constant the value of which does not depend from h and depends from the norm of functions and norm of their partial derivatives of the second order. For the expressions with multiple integrals we use the same approach: where x i = ih, h = (b(t j ) − a(t j ))/N , x ia = a(t j ), x i b = b(t j ), η l = lh; B 0 , B 1 and υ are real constants, ξ j ∈ [x j , x j+1 ] and η l ∈ [x l , x l+1 ] are given constants in the formula of numerical error of trapezoidal rule [17], T z (i a , i b , f (x i , t j ), x i ) -is a trapezoidal rule given by Eq. (37), A f is a positive constant the value of which does not depend from h and depends from the norm of functions and norm of their partial derivatives of the first and second order.