Well-Posedness and a Finite Difference Approximation for a Mathematical Model of HPV-Induced Cervical Cancer

Baoling Ma¹, Jeremy J. Thibodeaux²

¹Department of Mathematics, Millersville University, Millersville, PA, USA.
²Department of Mathematical Sciences, University of the Virgin Islands, St. Croix, USVI, USA.

DOI: 10.4236/am.2023.143009                

Abstract

We present a first-order finite difference scheme for approximating solutions of a mathematical model of cervical cancer induced by the human papillomavirus (HPV), which consists of four nonlinear partial differential equations and a nonlinear first-order ordinary differential equation. The scheme is analyzed and used to provide an existence-uniqueness result. Numerical simulations are performed in order to demonstrate the first-order rate of convergence. A sensitivity analysis was done in order to compare the effects of two drug types, those that increase the death rate of HPV-infected cells, and those that increase the death rate of the precancerous cell population. The model predicts that treatments that affect the precancerous cell population by directly increasing the corresponding death rate are far more effective than those that increase the death rate of HPV-infected cells.

Keywords

Age Structured, Cervical Cancer, Finite Difference, Human Papillomavirus, Mathematical Model

1. Introduction

Human pappilomavirus (HPV) is one of the most widespread sexually transmitted diseases worldwide. New infections in the United States alone are estimated to be between 1 and 5.5 million per year [1] . Over 200 types have been genetically distinguished.

In females, the virus targets basal epithelial cells in the cervix. With prolonged exposure, these infections can lead to the development of precancerous cells and, without treatment, cancerous cells [2] . Of the over 200 distinguished gentotypes, HPV types 16, 18, 31, and 45 are considered to be high-risk for inducing cervical cancer after prolonged infection.

Mathematical modeling techniques have been used to investigate the dynamics of cervical cancer, both from epidemiological [3] [4] [5] [6] and cellular [2] [7] [8] [9] [10] viewpoints. Our work here is in the latter context. The initial studies of within-host HPV infection leading to cervical cancer described the dynamics with systems of ordinary differential equations [2] [8] [10] . The models included the dynamics of the susceptible epithelial cells, infected cells, precancerous cells, cancerous cells, and viral particles. In [7] , the authors modified the previous models to incorporate the age of each of the cell types previously mentioned. This results in a more extensive description of the biological processes, with the added cost of a more complicated model. Taking the age of cells into account yields a model of a system of four first-order, nonlinear, partial differential equations and an ordinary differential equation. More recently, Sari et al. [9] proposed a modification of the model in [7] that eliminates the ordinary differential equation describing the viral population.

Systems of first-order partial differential equations have been applied to a wide range of biological processes, including algal growth, amphibian population growth, erythropoiesis, and fungal propagation [11] [12] [13] [14] [15] . The mathematical approaches to analyzing those systems have been just as varied. In some cases, such as in [7] [9] , stability analysis can be successfully applied to investigate the long-term behavior of solutions. If more knowledge of transient solution behavior or sensitivity analysis is desired, then numerical methods provide a useful tool.

Finite element methods, integrating along characteristics, upper-lower solution techniques, and finite difference schemes have all been successfully applied to various systems of first-order partial differential equations [11] [14] [16] [17] [18] . In this work, we follow the ideas presented in [18] , which were successfully applied to the nonlinear size-structured population model in [16] . The approach yields a simple first-order finite difference scheme that can first be used to provide existence-uniqueness results, and then can be implemented to approximate solutions of the model.

The paper is organized as follows. In Section 2, we define the various components of the model and the paraments that affect their behavior. In Section 3, we define weak solutions of the model and provide conditions on parameters that are required for the results that are presented in later sections. We also define the finite difference scheme in detail and the notation that will be used throughout the paper. In Section 4, we will provide the results necessary to establish an existence-uniqueness result. In Section 5, we present the main theoretical results of this work. Namely, the construction of sequences of functions that converge to the unique weak solution of the model. To our knowledge, this is the first well-posedness result for the model. In Section 6, we numerically analyze the finite difference scheme itself and show that it indeed provides first-order accuracy. In Section 7, we explore the model’s sensitivity to certain parameters that are associated with currently available drug treatments. Finally, in Section 8 we will make some concluding remarks.

2. The Model

The model that is the focus of this work was developed in [7] and takes the form of a system of four first-order partial differential equations and an ordinary differential equation. The model describes the population densities of the susceptible cell ( $S\left(a\mathrm{,}t\right)$ ), infected cell ( $I\left(a\mathrm{,}t\right)$ ), precancerous cell ( $P\left(a\mathrm{,}t\right)$ ), and cancerous cell ( $C\left(a\mathrm{,}t\right)$ ) populations, as well as the number of HPV viral particles ( $V\left(t\right)$ ), at time t and age a. The total populations of the susceptible, infected, precancerous, and cancerous cell types at time t are denoted by $N_S\left(t\right)={\displaystyle {\int }_0^{a_d}}S\left(a\mathrm{,}t\right)\text{d}a$ , $N_I\left(t\right)={\displaystyle {\int }_0^{a_d}}I\left(a\mathrm{,}t\right)\text{d}a$ , $N_P\left(t\right)={\displaystyle {\int }_0^{a_d}}P\left(a\mathrm{,}t\right)\text{d}a$ , and $N_C\left(t\right)={\displaystyle {\int }_0^{a_d}}C\left(a\mathrm{,}t\right)\text{d}a$ , respectively, where $a_d$ is the maximum lifespan of the cells. We denote their respective death rates by $d_s\mathrm{,}{d}_i\mathrm{,}{d}_p$ , and $d_c$ . The infection of healthy cells by free virus particles is modeled by the mass action term $\alpha V\left(t\right)S\left(a\mathrm{,}t\right)$ . Infected cells are partially moved into the precancerous cell population, modeled by the term $\delta I\left(t\right)$ . The rate at which precancerous cells move into the cancerous cell population is modeled dynamically with the bounded function $\mu \left(N_P\right)=\frac{\theta N_P\left(t\right)}{1+\kappa N_P\left(t\right)}$ , although the theoretical results we provide in this work are valid for any continuously differentiable function $\mu$ with $0<\mu \le 1$ . The viral population is modeled by means of an ordinary differential equation. It is assumed that there is a constant production of viral particles, at a rate of $\Lambda$ . There is a density-dependent death rate $d_v\left(V\right)$ . Since the total population of infected cells dies at a rate of $d_i{N}_I\left(t\right)$ , it is assumed that viral particles are produced at a rate proportional to this rate, with a constant of proportionality n. Non-cancerous cells are assumed to have a reproductive age range of $\left[a_r\mathrm{,}{a}_q\right]$ , while that of cancerous cells is $\left[a_c\mathrm{,}{a}_l\right]$ . The reproductive age range of cancer cells is shifted relative to that of the noncancerous cells [7] , and so it is assumed that $a_c<a_r$ and $a_l<a_q$ . Noncancerous cells are assumed to have a birth rate of ${\beta }_0$ , while cancerous cells have a birth rate of ${\beta }_c$ . Finally, the initial cell densities are given by $S\left(a,0\right)=S_0\left(a\right)$ , $I\left(a,0\right)=I_0\left(a\right)$ , $P\left(a,0\right)=P_0\left(a\right)$ , and $C\left(a\mathrm{,0}\right)=C_0\left(a\right)$ , and the initial viral load is $V\left(0\right)=V_0$ .

These definitions and assumptions lead to the following model:

$\begin{array}{ll}\frac{\partial S\left(a,t\right)}{\partial t}+\frac{\partial S\left(a,t\right)}{\partial a}=-\left(d_s+\alpha V\left(t\right)\right)S\left(a,t\right),\hfill & 0<t<T,0<a<a_d,\hfill \\ \frac{\partial I\left(a,t\right)}{\partial t}+\frac{\partial I\left(a,t\right)}{\partial a}=-\left(d_i+\delta \right)I\left(a,t\right)+\alpha V\left(t\right)S\left(a,t\right),\hfill & 0<t<T,0<a<a_d,\hfill \\ \frac{\partial P\left(a,t\right)}{\partial t}+\frac{\partial P\left(a,t\right)}{\partial a}=-d_{p}P\left(a,t\right)+\delta I\left(a,t\right)-\mu \left(N_p\right)P\left(a,t\right),\hfill & 0<t<T,0<a<a_d,\hfill \\ \frac{\partial C\left(a,t\right)}{\partial t}+\frac{\partial C\left(a,t\right)}{\partial a}=-d_{c}C\left(a,t\right)+\mu \left(N_p\right)P\left(a,t\right),\hfill & 0<t<T,0<a<a_d,\hfill \\ \frac{\text{d}V\left(t\right)}{\text{d}t}=\Lambda -d_v\left(V\right)V\left(t\right)+n{d}_i{\displaystyle {\int }_0^{a_d}}I\left(a,t\right)\text{d}a,\hfill & 0<t<T,\hfill \end{array}$

$\begin{array}{ll}S\left(0,t\right)={\beta }_0{\displaystyle {\int }_{a_r}^{a_q}}S\left(a,t\right)\text{d}a,\hfill & 0<t<T,\hfill \\ I\left(0,t\right)={\beta }_0{\displaystyle {\int }_{a_r}^{a_q}}I\left(a,t\right)\text{d}a,\hfill & 0<t<T,\hfill \\ P\left(0,t\right)={\beta }_0{\displaystyle {\int }_{a_r}^{a_q}}P\left(a,t\right)\text{d}a,\hfill & 0<t<T,\hfill \\ C\left(0,t\right)={\beta }_c{\displaystyle {\int }_{a_c}^{a_l}}C\left(a,t\right)\text{d}a,\hfill & 0<t<T,\hfill \end{array}$

$\begin{array}{ll}S\left(a,0\right)=S_0\left(a\right),\hfill & 0<a<a_d,\hfill \\ I\left(a,0\right)=I_0\left(a\right),\hfill & 0<a<a_d,\hfill \\ P\left(a,0\right)=P_0\left(a\right),\hfill & 0<a<a_d,\hfill \\ C\left(a,0\right)=C_0\left(a\right),\hfill & 0<a<a_d,\hfill \\ V\left(0\right)=V_0\hfill & \hfill \end{array}$ (2.1)

An uncommon aspect of this model is that the integrals in the boundary conditions are not over the entire interval $\left[\mathrm{0,}{a}_d\right]$ . This presents difficulties in both the implementation of the finite difference scheme and the proofs that will follow in later sections.

Remark 2.1 While the work in [9] proposes a modified model that eliminates the ordinary differential equation in (2.1), the ideas of the forthcoming numerical scheme can be adapted for such a case. Hence, we preform the analysis on the model proposed in [7] .

3. A Finite Difference Scheme

Let c be a sufficiently large positive constant. Throughout the discussion we impose the following regularity conditions on our model parameters in (2.1).

1) The function $d_v\mathrm{<mo>\; :\; </mo>}\left[\mathrm{0,}\infty \right)\to \left[\mathrm{0,}\infty \right)$ is continuous.

2) The function $\mu \mathrm{<mo>\; :\; </mo>}\left[\mathrm{0,}\infty \right)\to \left[\mathrm{0,}\infty \right)$ is continuous and satisfies $0<\mu \left(x\right)\le 1$ for all x in $\left[\mathrm{0,}\infty \right)$ .

3) The parameters $V_0\mathrm{,}{d}_s\mathrm{,}{d}_i\mathrm{,}{d}_p\mathrm{,}\alpha \mathrm{,}\delta \mathrm{,}n\mathrm{,}{a}_l\mathrm{,}{a}_c\mathrm{,}{a}_r\mathrm{,}{a}_q\mathrm{,}{\beta }_0\mathrm{,}{\beta }_c$ are nonnegative constants.

4) The initial conditions $S_0\left(a\right)\mathrm{,}{I}_0\left(a\right)\mathrm{,}{P}_0\left(a\right)\mathrm{,}{C}_0\left(a\right)$ are nonnegative with

${\Vert S_0\Vert }_{BV\left(\mathrm{0,}{a}_d\right)}\mathrm{,}{\Vert I^0\Vert }_{BV\left(\mathrm{0,}{a}_d\right)}\mathrm{,}{\Vert P_0\Vert }_{BV\left(\mathrm{0,}{a}_d\right)}\mathrm{,}{\Vert C_0\Vert }_{BV\left(\mathrm{0,}{a}_d\right)}\le c$

Now we give the definition of a weak solution to problem (2.1) as follows: A 5-tuple $\left(S\mathrm{,}I\mathrm{,}P\mathrm{,}C\mathrm{,}V\right)\in (BV\left(\left(\mathrm{0,}{a}_d\right)\times \left(\mathrm{0,}T\right)\right)\mathrm{,}BV\left(\left(\mathrm{0,}{a}_d\right)\times \left(\mathrm{0,}T\right)\right)\mathrm{,}BV\left(\left(\mathrm{0,}{a}_d\right)\times \left(\mathrm{0,}T\right)\right)\mathrm{,}$ $BV\left(\left(\mathrm{0,}{a}_d\right)\times \left(\mathrm{0,}T\right)\right)\mathrm{,}C\left[\mathrm{0,}T\right])$ is called a weak solution to the problem (2.1) if it satisfies:

$\begin{array}{l}{\displaystyle {\int }_0^{a_d}}S\left(\nu ,t\right){\zeta }_1\left(\nu ,t\right)\text{d}\nu \\ ={\displaystyle {\int }_0^{a_d}}{S}_0\left(\nu \right){\zeta }_1\left(\nu ,0\right)\text{d}\nu -{\displaystyle {\int }_0^t}S\left(a_d^{-},\tau \right){\zeta }_1\left(a_d,\tau \right)\text{d}\tau \\ +{\beta }_0{\displaystyle {\int }_0^t}\left({\displaystyle {\int }_{a_r}^{a_q}}S\left(a,\tau \right)\text{d}a\right){\zeta }_1\left(0,\tau \right)\text{d}\tau \\ +{\displaystyle {\int }_0^t}{\displaystyle {\int }_0^{a_d}}\left(\frac{\partial {\zeta }_1}{\partial \tau }\left(\nu ,\tau \right)+\frac{\partial {\zeta }_1}{\partial \nu }\left(\nu ,\tau \right)\right)S\left(\nu ,\tau \right)\text{d}\nu \text{d}\tau \\ -{\displaystyle {\int }_0^t}{\displaystyle {\int }_0^{a_d}}\left(d_s+\alpha V\left(\tau \right)\right)S\left(\nu ,\tau \right){\zeta }_1\left(\nu ,\tau \right)\text{d}\nu \text{d}\tau ,\end{array}$

$\begin{array}{l}{\displaystyle {\int }_0^{a_d}}I\left(\nu ,t\right){\zeta }_2\left(\nu ,t\right)\text{d}\nu \\ ={\displaystyle {\int }_0^{a_d}}{I}_0\left(\nu \right){\zeta }_2\left(\nu ,0\right)\text{d}\nu -{\displaystyle {\int }_0^t}I\left(a_d^{-},\tau \right){\zeta }_2\left(a_d,\tau \right)\text{d}\tau \\ +{\beta }_0{\displaystyle {\int }_0^t}\left({\displaystyle {\int }_{a_r}^{a_q}}I\left(a,\tau \right)\text{d}a\right){\zeta }_2\left(0,\tau \right)\text{d}\tau \\ +{\displaystyle {\int }_0^t}{\displaystyle {\int }_0^{a_d}}\left(\frac{\partial {\zeta }_2}{\partial \tau }\left(\nu ,\tau \right)+\frac{\partial {\zeta }_2}{\partial \nu }\left(\nu ,\tau \right)\right)I\left(\nu ,\tau \right)\text{d}\nu \text{d}\tau \\ +{\displaystyle {\int }_0^t}{\displaystyle {\int }_0^{a_d}}\left(-\left(d_i+\delta \right)I\left(\nu ,\tau \right)+\alpha V\left(t\right)S\left(\nu ,\tau \right)\right){\zeta }_2\left(\nu ,\tau \right)\text{d}\nu \text{d}\tau ,\end{array}$

$\begin{array}{l}{\displaystyle {\int }_0^{a_d}}P\left(\nu ,t\right){\zeta }_3\left(\nu ,t\right)\text{d}\nu \\ ={\displaystyle {\int }_0^{a_d}}{P}_0\left(\nu \right){\zeta }_3\left(\nu ,0\right)\text{d}\nu -{\displaystyle {\int }_0^t}P\left(a_d^{-},\tau \right){\zeta }_3\left(a_d,\tau \right)\text{d}\tau \\ +{\beta }_0{\displaystyle {\int }_0^t}\left({\displaystyle {\int }_{a_r}^{a_q}}P\left(a,\tau \right)\text{d}a\right){\zeta }_3\left(0,\tau \right)\text{d}\tau \\ +{\displaystyle {\int }_0^t}{\displaystyle {\int }_0^{a_d}}\left(\frac{\partial {\zeta }_3}{\partial \tau }\left(\nu ,\tau \right)+\frac{\partial {\zeta }_3}{\partial \nu }\left(\nu ,\tau \right)\right)P\left(\nu ,\tau \right)\text{d}\nu \text{d}\tau \\ +{\displaystyle {\int }_0^t}{\displaystyle {\int }_0^{a_d}}\left(-d_{p}P\left(\nu ,\tau \right)+\delta I\left(\nu ,\tau \right)-\mu \left(N_p\right)P\left(\nu ,\tau \right)\right){\zeta }_3\left(\nu ,\tau \right)\text{d}\nu \text{d}\tau ,\end{array}$

$\begin{array}{l}{\displaystyle {\int }_0^{a_d}}C\left(\nu ,t\right){\zeta }_4\left(\nu ,t\right)\text{d}\nu \\ ={\displaystyle {\int }_0^{a_d}}{C}_0\left(\nu \right){\zeta }_4\left(\nu ,0\right)\text{d}\nu -{\displaystyle {\int }_0^t}C\left(a_d^{-},\tau \right){\zeta }_4\left(a_d,\tau \right)\text{d}\tau \\ +{\beta }_c{\displaystyle {\int }_0^t}\left({\displaystyle {\int }_{a_r}^{a_q}}C\left(a,\tau \right)\text{d}a\right){\zeta }_4\left(0,\tau \right)\text{d}\tau \\ +{\displaystyle {\int }_0^t}{\displaystyle {\int }_0^{a_d}}\left(\frac{\partial {\zeta }_4}{\partial \tau }\left(\nu ,\tau \right)+\frac{\partial {\zeta }_4}{\partial \nu }\left(\nu ,\tau \right)\right)C\left(\nu ,\tau \right)\text{d}\nu \text{d}\tau \\ +{\displaystyle {\int }_0^t}{\displaystyle {\int }_0^{a_d}}\left(-d_{c}C\left(\nu ,\tau \right)+\mu \left(N_p\right)P\left(\nu ,\tau \right)\right){\zeta }_4\left(\nu ,\tau \right)\text{d}\nu \text{d}\tau ,\end{array}$

$V\left(t\right)=V_0+{\displaystyle {\int }_0^t}\left(\Lambda -d_v\left(V\right)V\left(\tau \right)+n{d}_i{\displaystyle {\int }_0^{a_d}}I\left(a,\tau \right)\text{d}a\right)\text{d}\tau ,$

for each $t\in \left(\mathrm{0,}T\right)$ , every test function ${\zeta }_i\in C^1\left(\left[0,a_d\right]\times \left[0,T\right]\right),i=1,2,3,4$ . Here, $S\left(a_d^{-}\mathrm{,}\tau \right)={\lim }_{\nu \to a_d}S\left(\nu \mathrm{,}\tau \right)$ , $I\left(a_d^{-}\mathrm{,}\tau \right)={\lim }_{\nu \to a_d}I\left(\nu \mathrm{,}\tau \right)$ , $P\left(a_d^{-}\mathrm{,}\tau \right)={\lim }_{\nu \to a_d}P\left(\nu \mathrm{,}\tau \right)$ , and $C\left(a_d^{-}\mathrm{,}\tau \right)={\lim }_{\nu \to a_d}C\left(\nu \mathrm{,}\tau \right)$ . First, we divide the intervals $\left[\mathrm{0,}{a}_d\right]$ and $\left[\mathrm{0,}T\right]$ into m and K subintervals, respectively, and use the following notation throughout the paper: $\Delta a=a_d/m$ , and $\Delta t=T/K$ denote the age and time mesh sizes, respectively. The mesh points are given by: $a_j=j\Delta a$ , $j=0,1,\cdots ,m$ and $t_k=k\Delta t$ , $k=0,1,\cdots ,K$ . We denote by $S_j^k$ , $I_j^k$ , $P_j^k$ , and $C_j^k$ , the finite difference approximations of $S\left(a_j\mathrm{,}{t}_k\right)$ , $I\left(a_j\mathrm{,}{t}_k\right)$ , $P\left(a_j\mathrm{,}{t}_k\right)$ , and $C\left(a_j\mathrm{,}{t}_k\right)$ , respectively. Further, we denote $V\left(t_k\right)$ by $V^k$ , $\mu \left(N_P\left(t_k\right)\right)$ by ${\mu }_P^k$ , and $d_v\left(V^k\right)$ by $d_v^k$ . We approximate the total cell populations with

$N_X^k={\displaystyle \underset{l=1}{\overset{m}{\sum }}}{X}_l^k\Delta a,$

where $X=S\mathrm{,}I\mathrm{,}P\mathrm{,}C$ . Since the parameters $a_r\mathrm{,}{a}_q\mathrm{,}{a}_c$ and $a_l$ will not generally coincide with our mesh points, we define the following:

$\begin{array}{l}{j}_r^{+}=\min \left\{j|a_j>a_r\right\},j_q^{-}=\max \left\{j|a_j<a_q\right\},\\ j_c^{+}=\min \left\{j|a_j>a_c\right\},j_l^{-}=\max \left\{j|a_j<a_l\right\}.\end{array}$

The total variation, $l_1$ norm, and $l_{\infty }$ norm are defined as

$TV\left(u^k\right)={\displaystyle \underset{j=0}{\overset{m-1}{\sum }}}\left|u_{j+1}^k-u_j^k\right|,$

${\Vert u^k\Vert }_1={\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|u_j^k\right|,{\Vert u^k\Vert }_{\infty }=\underset{0\le j\le m}{\max }\left|u_j^k\right|,$

respectively. The first-order scheme that will be considered here is given by:

$\begin{array}{ll}\frac{S_j^{k+1}-S_j^k}{\Delta t}+\frac{S_j^k-S_{j-1}^k}{\Delta a}=-\left(d_s+\alpha V^k\right)S_j^{k+1},\hfill & 1\le j\le m,0\le k\le K-1\hfill \\ \frac{I_j^{k+1}-I_j^k}{\Delta t}+\frac{I_j^k-I_{j-1}^k}{\Delta a}=-\left(d_i+\delta \right)I_j^{k+1}+\alpha V^k{S}_j^{k+1},\hfill & 1\le j\le m,0\le k\le K-1\hfill \\ \frac{P_j^{k+1}-P_j^k}{\Delta t}+\frac{P_j^k-P_{j-1}^k}{\Delta a}=-d_p{P}_j^{k+1}+\delta I_j^{k+1}-{\mu }_P^k{P}_j^{k+1},\hfill & 1\le j\le m,0\le k\le K-1\hfill \\ \frac{C_j^{k+1}-C_j^k}{\Delta t}+\frac{C_j^k-C_{j-1}^k}{\Delta a}=-d_c{C}_j^{k+1}+{\mu }_P^k{P}_j^{k+1},\hfill & 1\le j\le m,0\le k\le K-1\hfill \\ \frac{V^{k+1}-V^k}{\Delta t}=\Lambda -d_v^k{V}^{k+1}+n{d}_i{N}_I^k,\hfill & 0\le k\le K-1.\hfill \end{array}$ (3.1)

The boundary conditions are computed as follows:

$\begin{array}{ll}{S}_0^{k+1}={\beta }_0\left(S\left(a_{j_r^{+}},t_{k+1}\right)\left(a_{j_r^{+}}-a_r\right)+{\displaystyle \underset{j=j_r^{+}}{\overset{j_q^{-}}{\sum }}}S\left(a_j,t_{k+1}\right)\Delta a+S\left(a_{j_q^{-}},t_{k+1}\right)\left(a_q-a_{j_q^{-}}\right)\right),\hfill & 0\le k\le K-1\hfill \\ I_0^{k+1}={\beta }_0\left(I\left(a_{j_r^{+}},t_{k+1}\right)\left(a_{j_r^{+}}-a_r\right)+{\displaystyle \underset{j=j_r^{+}}{\overset{j_q^{-}}{\sum }}}I\left(a_j,t_{k+1}\right)\Delta a+I\left(a_{j_q^{-}},t_{k+1}\right)\left(a_q-a_{j_q^{-}}\right)\right),\hfill & 0\le k\le K-1\hfill \\ P_0^{k+1}={\beta }_0\left(P\left(a_{j_r^{+}},t_{k+1}\right)\left(a_{j_r^{+}}-a_r\right)+{\displaystyle \underset{j=j_r^{+}}{\overset{j_q^{-}}{\sum }}}P\left(a_j,t_{k+1}\right)\Delta a+P\left(a_{j_q^{-}},t_{k+1}\right)\left(a_q-a_{j_q^{-}}\right)\right),\hfill & 0\le k\le K-1\hfill \\ C_0^{k+1}={\beta }_c\left(C\left(a_{j_c^{+}},t_{k+1}\right)\left(a_{j_c^{+}}-a_c\right)+{\displaystyle \underset{j=j_c^{+}}{\overset{j_l^{-}}{\sum }}}C\left(a_j,t_{k+1}\right)\Delta a+C\left(a_{j_l^{-}},t_{k+1}\right)\left(a_l-a_{j_l^{-}}\right)\right),\hfill & 0\le k\le K-1.\hfill \end{array}$ (3.2)

The initial conditions for p and m are computed as follows:

$S_i^0=S^0\left(a_j\right),I_i^0=I^0\left(a_j\right),P_i^0=P^0\left(a_j\right),C_i^0=C^0\left(a_j\right),j=0,\cdots ,m$

Equation (3.1) is more conveniently used in the following form:

$\begin{array}{lll}{S}_j^{k+1}=\frac{S_j^k\left(1-\frac{\Delta t}{\Delta a}\right)+\frac{\Delta t}{\Delta a}{S}_{j-1}^k}{1+\Delta t\left(d_s+\alpha V^k\right)},\hfill & 1\le j\le m,\hfill & 0\le k\le K-1\hfill \\ I_j^{k+1}=\frac{I_j^k\left(1-\frac{\Delta t}{\Delta a}\right)+\frac{\Delta t}{\Delta a}{I}_{j-1}^k+\Delta t\alpha V^k{S}_j^{k+1}}{1+\Delta t\left(d_i+\delta \right)},\hfill & 1\le j\le m,\hfill & 0\le k\le K-1\hfill \\ P_j^{k+1}=\frac{P_j^k\left(1-\frac{\Delta t}{\Delta a}\right)+\frac{\Delta t}{\Delta a}{P}_{j-1}^k+\Delta t\delta I_j^{k+1}}{1+\Delta t\left(d_p+{\mu }_P^k\right)},\hfill & 1\le j\le m,\hfill & 0\le k\le K-1\hfill \end{array}$

$\begin{array}{lll}{C}_j^{k+1}=\frac{C_j^k\left(1-\frac{\Delta t}{\Delta a}\right)+\frac{\Delta t}{\Delta a}{C}_{j-1}^k+\Delta t{\mu }_P^k{P}_j^{k+1}}{1+\Delta t{d}_c},\hfill & 1\le j\le m,\hfill & 0\le k\le K-1\hfill \\ V^{k+1}=\frac{V^k+\Delta t\left(\Lambda +n{d}_i{N}_I^k\right)}{1+\Delta t{d}_v^k},\hfill & \hfill & 0\le k\le K-1\hfill \end{array}$ (3.3)

It is clear from the assumptions on the initial conditions and Equation (3.3) that the numerical solutions will remain nonnegative as long as $\frac{\Delta t}{\Delta a}<1$ .

4. Estimates for the Finite Difference Approximations

The lemma below shows that the numerical approximations are bounded in the $l^1$ norm.

Lemma 4.1 There exists a positive constant $M_1$ , independent of $\Delta a$ and $\Delta t$ , such that

${\Vert S^k\Vert }_1+{\Vert I^k\Vert }_1+{\Vert P^k\Vert }_1+{\Vert C^k\Vert }_1+\left|V^k\right|\le M_1\mathrm{,}$

for $k=0,1,\cdots ,K$ .

Proof. Multiply the first four equations of (3.1) on both sides by $\Delta t$ and $\Delta a$ and sum over $j=1,2,\cdots ,m$

$\begin{array}{l}{\Vert S^{k+1}\Vert }_1-{\Vert S^k\Vert }_1+S_m^k\Delta t-S_0^k\Delta t=-\left(d_s+\alpha V^k\right){\Vert S^{k+1}\Vert }_1\Delta t,\\ {\Vert I^{k+1}\Vert }_1-{\Vert I^k\Vert }_1+I_m^k\Delta t-I_0^k\Delta t=-\left(d_i+\delta \right){\Vert I^{k+1}\Vert }_1\Delta t+\alpha V^k{\Vert S^{k+1}\Vert }_1\Delta t,\\ {\Vert P^{k+1}\Vert }_1-{\Vert P^k\Vert }_1+P_m^k\Delta t-P_0^k\Delta t=-d_p{\Vert P^{k+1}\Vert }_1\Delta t+\delta {\Vert I^{k+1}\Vert }_1\Delta t-{\mu }_p^k{\Vert P^{k+1}\Vert }_1\Delta t,\\ {\Vert C^{k+1}\Vert }_1-{\Vert C^k\Vert }_1+C_m^k\Delta t-C_0^k\Delta t=-d_c{\Vert C^{k+1}\Vert }_1\Delta t+{\mu }_p^k{\Vert P^{k+1}\Vert }_1\Delta t\end{array}$ (4.1)

Multiply the fifth equation of (3.1) on both sides by $\Delta t$

$V^{k+1}-V^k=\Lambda \Delta t-d_v^k{V}^{k+1}\Delta t+n{d}_i{\Vert I^k\Vert }_1\Delta t$ (4.2)

Adding up the five equations in (4.1) and (4.2) we have

$\begin{array}{l}\left(1+d_s\Delta t\right){\Vert S^{k+1}\Vert }_1+\left(1+d_i\Delta t\right){\Vert I^{k+1}\Vert }_1+\left(1+d_p\Delta t\right){\Vert P^{k+1}\Vert }_1\\ +\left(1+d_c\Delta t\right){\Vert C^{k+1}\Vert }_1+\left(1+d_v^k\Delta t\right)\left|V^{k+1}\right|\\ ={\Vert S^k\Vert }_1+\left(1+n{d}_i\Delta t\right){\Vert I^k\Vert }_1+{\Vert P^k\Vert }_1+{\Vert C^k\Vert }_1+\left|V^k\right|\\ +\left(S_0^k+I_0^k+P_0^k+C_0^k+\Lambda \right)\Delta t-\left(S_m^k+I_m^k+P_m^k+C_m^k\right)\Delta t\end{array}$

From the boundary condition, we have

$S_0^k\le {\beta }_0{\Vert S^k\Vert }_1\mathrm{,}{I}_0^k\le {\beta }_0{\Vert I^k\Vert }_1\mathrm{,}{P}_0^k\le {\beta }_0{\Vert P^k\Vert }_1\mathrm{,}{C}_0^k\le {\beta }_c{\Vert C^k\Vert }_1$

Therefore,

$\begin{array}{l}\left(1+d_s\Delta t\right){\Vert S^{k+1}\Vert }_1+\left(1+d_i\Delta t\right){\Vert I^{k+1}\Vert }_1+\left(1+d_p\Delta t\right){\Vert P^{k+1}\Vert }_1\\ +\left(1+d_c\Delta t\right){\Vert C^{k+1}\Vert }_1+\left(1+d_v^k\Delta t\right)\left|V^{k+1}\right|\\ \le \left(1+{\beta }_0\Delta t\right)\left({\Vert S^k\Vert }_1+{\Vert P^k\Vert }_1\right)+\left(1+{\beta }_c\Delta t\right){\Vert C^k\Vert }_1\\ +\left(1+{\beta }_0\Delta t+n{d}_i\Delta t\right){\Vert I^k\Vert }_1+\left|V^k\right|+\Lambda \Delta t\end{array}$

Let $\tilde{C}=\max \left\{{\beta }_0+n{d}_i,{\beta }_c\right\}$ . Then

$\begin{array}{l}{\Vert S^{k+1}\Vert }_1+{\Vert I^{k+1}\Vert }_1+{\Vert P^{k+1}\Vert }_1+{\Vert C^{k+1}\Vert }_1+\left|V^{k+1}\right|\\ \le \left(1+\tilde{C}\Delta t\right)\left({\Vert S^k\Vert }_1+{\Vert I^k\Vert }_1+{\Vert P^k\Vert }_1+{\Vert C^k\Vert }_1+\left|V^k\right|\right)+\Lambda \Delta t\mathrm{,}\end{array}$

which implies the result.

Next we establish a bound on the infinity norm of the numerical approximations.

Lemma 4.2 There exists a positive constant $M_2$ , independent of $\Delta a$ and $\Delta t$ , such that

${\Vert S^k\Vert }_{\infty }+{\Vert I^k\Vert }_{\infty }+{\Vert P^k\Vert }_{\infty }+{\Vert C^k\Vert }_{\infty }\le M_2\mathrm{,}$

for $k=0,1,\cdots ,K$ .

Proof. If ${\Vert S^{k+1}\Vert }_{\infty }=S_0^{k+1}$ , then

${\Vert S^{k+1}\Vert }_{\infty }\le {\beta }_0{\Vert S^{k+1}\Vert }_1\le {\beta }_0{M}_1\mathrm{.}$

If ${\Vert S^{k+1}\Vert }_{\infty }$ is not obtained on the boundary, then there exists a $j_0$ with $1\le j_0\le m$ such that ${\Vert S^{k+1}\Vert }_{\infty }=S_{j0}^{k+1}$ .

From the first equation of (3.1)

$S_{j_0}^{k+1}=S_{j_0}^k-\frac{\Delta t}{\Delta a}\left(S_{j_0}^k-S_{j_0-1}^k\right)-\left(d_s+\alpha V^k\right)S_{j_0}^{k+1}\Delta t,$

That is,

$S_{j_0}^{k+1}\left[1+\left(d_s+\alpha V^k\right)\Delta t\right]=S_{j_0}^k\left(1-\frac{\Delta t}{\Delta a}\right)+S_{j_0-1}^k\frac{\Delta t}{\Delta a}.$

Thus,

${\Vert S^{k+1}\Vert }_{\infty }\le {\Vert S^k\Vert }_{\infty }\mathrm{.}$

This implies that ${\Vert S^{k+1}\Vert }_{\infty }$ is bounded since ${\Vert S^0\Vert }_{\infty }$ is bounded. Similarly, if ${\Vert I^{k+1}\Vert }_{\infty }=I_0^{k+1}$ , then

${\Vert I^{k+1}\Vert }_{\infty }\le {\beta }_0{\Vert I^{k+1}\Vert }_1\le {\beta }_0{M}_1\mathrm{.}$

If ${\Vert I^{k+1}\Vert }_{\infty }$ is not obtained on the boundary, then there exists a $j_0$ with $1\le j_0\le m$ such that ${\Vert I^{k+1}\Vert }_{\infty }=I_{j0}^{k+1}$ .

From the second equation of (3.1)

$I_{j_0}^{k+1}=I_{j_0}^k-\frac{\Delta t}{\Delta a}\left(I_{j_0}^k-I_{j_0-1}^k\right)-\left(d_i+\delta \right)I_{j_0}^{k+1}\Delta t+\alpha V^k{S}_{j_0}^{k+1}\Delta t.$

Therefore,

$\left(1+\left(d_i+\delta \right)\Delta t\right){\Vert I^{k+1}\Vert }_{\infty }\le {\Vert I^k\Vert }_{\infty }+\alpha V^k{\Vert S^{k+1}\Vert }_{\infty }\Delta t\mathrm{.}$

${\Vert I^{k+1}\Vert }_{\infty }\le {\Vert I^k\Vert }_{\infty }+\alpha M_1{\Vert S^{k+1}\Vert }_{\infty }\Delta t\mathrm{.}$

Similarly, if ${\Vert P^{k+1}\Vert }_{\infty }=P_0^{k+1}$ , then

${\Vert P^{k+1}\Vert }_{\infty }\le {\beta }_0{\Vert P^{k+1}\Vert }_1\le {\beta }_0{M}_1\mathrm{.}$

If ${\Vert P^{k+1}\Vert }_{\infty }$ is not obtained on the boundary, then there exists a $j_0$ with $1\le j_0\le m$ such that ${\Vert P^{k+1}\Vert }_{\infty }=P_{j_0}^{k+1}$ .

From the third equation of (3.1)

$P_{j_0}^{k+1}\left(1+\Delta t\left(d_p+{\mu }_P^k\right)\right)=\left(1-\frac{\Delta t}{\Delta a}\right)P_{j_0}^k+\frac{\Delta t}{\Delta a}{P}_{j_0-1}^k+\delta I_{j_0}^{k+1}\Delta t\mathrm{.}$

Therefore,

${\Vert P^{k+1}\Vert }_{\infty }\le {\Vert P^k\Vert }_{\infty }+\delta {\Vert I^{k+1}\Vert }_{\infty }\Delta t\mathrm{.}$

Similarly, if ${\Vert C^{k+1}\Vert }_{\infty }=C_0^{k+1}$ , then

${\Vert C^{k+1}\Vert }_{\infty }\le {\beta }_c{\Vert C^{k+1}\Vert }_1\le {\beta }_c{M}_1\mathrm{.}$

If ${\Vert C^{k+1}\Vert }_{\infty }$ is not obtained on the boundary, then there exists a $j_0$ with $1\le j_0\le m$ such that ${\Vert C^{k+1}\Vert }_{\infty }=C_{j_0}^{k+1}$ .

From the third equation of (3.1)

$C_{j_0}^{k+1}\left(1+d_c\Delta t\right)=\left(1-\frac{\Delta t}{\Delta a}\right)C_{j_0}^k+\frac{\Delta t}{\Delta a}{C}_{j_0-1}^k+{\mu }_p^k{P}_{j_0}^{k+1}\Delta t\mathrm{.}$

Therefore,

${\Vert C^{k+1}\Vert }_{\infty }\le {\Vert C^k\Vert }_{\infty }+{\mu }_p^k{\Vert P^{k+1}\Vert }_{\infty }\Delta t\le {\Vert C^k\Vert }_{\infty }+{\Vert P^{k+1}\Vert }_{\infty }\Delta t\mathrm{.}$

In the next lemma we show that the difference approximations are of bounded total variation.

Lemma 4.3 There exists a positive constant $M_3$ , independent of $\Delta a$ and $\Delta t$ , such that

$TV\left(S^k\right)+TV\left(I^k\right)+TV\left(P^k\right)+TV\left(C^k\right)\le M_3$

for $k=0,1,\cdots ,K$ .

Proof. From the first equation of (3.2) we have

$\begin{array}{l}{S}_j^{k+1}=\frac{S_j^k\left(1-\frac{\Delta t}{\Delta a}\right)+\frac{\Delta t}{\Delta a}{S}_{j-1}^k}{1+\Delta t\left(d_s+\alpha V^k\right)}\\ S_j^{k+1}\left(1+\Delta t\left(d_s+\alpha V^k\right)\right)=S_j^k\left(1-\frac{\Delta t}{\Delta a}\right)+\frac{\Delta t}{\Delta a}{S}_{j-1}^k\end{array}$

Thus,

$S_{j+1}^{k+1}\left(1+\Delta t\left(d_s+\alpha V^k\right)\right)=S_{j+1}^k\left(1-\frac{\Delta t}{\Delta a}\right)+\frac{\Delta t}{\Delta a}{S}_j^k$

Subtracting the above two equations and summing over $j=1,2,\cdots ,m-1$ :

$\begin{array}{l}{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left(S_{j+1}^{k+1}-S_j^{k+1}\right)\left(1+\Delta t\left(d_s+\alpha V^k\right)\right)\\ ={\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left(S_{j+1}^k-S_j^k\right)\left(1-\frac{\Delta t}{\Delta a}\right)+{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left(S_j^k-S_{j-1}^k\right)\frac{\Delta t}{\Delta a}\end{array}$

By the first equation of (3.2):

$S_1^{k+1}\left(1+\Delta t\left(d_s+\alpha V^k\right)\right)=S_1^k\left(1-\frac{\Delta t}{\Delta a}\right)+\frac{\Delta t}{\Delta a}{S}_0^k$

Combining:

$\begin{array}{l}TV\left(S^{k+1}\right)\left(1+\Delta t\left(d_s+\alpha V^k\right)\right)\\ =\left(\left|S_1^{k+1}-S_0^{k+1}\right|+{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left|S_{j+1}^{k+1}-S_j^{k+1}\right|\right)\left(1+\Delta t\left(d_s+\alpha V^k\right)\right)\end{array}$

Here,

$\begin{array}{l}\left|S_1^{k+1}-S_0^{k+1}\right|\left(1+\Delta t\left(d_s+\alpha V^k\right)\right)\\ =\left|S_1^k\left(1-\frac{\Delta t}{\Delta a}\right)+\frac{\Delta t}{\Delta a}{S}_0^k-S_0^{k+1}\left(1+\Delta t\left(d_s+\alpha V^k\right)\right)\right|\\ =\left|\frac{\Delta t}{\Delta a}\left(S_0^k-S_1^k\right)+S_1^k-S_0^{k+1}-S_0^{k+1}\Delta t\left(d_s+\alpha V^k\right)\right|\\ =\left|\left(1-\frac{\Delta t}{\Delta a}\right)\left(S_1^k-S_0^k\right)-\left(S_0^{k+1}-S_0^k\right)-S_0^{k+1}\Delta t\left(d_s+\alpha V^k\right)\right|\end{array}$

From the boundary condition,

$\begin{array}{l}\left|S_0^{k+1}-S_0^k\right|\\ ={\beta }_0|\left(S\left(a_{j_r^{+}},t_{k+1}\right)\left(a_{j_r^{+}}-a_r\right)+{\displaystyle \underset{j=j_r^{+}}{\overset{j_q^{-}}{\sum }}}{S}_j^{k+1}\Delta a+S\left(a_{j_q^{-}},t_{k+1}\right)\left(a_q-a_{j_q^{-}}\right)\right)\\ -\left(S\left(a_{j_r^{+}},t_k\right)\left(a_{j_r^{+}}-a_r\right)+{\displaystyle \underset{j=j_r^{+}}{\overset{j_q^{-}}{\sum }}}{S}_j^k\Delta a+S\left(a_{j_q^{-}},t_k\right)\left(a_q-a_{j_q^{-}}\right)\right)|\\ ={\beta }_0|\left(S\left(a_{j_r^{+}},t_{k+1}\right)-S\left(a_{j_r^{+}},t_k\right)\right)\left(a_{j_r^{+}}-a_r\right)+{\displaystyle \underset{j=j_r^{+}}{\overset{j_q^{-}}{\sum }}}\left(S_j^{k+1}-S_j^k\right)\Delta a\\ \begin{array}{c}\stackrel{}{}\\ \end{array}+\left(S(a_{j_q^{-}},t_{k+1}-S(a_{j_q^{-}},t_k\right)\left(a_q-a_{j_q^{-}}\right)|\end{array}$

$\begin{array}{l}\le {\beta }_0(\left|S\left(a_{j_r^{+}},t_{k+1}\right)-S\left(a_{j_r^{+}},t_k\right)\right|\left(a_{j_r^{+}}-a_r\right)+{\displaystyle \underset{j=j_r^{+}}{\overset{j_q^{-}}{\sum }}}\left|S_j^{k+1}-S_j^k\right|\Delta a\\ +\left|S\left(a_{j_q^{-}},t_{k+1}\right)-S\left(a_{j_q^{-}},t_k\right)\right|\left(a_q-a_{j_q^{-}}\right))\\ \le {\beta }_0{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|S_j^{k+1}-S_j^k\right|\Delta a\end{array}$

From the first equation in (3.1),

$\left(S_j^{k+1}-S_j^k\right)\Delta a=-\left(S_j^k-S_{j-1}^k\right)\Delta t-\left(d_s+\alpha V^k\right)S_j^{k+1}\Delta t\Delta a$

Thus,

$\begin{array}{c}\left|S_0^{k+1}-S_0^k\right|\le {\beta }_0{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|S_j^{k+1}-S_j^k\right|\Delta a\\ \le {\beta }_0{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|S_j^k-S_{j-1}^k\right|\Delta t+{\beta }_0{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left(d_s+\alpha V^k\right)S_j^{k+1}\Delta t\Delta a\\ \le {\beta }_{0}TV\left(S^k\right)\Delta t+{\beta }_0\left(d_s+\alpha M_2\right){\Vert S^{k+1}\Vert }_1\Delta t\end{array}$

Therefore,

$\begin{array}{c}TV\left(S^{k+1}\right)\le TV\left(S^{k+1}\right)\left(1+\Delta t\left(d_s+\alpha V^k\right)\right)\\ =\left(\left|S_1^{k+1}-S_0^{k+1}\right|+{\displaystyle \underset{j\mathrm{<mo>\; =\; </mo>1}}{\overset{m-1}{\sum }}}\left|S_{j+1}^{k+1}-S_j^{k+1}\right|\right)\left(1+\Delta t\left(d_s+\alpha V^k\right)\right)\\ \le \left(1-\frac{\Delta t}{\Delta a}\right)TV\left(S^k\right)+S_0^{k+1}\Delta t\left(d_s+\alpha V^k\right)\\ +\left|S_0^{k+1}-S_0^k\right|+TV\left(S^k\right)\frac{\Delta t}{\Delta a}\end{array}$

$\begin{array}{c}\le TV\left(S^k\right)+{\beta }_0{\Vert S^{k+1}\Vert }_1\Delta t\left(d_s+\alpha M_2\right)+{\beta }_{0}TV\left(S^k\right)\Delta t\\ +{\beta }_0\left(d_s+\alpha M_2\right){\Vert S^{k+1}\Vert }_1\Delta t\\ =\left(1+{\beta }_0\Delta t\right)TV\left(S^k\right)+({\beta }_0{\Vert S^{k+1}\Vert }_1\left(d_s+\alpha M_2\right)\\ +{\beta }_0\left(d_s+\alpha M_2\right){\Vert S^{k+1}\Vert }_1)\Delta t\end{array}$

Similarly,

$\begin{array}{l}{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left(I_{j+1}^{k+1}-I_j^{k+1}\right)\left(1+\Delta t\left(d_i+\delta \right)\right)\\ ={\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left(I_{j+1}^k-I_j^k\right)\left(1-\frac{\Delta t}{\Delta a}\right)+\frac{\Delta t}{\Delta a}{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left(I_j^k-I_{j-1}^k\right)+\Delta t\alpha V^k{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left(S_{j+1}^{k+1}-S_j^{k+1}\right)\end{array}$

$\begin{array}{l}\left|I_1^{k+1}-I_0^{k+1}\right|\left(1+\Delta t\left(d_i+\delta \right)\right)\\ \le \left(1-\frac{\Delta t}{\Delta a}\right)\left|I_1^k-I_0^k\right|+\left|I_0^{k+1}-I_0^k\right|+\Delta t\alpha V^k{S}_1^{k+1}\\ \le \left(1-\frac{\Delta t}{\Delta a}\right)\left|I_1^k-I_0^k\right|+{\beta }_0{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|I_j^{k+1}-I_j^k\right|\Delta a+\Delta t\alpha V^k{S}_1^{k+1}\end{array}$

$\begin{array}{l}\le \left(1-\frac{\Delta t}{\Delta a}\right)\left|I_1^k-I_0^k\right|+{\beta }_0({\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|I_j^k-I_{j-1}^k\right|\Delta t+{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left(d_i+\delta \right)I_j^{k+1}\Delta a\Delta t\\ +{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\alpha V^k{S}_j^{k+1}\Delta a\Delta t)+\Delta t\alpha V^k{S}_1^{k+1}\\ \le \left(1-\frac{\Delta t}{\Delta a}\right)\left|I_1^k-I_0^k\right|+{\beta }_0\left(TV\left(I^k\right)\Delta t+\left(d_i+\delta \right){\Vert I^{k+1}\Vert }_1\Delta t+\alpha V^k{\Vert S^{k+1}\Vert }_1\Delta t\right)\\ +\Delta t\alpha V^k{S}_1^{k+1}\end{array}$

Therefore,

$\begin{array}{c}TV\left(I^{k+1}\right)\le TV\left(I^{k+1}\right)\left(1+\Delta t\left(d_i+\delta \right)\right)\\ =\left(\left|I_1^{k+1}-I_0^{k+1}\right|+{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left|I_{j+1}^{k+1}-I_j^{k+1}\right|\right)\left(1+\Delta t\left(d_i+\delta \right)\right)\\ \le \left(1+{\beta }_0\Delta t\right)TV\left(I^k\right)+\Delta t(\alpha V^{k}TV\left(S^{k+1}\right)+{\beta }_0\left(d_i+\delta \right){\Vert I^{k+1}\Vert }_1\\ +\alpha V^k{\Vert S^{k+1}\Vert }_1+\alpha V^k{\Vert S^{k+1}\Vert }_{\infty })\end{array}$

Similarly, by the third equation of (3.1) and (3.3),

$\begin{array}{l}{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left|P_{j+1}^{k+1}-P_j^{k+1}\right|\left(1+\Delta t\left(d_p+{\mu }_p^k\right)\right)\\ \le \left(1-\frac{\Delta t}{\Delta a}\right){\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left|P_{j+1}^k-P_j^k\right|+\frac{\Delta t}{\Delta a}{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left|P_j^k-P_{j-1}^k\right|+\Delta t\delta {\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left|I_{j+1}^{k+1}-I_j^{k+1}\right|\end{array}$

and

$\begin{array}{l}\left|P_1^{k+1}-P_0^{k+1}\right|\left(1+\Delta t\left(d_p+{\mu }_p^k\right)\right)\\ =\left|\left(1-\frac{\Delta t}{\Delta a}\right)\left(P_1^k-P_0^k\right)-\left(P_0^{k+1}-P_0^k\right)-P_0^{k+1}\Delta t\left(d_p+{\mu }_p^k\right)\right|\\ \le \left(1-\frac{\Delta t}{\Delta a}\right)\left|P_1^k-P_0^k\right|+{\beta }_0{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|P_j^{k+1}-P_j^k\right|\Delta a+P_0^{k+1}\Delta t\left(d_p+{\mu }_p^k\right)\\ \le \left(1-\frac{\Delta t}{\Delta a}\right)\left|P_1^k-P_0^k\right|+{\beta }_0\Delta t\left(TV\left(P^k\right)+\left(d_p+{\mu }_p^k\right){\Vert P^{k+1}\Vert }_1+\delta {\Vert I^{k+1}\Vert }_1\right)\\ +P_0^{k+1}\Delta t\left(d_p+{\mu }_p^k\right)\end{array}$

Therefore,

$\begin{array}{c}TV\left(P^{k+1}\right)\le TV\left(P^{k+1}\right)\left(1+\Delta t\left(d_p+{\mu }_p^k\right)\right)\\ ={\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left|P_{j+1}^k-P_j^k\right|\left(1+\Delta t\left(d_p+{\mu }_p^k\right)\right)+\left|P_1^{k+1}-P_0^{k+1}\right|\left(1+\Delta t\left(d_p+{\mu }_p^k\right)\right)\\ \le \left(1+{\beta }_0\Delta t\right)TV\left(P^k\right)+\Delta t(\delta TV\left(I^{k+1}\right)+{\Vert P^{k+1}\Vert }_{\infty }\left(d_p+{\mu }_p^k\right)\\ +{\beta }_0\left(d_p+{\mu }_p^k\right){\Vert P^{k+1}\Vert }_1+{\beta }_0\delta {\Vert I^{k+1}\Vert }_1)\end{array}$

Similarly, by the fourth equation of (3.1) and (3.3),

$\begin{array}{l}{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left|C_{j+1}^{k+1}-C_j^{k+1}\right|\left(1+\Delta t{d}_c\right)\\ \le \left(1-\frac{\Delta t}{\Delta a}\right){\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left|C_{j+1}^k-C_j^k\right|+\frac{\Delta t}{\Delta a}{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left|C_j^k-C_{j-1}^k\right|+\Delta t{\mu }_p^k{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left|P_{j+1}^{k+1}-P_j^{k+1}\right|\end{array}$

and

$\begin{array}{l}\left|C_1^{k+1}-C_0^{k+1}\right|\left(1+\Delta t{d}_c\right)\\ =\left|\left(1-\frac{\Delta t}{\Delta a}\right)\left(C_1^k-C_0^k\right)-\left(C_0^{k+1}-C_0^k\right)-C_0^{k+1}\Delta t{d}_c+\Delta t{\mu }_p^k{P}_1^{k+1}\right|\\ \le \left(1-\frac{\Delta t}{\Delta a}\right)\left|C_1^k-C_0^k\right|+\left|C_0^{k+1}-C_0^k\right|+\Delta t\left(d_c{\Vert C^{k+1}\Vert }_{\infty }+{\mu }_p^k{\Vert P^{k+1}\Vert }_{\infty }\right)\\ \le \left(1-\frac{\Delta t}{\Delta a}\right)\left|C_1^k-C_0^k\right|+{\beta }_c{\displaystyle \underset{i=1}{\overset{m}{\sum }}}\left|C_j^{k+1}-C_j^k\right|\Delta a+\Delta t\left(d_c{\Vert C^{k+1}\Vert }_{\infty }+{\mu }_p^k{\Vert P^{k+1}\Vert }_{\infty }\right)\\ \le \left(1-\frac{\Delta t}{\Delta a}\right)\left|C_1^k-C_0^k\right|+{\beta }_c\Delta t\left(TV\left(C^k\right)+d_c{\Vert C^{k+1}\Vert }_1+{\mu }_p^k{\Vert P^{k+1}\Vert }_1\right)\\ +\Delta t\left(d_c{\Vert C^{k+1}\Vert }_{\infty }+{\mu }_p^k{\Vert P^{k+1}\Vert }_{\infty }\right)\end{array}$

Therefore,

$\begin{array}{c}TV\left(C^{k+1}\right)\le TV\left(C^{k+1}\right)\left(1+\Delta t{d}_c\right)\\ ={\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left|C_{j+1}^k-C_j^k\right|\left(1+\Delta t{d}_c\right)+\left|C_1^{k+1}-C_0^{k+1}\right|\left(1+\Delta t{d}_c\right)\\ \le \left(1+{\beta }_c\Delta t\right)TV\left(C^k\right)+\Delta t({\mu }_p^{k}TV\left(P^{k+1}\right)+{\beta }_c(TV\left(C^k\right)\\ +d_c{\Vert C^{k+1}\Vert }_1+{\mu }_p^k{\Vert P^{k+1}\Vert }_1)+d_c{\Vert C^{k+1}\Vert }_{\infty }+{\mu }_p^k{\Vert P^{k+1}\Vert }_{\infty })\end{array}$

The next lemma shows the numerical approximations satisfy an ${\mathcal{l}}^1$ Lipschitz-type condition in t.

Lemma 4.4 There exists $M_4>0$ , independent of $\Delta a$ and $\Delta t$ , such that for any integers $N_2>N_1\ge 0$ , the following estimates hold:

$\begin{array}{l}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|\frac{S_j^{N_2}-S_j^{N_1}}{\Delta t}\right|\Delta a\le M_4\left(N_2-N_1\right),\\ {\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|\frac{I_j^{N_2}-I_j^{N_1}}{\Delta t}\right|\Delta a\le M_4\left(N_2-N_1\right),\\ {\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|\frac{P_j^{N_2}-P_j^{N_1}}{\Delta t}\right|\Delta a\le M_4\left(N_2-N_1\right),\\ {\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|\frac{C_j^{N_2}-C_j^{N_1}}{\Delta t}\right|\Delta a\le M_4\left(N_2-N_1\right),\end{array}$

Proof. From the first equation of (3.1),

${\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|\frac{S_j^{k+1}-S_j^k}{\Delta t}\right|\Delta a\le TV\left(S^k\right)+\left(d_s+\alpha V^k\right){\Vert S^{k+1}\Vert }_1\le {\bar{M}}_1$

for some positive constant ${\bar{M}}_1$ . Hence,

${\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|\frac{S_j^{N_2}-S_j^{N_1}}{\Delta t}\right|\Delta a\le {\displaystyle \underset{k=N_1}{\overset{N_2-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|\frac{S_j^{k+1}-S_j^k}{\Delta t}\right|\Delta a\le {\bar{M}}_1\left(N_2-N_1\right)$

Similarly, there exist positive integers ${\bar{M}}_2\mathrm{,}{\bar{M}}_3\mathrm{,}{\bar{M}}_4$ such that

$\begin{array}{c}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|\frac{I_j^{N_2}-I_j^{N_1}}{\Delta t}\right|\Delta a\le {\displaystyle \underset{k=N_1}{\overset{N_2-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|\frac{I_j^{k+1}-I_j^k}{\Delta t}\right|\Delta a\\ \le \left(TV\left(I^k\right)+\left(d_i+\delta \right){\Vert I^{k+1}\Vert }_1+\alpha V^k{\Vert S^{k+1}\Vert }_1\right)\left(N_2-N_1\right)\\ \le {\bar{M}}_2\left(N_2-N_1\right)\mathrm{,}\end{array}$

$\begin{array}{c}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|\frac{P_j^{N_2}-P_j^{N_1}}{\Delta t}\right|\Delta a\le {\displaystyle \underset{k=N_1}{\overset{N_2-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|\frac{P_j^{k+1}-P_j^k}{\Delta t}\right|\Delta a\\ \le \left(TV\left(P^k\right)+\left(d_p+{\mu }_p^k\right){\Vert P^{k+1}\Vert }_1+\delta {\Vert I^{k+1}\Vert }_1\right)\left(N_2-N_1\right)\\ \le {\bar{M}}_3\left(N_2-N_1\right),\end{array}$

and

$\begin{array}{c}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|\frac{C_j^{N_2}-C_j^{N_1}}{\Delta t}\right|\Delta a\le {\displaystyle \underset{k=N_1}{\overset{N_2-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left|\frac{C_j^{k+1}-C_j^k}{\Delta t}\right|\Delta a\\ \le \left(TV\left(C^k\right)+d_c{\Vert C^{k+1}\Vert }_1+{\mu }_p^k{\Vert P^{k+1}\Vert }_1\right)\left(N_2-N_1\right)\\ \le {\bar{M}}_4\left(N_2-N_1\right)\mathrm{.}\end{array}$

Set $M_4={\max }_{i=1,2,3,4}\left\{{\bar{M}}_i\right\}$ and the results follow.

5. Convergence of the Difference Scheme and Existence of a Unique Weak Solution

Following similar notation as in [18] [19] , we define a family of functions $\left\{{\mathbb{U}}_1{}_{\left(\Delta a\mathrm{,}\Delta t\right)}\right\}$ , $\left\{{\mathbb{U}}_2{}_{\left(\Delta a\mathrm{,}\Delta t\right)}\right\}$ , $\left\{{\mathbb{U}}_3{}_{\left(\Delta a\mathrm{,}\Delta t\right)}\right\}$ , and $\left\{{\mathbb{U}}_4{}_{\left(\Delta a\mathrm{,}\Delta t\right)}\right\}$ by

$\begin{array}{l}{\mathbb{U}}_1{}_{\left(\Delta a,\Delta t\right)}\left(a,t\right)=S_i^k,{\mathbb{U}}_2{}_{\left(\Delta a,\Delta t\right)}\left(a,t\right)=I_j^k,\\ {\mathbb{U}}_3{}_{\left(\Delta a,\Delta t\right)}\left(a,t\right)=P_j^k,{\mathbb{U}}_4{}_{\left(\Delta a,\Delta t\right)}\left(a,t\right)=C_j^k,\end{array}$

and ${\mathbb{V}}_{\Delta t}\left(t\right)=V^{k-1}+\frac{V^k-V^{k-1}}{\Delta t}\left(t-t_{k-1}\right)$ , for $a\in \left[a_{j-1},a_j\right)$ , $t\in \left[t_{k-1},t_k\right)$ , $j=1,\cdots ,m$ , $k=1,\cdots ,K$ . Then by Lemmas 4.1 - 4.4, the each of the functions $\left\{{\mathbb{U}}_1{}_{\left(\Delta a\mathrm{,}\Delta t\right)}\right\}$ , $\left\{{\mathbb{U}}_2{}_{\left(\Delta a\mathrm{,}\Delta t\right)}\right\}$ , $\left\{{\mathbb{U}}_3{}_{\left(\Delta a\mathrm{,}\Delta t\right)}\right\}$ , $\left\{{\mathbb{U}}_4{}_{\left(\Delta a\mathrm{,}\Delta t\right)}\right\}$ is compact in the topology of ${\mathcal{L}}^1\left(\left(\mathrm{0,}{a}_d\right)\times \left(\mathrm{0,}T\right)\right)$ . Also, $\left\{{\mathbb{V}}_{\Delta t}\right\}$ is compact in the topology of $C\left[\mathrm{0,}T\right]$ .

Theorem 5.1 There exists a subsequence of functions $\left\{{\mathbb{U}}_1{}_{\left(\Delta a_r,\Delta t_r\right)}\right\}\subset \left\{{\mathbb{U}}_1{}_{\left(\Delta a\mathrm{,}\Delta t\right)}\right\}$ , $\left\{{\mathbb{U}}_2{}_{\left(\Delta a_r,\Delta t_r\right)}\right\}\subset \left\{{\mathbb{U}}_2{}_{\left(\Delta a\mathrm{,}\Delta t\right)}\right\}$ , $\left\{{\mathbb{U}}_3{}_{\left(\Delta a_r,\Delta t_r\right)}\right\}\subset \left\{{\mathbb{U}}_3{}_{\left(\Delta a\mathrm{,}\Delta t\right)}\right\}$ , $\left\{{\mathbb{U}}_4{}_{\left(\Delta a_r,\Delta t_r\right)}\right\}\subset \left\{{\mathbb{U}}_4{}_{\left(\Delta a\mathrm{,}\Delta t\right)}\right\}$ , and $\left\{{\mathbb{V}}_{\Delta t_r}\right\}\subset \left\{{\mathbb{V}}_{\Delta t}\right\}$ , which converge to functions $S\left(a\mathrm{,}t\right)\mathrm{,}I\left(a\mathrm{,}t\right)\mathrm{,}P\left(a\mathrm{,}t\right)\mathrm{,}C\left(a\mathrm{,}t\right)\in BV\left(\left[\mathrm{0,}{a}_d\right]\times \left[\mathrm{0,}T\right]\right)$ , and $V\left(t\right)\in C\left[\mathrm{0,}T\right]$ , respectively, in the sense that for all $t>0$

${\displaystyle {\int }_0^{a_d}}\left|{\mathbb{U}}_1{}_{\left(\Delta a_r,\Delta t_r\right)}\left(a\mathrm{,}t\right)-S\left(a\mathrm{,}t\right)\right|\text{d}a\to \mathrm{0,}$

${\displaystyle {\int }_0^{a_d}}\left|{\mathbb{U}}_2{}_{\left(\Delta a_r,\Delta t_r\right)}\left(a\mathrm{,}t\right)-I\left(a\mathrm{,}t\right)\right|\text{d}a\to 0$

${\displaystyle {\int }_0^{a_d}}\left|{\mathbb{U}}_3{}_{\left(\Delta a_r,\Delta t_r\right)}\left(a\mathrm{,}t\right)-P\left(a\mathrm{,}t\right)\right|\text{d}a\to \mathrm{0,}$

${\displaystyle {\int }_0^{a_d}}\left|{\mathbb{U}}_4{}_{\left(\Delta a_r,\Delta t_r\right)}\left(a\mathrm{,}t\right)-C\left(a\mathrm{,}t\right)\right|\text{d}a\to \mathrm{0,}$

${\displaystyle {\int }_0^T}{\displaystyle {\int }_0^{a_d}}\left|{\mathbb{U}}_1{}_{\left(\Delta a_r,\Delta t_r\right)}\left(a\mathrm{,}t\right)-S\left(a\mathrm{,}t\right)\right|\text{d}a\text{d}t\to 0,$

${\displaystyle {\int }_0^T}{\displaystyle {\int }_0^{a_d}}\left|{\mathbb{U}}_2{}_{\left(\Delta a_r,\Delta t_r\right)}\left(a\mathrm{,}t\right)-I\left(a\mathrm{,}t\right)\right|\text{d}a\text{d}t\to 0,$

${\displaystyle {\int }_0^T}{\displaystyle {\int }_0^{a_d}}\left|{\mathbb{U}}_3{}_{\left(\Delta a_r,\Delta t_r\right)}\left(a\mathrm{,}t\right)-P\left(a\mathrm{,}t\right)\right|\text{d}a\text{d}t\to 0,$

${\displaystyle {\int }_0^T}{\displaystyle {\int }_0^{a_d}}\left|{\mathbb{U}}_4{}_{\left(\Delta a_r,\Delta t_r\right)}\left(a\mathrm{,}t\right)-C\left(a\mathrm{,}t\right)\right|\text{d}a\text{d}t\to 0,$

$\underset{t\in \left[\mathrm{0,}T\right]}{\max }\left|{\mathbb{V}}_{\Delta t_r}\left(t\right)-V\left(t\right)\right|\to \mathrm{0,}$

as $r\to \infty$ (i.e., $\Delta a_r\mathrm{,}\Delta t_r\to 0$ ). Furthermore, there exist constants $M_5$ (depending on ${\Vert S_0\Vert }_{BV\left[\mathrm{0,}{a}_d\right]}$ , ${\Vert I_0\Vert }_{BV\left[\mathrm{0,}{a}_d\right]}$ , ${\Vert P_0\Vert }_{BV\left[\mathrm{0,}{a}_d\right]}$ , and ${\Vert C_0\Vert }_{BV\left[\mathrm{0,}{a}_d\right]}$ ) such that the limit function satisfy ${\Vert S\Vert }_{BV\left(\left[\mathrm{0,}{a}_d\right]\times \left[\mathrm{0,}T\right]\right)}\le M_5$ , ${\Vert V\Vert }_{BV\left(\left[\mathrm{0,}{a}_d\right]\times \left[\mathrm{0,}T\right]\right)}\le M_5$ , ${\Vert P\Vert }_{BV\left(\left[\mathrm{0,}{a}_d\right]\times \left[\mathrm{0,}T\right]\right)}\le M_5$ , ${\Vert C\Vert }_{BV\left(\left[\mathrm{0,}{a}_d\right]\times \left[\mathrm{0,}T\right]\right)}\le M_5$ , and ${\Vert V\Vert }_{C\left[\mathrm{0,}T\right]}\le M_5$ .

Proof. The results for $S\left(a\mathrm{,}t\right)$ , $I\left(a\mathrm{,}t\right)$ , $P\left(a\mathrm{,}t\right)$ , and $C\left(a\mathrm{,}t\right)$ follow from the proof of Lemma 16.7 on P. 276 in [18] . The results for $V\left(t\right)$ follows from Theorem I.28 (Ascoli’s Theorem) in [20] .

We show in the next theorem that the set of limit functions $\left\{S\left(a\mathrm{,}t\right)\mathrm{,}I\left(a\mathrm{,}t\right)\mathrm{,}P\left(a\mathrm{,}t\right)\mathrm{,}C\left(a\mathrm{,}t\right)\mathrm{,}V\left(t\right)\right\}$ constructed by the finite difference scheme is a weak solution to problem (2.1).

Theorem 5.2 The set of limit functions $\left\{S\left(a\mathrm{,}t\right)\mathrm{,}I\left(a\mathrm{,}t\right)\mathrm{,}P\left(a\mathrm{,}t\right)\mathrm{,}C\left(a\mathrm{,}t\right)\mathrm{,}V\left(t\right)\right\}$ defined in Theorem 5.1 is a weak solution of problem (2.1) and satisfies

${\Vert S\left(\cdot \mathrm{,}t\right)\Vert }_1+{\Vert I\left(\cdot \mathrm{,}t\right)\Vert }_1+{\Vert P\left(\cdot \mathrm{,}t\right)\Vert }_1+{\Vert C\left(\cdot \mathrm{,}t\right)\Vert }_1+\left|V\left(t\right)\right|\le M_1\mathrm{,}$

and

$\begin{array}{l}{\Vert S\Vert }_{L^{\infty }\left(\left(\mathrm{0,}{a}_d\right)\times \left(\mathrm{0,}T\right)\right)}+{\Vert I\Vert }_{L^{\infty }\left(\left(\mathrm{0,}{a}_d\right)\times \left(\mathrm{0,}T\right)\right)}+{\Vert P\Vert }_{L^{\infty }\left(\left(\mathrm{0,}{a}_d\right)\times \left(\mathrm{0,}T\right)\right)}+{\Vert C\Vert }_{L^{\infty }\left(\left(\mathrm{0,}{a}_d\right)\times \left(\mathrm{0,}T\right)\right)}+{\Vert V\Vert }_{c\left[\mathrm{0,}T\right]}\\ \le M_1+M_2\mathrm{.}\end{array}$

Proof. Let ${\zeta }_i\in C^1\left(\left[\mathrm{0,}{a}_d\right]\times \left[\mathrm{0,}T\right]\right)$ for $i=1,2,3,4$ . Denote the finite difference approximations ${\zeta }_i\left({\nu }_j\mathrm{,}{t}_k\right)$ by ${\left({\zeta }_i\right)}_j^k$ . Multiplying the first equation of the finite difference scheme (3.1) by ${\left({\zeta }_1\right)}_j^{k+1}$ and rearranging some of the terms, we obtain

$\begin{array}{l}{S}_j^{k+1}{\left({\zeta }_1\right)}_j^{k+1}-S_j^k{\left({\zeta }_1\right)}_j^k\\ =S_j^k\left({\left({\zeta }_1\right)}_j^{k+1}-{\left({\zeta }_1\right)}_j^k\right)-\frac{\Delta t}{\Delta a}\left({\left({\zeta }_1\right)}_j^{k+1}{S}_j^k-{\left({\zeta }_1\right)}_{j-1}^{k+1}{S}_{j-1}^k\right)\\ +\frac{\Delta t}{\Delta a}{S}_{j-1}^k\left({\left({\zeta }_1\right)}_j^{k+1}-{\left({\zeta }_1\right)}_{j-1}^{k+1}\right)-\left(d_s+\alpha V^k\right)S_j^{k+1}{\left({\zeta }_1\right)}_j^{k+1}\Delta t\end{array}$

Multiplying the above equation by $\Delta a$ , summing over $j=1,2,\cdots ,m$ , $k=0,1,\cdots ,K-1$ , and using the boundary condition we have

$\begin{array}{l}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left(S_j^K{\left({\zeta }_1\right)}_j^K-S_j^0{\left({\zeta }_1\right)}_j^0\right)\Delta a\\ ={\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left(S_j^k\left({\left({\zeta }_1\right)}_j^{k+1}-{\left({\zeta }_1\right)}_j^k\right)\Delta a-\left(d_s+\alpha V^k\right)S_j^{k+1}{\left({\zeta }_1\right)}_j^{k+1}\Delta a\Delta t\right)-{\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{S}_m^k{\left({\zeta }_1\right)}_m^{k+1}\Delta t\\ +{\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{\left({\zeta }_1\right)}_0^{k+1}{\beta }_0\left(S\left(a_{j_r^{+}},t_k\right)\left(a_{j_r^{+}}-a_r\right)+{\displaystyle \underset{j=j_r^{+}}{\overset{j_q^{-}}{\sum }}}{S}_j^k\Delta a+S\left(a_{j_q^{-}},t_k\right)\left(a_q-a_{j_q^{-}}\right)\right)\Delta t\\ +{\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{S}_{j-1}^k\left({\left({\zeta }_1\right)}_j^{k+1}-{\left({\zeta }_1\right)}_{j-1}^{k+1}\right)\Delta t,\end{array}$ (5.1)

Similarly,

$\begin{array}{l}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left(I_j^K{\left({\zeta }_2\right)}_j^K-I_j^0{\left({\zeta }_2\right)}_j^0\right)\Delta a\\ ={\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left(I_j^k\left({\left({\zeta }_2\right)}_j^{k+1}-{\left({\zeta }_2\right)}_j^k\right)\Delta a+I_{j-1}^k\left({\left({\zeta }_2\right)}_j^{k+1}-{\left({\zeta }_2\right)}_{j-1}^{k+1}\right)\right)\Delta t-{\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{I}_m^k{\left({\zeta }_2\right)}_m^{k+1}\Delta t\\ +{\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{\left({\zeta }_2\right)}_0^{k+1}{\beta }_0\left(I\left(a_{j_r^{+}},t_k\right)\left(a_{j_r^{+}}-a_r\right)+{\displaystyle \underset{j=j_r^{+}}{\overset{j_q^{-}}{\sum }}}{I}_j^k\Delta a+I\left(a_{j_q^{-}},t_k\right)\left(a_q-a_{j_q^{-}}\right)\right)\Delta t\\ +{\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left(-\left(d_i+\delta \right)I_j^{k+1}+\alpha V^k{S}_j^{k+1}\right){\left({\zeta }_2\right)}_j^{k+1}\Delta a\Delta t,\end{array}$ (5.2)

$\begin{array}{l}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left(P_j^K{\left({\zeta }_3\right)}_j^K-P_j^0{\left({\zeta }_3\right)}_j^0\right)\Delta a\\ ={\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left(P_j^k\left({\left({\zeta }_3\right)}_j^{k+1}-{\left({\zeta }_3\right)}_j^k\right)\Delta a+P_{j-1}^k\left({\left({\zeta }_3\right)}_j^{k+1}-{\left({\zeta }_3\right)}_{j-1}^{k+1}\right)\right)\Delta t-{\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{P}_m^k{\left({\zeta }_3\right)}_m^{k+1}\Delta t\\ +{\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{\left({\zeta }_3\right)}_0^{k+1}{\beta }_0\left(P\left(a_{j_r^{+}},t_k\right)\left(a_{j_r^{+}}-a_r\right)+{\displaystyle \underset{j=j_r^{+}}{\overset{j_q^{-}}{\sum }}}{P}_j^k\Delta a+P\left(a_{j_q^{-}},t_k\right)\left(a_q-a_{j_q^{-}}\right)\right)\Delta t\\ +{\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left(-d_p{P}_j^{k+1}+\delta I_j^{k+1}-{\mu }_P^k{P}_j^{k+1}\right){\left({\zeta }_3\right)}_j^{k+1}\Delta a\Delta t,\end{array}$ (5.3)

and

$\begin{array}{l}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left(C_j^K{\left({\zeta }_4\right)}_j^K-C_j^0{\left({\zeta }_4\right)}_j^0\right)\Delta a\\ ={\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left(C_j^k\left({\left({\zeta }_4\right)}_j^{k+1}-{\left({\zeta }_4\right)}_j^k\right)\Delta a+C_{j-1}^k\left({\left({\zeta }_4\right)}_j^{k+1}-{\left({\zeta }_4\right)}_{j-1}^{k+1}\right)\right)\Delta t-{\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{C}_m^k{\left({\zeta }_4\right)}_m^{k+1}\Delta t\\ +{\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{\left({\zeta }_4\right)}_0^{k+1}{\beta }_c\left(C\left(a_{j_r^{+}},t_k\right)\left(a_{j_r^{+}}-a_r\right)+{\displaystyle \underset{j=j_r^{+}}{\overset{j_q^{-}}{\sum }}}{C}_j^k\Delta a+C\left(a_{j_q^{-}},t_k\right)\left(a_q-a_{j_q^{-}}\right)\right)\Delta t\\ +{\displaystyle \underset{k=0}{\overset{K-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}\left(-d_c{C}_j^{k+1}+{\mu }_P^k{P}_j^{k+1}\right){\left({\zeta }_4\right)}_j^{k+1}\Delta a\Delta t.\end{array}$ (5.4)

Using the above fact and following a similar argument to that used in the proof of Lemma 16.9 on page 280 of [18] , it can be shown that the limit of the difference approximations in Theorem 5.1 is a weak solution to the problem (2.1) by letting $m\mathrm{,}K\to \infty$ . The bounds are obtained by taking the limit in the bounds of the difference approximations in Lemmas 4.1 and 4.2.

Remark 1 Uniqueness of the weak solution to the model (2.1) follows from similar arguments to those used in [13] . Thus, from Theorems 5.1 and 5.2 and the uniqueness of the weak solution, we have that the finite difference approximation (3.1) converges to the unique weak solution of system (2.1) in the sense given in Theorem 5.1.

6. Numerical Results

In order to check our code and confirm first order convergence, we chose a set of functions as an exact solution set and then added terms in order to ensure that they did indeed form a solution. To be more precise, we chose the solution set (6.1) and substituted it into (6.2) and then solved for the forcing functions $f_i$ , $i=1,2,\cdots ,5$ .

$\begin{array}{lll}S\left(a,t\right)={\text{e}}^{a+t},\hfill & I\left(a,t\right)={\text{e}}^{a-t},\hfill & P\left(a,t\right)={\text{e}}^{a-t^2}\hfill \\ C\left(a,t\right)={\text{e}}^{2a-t},\hfill & V\left(t\right)={\text{e}}^{-t}\hfill & \hfill \end{array}$ (6.1)

$\begin{array}{lll}\frac{\partial S\left(a,t\right)}{\partial t}+\frac{\partial S\left(a,t\right)}{\partial a}=-\left(d_s+\alpha V\left(t\right)\right)S\left(a,t\right)+f_1\left(a,t\right),\hfill & 0<t<T,\hfill & 0<a<a_d,\hfill \\ \frac{\partial I\left(a,t\right)}{\partial t}+\frac{\partial I\left(a,t\right)}{\partial a}=-\left(d_i+\delta \right)I\left(a,t\right)+\alpha V\left(t\right)S\left(a,t\right)+f_2\left(a,t\right),\hfill & 0<t<T,\hfill & 0<a<a_d,\hfill \\ \frac{\partial P\left(a,t\right)}{\partial t}+\frac{\partial P\left(a,t\right)}{\partial a}=-d_{p}P\left(a,t\right)+\delta I\left(a,t\right)-\mu \left(N_p\right)P\left(a,t\right)+f_3\left(a,t\right),\hfill & 0<t<T,\hfill & 0<a<a_d,\hfill \end{array}$

$\begin{array}{lll}\frac{\partial C\left(a,t\right)}{\partial t}+\frac{\partial C\left(a,t\right)}{\partial a}=-d_{c}C\left(a,t\right)+\mu \left(N_p\right)P\left(a,t\right)+f_4\left(a,t\right),\hfill & 0<t<T,\hfill & 0<a<a_d,\hfill \\ \frac{\text{d}V\left(t\right)}{\text{d}t}=\Lambda -d_v\left(V\right)V\left(t\right)+n{d}_i{\displaystyle {\int }_0^{a_d}}I\left(a,t\right)da+f_5\left(a,t\right),\hfill & 0<t<T,\hfill & \hfill \end{array}$ (6.2)

It is easily verified that this results in

$\begin{array}{l}{f}_1\left(a,t\right)={\text{e}}^{a+t}\left(2+d_s\right)+\alpha {\text{e}}^a\\ f_2\left(a,t\right)={\text{e}}^{a-t}\left(d_i+\delta \right)-\alpha {\text{e}}^a\\ f_3\left(a,t\right)={\text{e}}^{a-t^2}\left[1-2t+d_p+\mu \left(N_p\left(t\right)\right)\right]-\delta {\text{e}}^{a-t}\\ f_4\left(a,t\right)={\text{e}}^{2a-t}\left(1+d_c\right)-\mu \left(N_p\left(t\right)\right){\text{e}}^{a-t^2}\\ f_5\left(a,t\right)={\text{e}}^{-t}\left(d_v\left(V\left(t\right)\right)-1\right)-\Lambda -n{d}_i{N}_i\left(t\right),\end{array}$ (6.3)

where

$N_p\left(t\right)={\displaystyle {\int }_0^{a_d}}P\left(a,t\right)\text{d}a={\text{e}}^{a_d-t^2}-{\text{e}}^{-t^2}$

and

$N_i\left(t\right)={\displaystyle {\int }_0^{a_d}}I\left(a,t\right)\text{d}a={\text{e}}^{a_d-t}-{\text{e}}^{-t}.$

Upon inspecting the boundary conditions, we see that once $a_r\mathrm{,}{a}_q\mathrm{,}{a}_c$ , and $a_l$ are chosen, the parameters ${\beta }_0$ and ${\beta }_c$ must satisfy

${\beta }_0=\frac{1}{{\text{e}}^{a_q}-{\text{e}}^{a_r}}\text{and}{\beta }_c=\frac{2}{{\text{e}}^{2a_l}-{\text{e}}^{2a_c}}.$ (6.4)

The parameter values, except for ${\beta }_0$ and ${\beta }_c$ , were taken from [7] and are given in Table 1. The values of ${\beta }_0$ and ${\beta }_c$ that are implied by (6.4) are also provided there.

We ran five simulations, with the step sizes $\Delta t$ and $\Delta a$ being halved each time. Once simulations were done, we calculated the following error

${\Vert N_s-N_s^{\Delta t}\Vert }_{\infty }=\underset{0\le j\le m}{\max }\left|N_s\left(t_j\right)-N_s^{\Delta t}\left(j\right)\right|,$

where $N_s^{\Delta t}\left(j\right)$ is the numerical approximation of $N_s$ at $t_j$ with step size $\Delta t$ . This error was also calculated for the other components of the model, i.e., $N_i$ , $N_p$ , $N_c$ , and V. Finally, we calculated log₂ of the ratios of consecutive errors in order to determine the order of accuracy. The results in Table 2 and Table 3 demonstrate that the method is achieving first-order accuracy.

7. Simulations of Treatment

In this section we present numerical simulations of treatment with specific drugs. We follow the approach taken in [2] , but present more details. As we are interested in arresting the further development of unchecked levels of cancerous and precancerous cells, we chose the parameter values for case (iv) in [7] . Since the precancerous and cancerous cell populations grow without bound in this case, the model’s behavior corresponds to dysplasia and the onset of metastasis. The base parameters are given in Table 4.

Table 1. Parameter values for error simulations.

Table 2. Convergence of the first-order method.

Table 3. Convergence of the first-order method.

Table 4. Parameter values for error simulations.

The initial conditions used were $S_0\left(a\right)=2{\exp }^{-a}$ , $V_0=0$ , $I_0\left(a\right)=P_0\left(a\right)=C_0\left(a\right)=0$ .

One approach to treatment would be to increase the death rate ( $d_i$ ) of HPV-infected cells, thereby indirectly reducing the potential precancerous cell population. There are antiviral treatments such as Cidovofir that could facilitate this outcome, as noted in [2] [21] [22] . In order to numerically experiment with this scenario, we performed a simulation with the parameters given in Table 4, and then compared it to simulations in which this parameter is increase by 10%, 25%, and 50%. The results of these simulations are presented in Figure 1.

For this parameter set, the model is not particularly sensitive to the parameter $d_i$ . Indeed, an increase in the parameter $d_i$ by 50% does not reduce the precancerous and cancerous cell counts by 50%.

Another approach to treatment would be to increase the death rate of precancerous cells ( $d_p$ ). It was noted in [2] that the drugs 5-fluorouracil and bleomycin have shown success in this respect [23] . After experimenting with the parameter $d_p$ , we found that the model was quite sensitive to this parameter. Small percent changes cause large changes in model input. We settled on presenting the output resulting from the increasing $d_p$ by 3%, 5%, and 10%. The results are shown in Figure 2.

According to the model, this approach is vastly more efficient than the one of increasing the death rate of HPV-infected cells. Treatments of this type appear to have the ability to reduce unbounded precancerous and cancerous cell populations to stable ones that are more easily controlled.

There are, of course, other parameters in the model with which one can experiment. Here we have decided to focus on the two that can be affected by drugs which are currently available.

Figure 1. Model sensitivity to the parameter $d_i$ .

Figure 2. Model sensitivity to the parameter $d_p$ .

8. Conclusions

In this work, we have expanded on previous mathematical analyses of HPV-induced cervical cancer [2] [7] [8] [9] . We present a more general model in the sense that the transition rate $\mu$ of precancerous cells to cancerous cells can be any continuously differentiable function that satisfies $0<\mu \le 1$ . We also allow for any initial conditions that are of bounded variation on the domain $\mathrm{<mo\; stretchy="false">\; [\; </mo>0,}{a}_d\mathrm{<mo\; stretchy="false">\; ]\; </mo>}$ .

After presenting a first-order finite difference scheme for approximating solutions, we use said schemed to provide an existence-uniqueness result for the generalized model. Further, we provide numerical evidence that our numerical scheme is indeed of first-order.

Once the numerical scheme was implemented, we utilized it along with ideas from [2] [7] in order to examine treatment strategies that can be implemented with existing medications. More specifically, we examined the model’s sensitivity to increasing the death rates of HPV-infected cells ( $d_i$ ) and precancerous cells ( $d_p$ ). Of those two approaches, the model predicts that using drugs that increase the death rate of precancerous cells is the far more efficient treatment.

Further work on this subject may include applying our numerical scheme to the recently developed model in [9] .

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Burd, E.M. (2003) Human Papillomavirus and Cerivcal Cancer. Clinical Microbiology Reviews, 16, 1-17. https://doi.org/10.1128/CMR.16.1.1-17.2003
[2]	Asih, T.S.N., Lenhart, S., Wise, S., Aryati, L., Adi-Kusumo, F., Hardianti, M.S. and Forde, J. (2016) The Dynamics of HPV Infection and Cervical Cancer Cells. Bulletin of Mathematical Biology, 78, 4-20. https://doi.org/10.1007/s11538-015-0124-2
[3]	Barnabas, R.V., Laukkanen, P., Koskela, P., Kontula, O., Lehtinen, M. and Garnett, G.P. (2006) Epidemiology of HPV 16 and Cervical Cancer in Finlandand the Potential Impact of Vaccination: Mathematical Modelling and Analyses. PLOS Medicine, 3, e138. https://doi.org/10.1371/journal.pmed.0030138
[4]	Baussano, I., Ronco, G., Segnan, N., French, K., Vineis, P. and Garnett, G.P. (2010) HPV-16 Infection and Cervical Cancer: Modeling the Influence of Duration of Infection and Precancerous Lesions. Epidemics, 2, 21-28. https://doi.org/10.1016/j.epidem.2010.02.002
[5]	Brown, V.L. and Jane White, K.A. (2011) The Role of Optimal Control in Assessing the Most Cost-Effective Implementation of a Vaccination Programme: HPV as a Case Study. Mathematical Biosciences, 231, 126-134. https://doi.org/10.1016/j.mbs.2011.02.009
[6]	Lee, S.L. and Tameru, A.M. (2012) A Mathematical Model of Human Papillomavirus (HPV) in the United States and Its Impact on Cervical Cancer. Journal of Cancer, 3, 262-268. https://doi.org/10.7150/jca.4161
[7]	Akimeko, V.V. and Adi-Kusumo, F. (2021) Stability Analysis of an Age-Structured Model of Cervical Cancer Cells and HPV Dynamics. Mathematical Biosciences and Engineering, 18, 6155-6177. https://doi.org/10.3934/mbe.2021308
[8]	Aryati, L., Noor-Asih, T.S., Adi-Kusumo, F. and Hardianti, M.S. (2018) Global Stability of the Disease Free Equilibrium in a Cervical Cancer Model: A Chance to Recover. Far East Journal of Mathematical Sciences, 103, 1535-1546. https://doi.org/10.17654/MS103101535
[9]	Sari, E.R., Adi-Kusumo, F. and Aryati, L. (2022) Mathematical Analysis of a SIPC Age-Structured Model of Cervical Cancer. Mathematical Biosciences and Engineering, 19, 6013-6039. https://doi.org/10.3934/mbe.2022281
[10]	Woodman, C.B.J., Collins, S.I. and Young, L.S. (2007) The Natural History of Cervical HPV Infection: Unresolved Issues. Nature Reviews Cancer, 7, 11-22. https://doi.org/10.1038/nrc2050
[11]	Ackleh, A.S. and Deng, K. (2003) On a First-Order Hyperbolic Coagulation Model, Mathematical Methods in the Applied Sciences, 26, 703-715. https://doi.org/10.1002/mma.395
[12]	Ackleh, A.S., Deng, K. and Huang, Q. (2010) Existence-Uniqueness Results and Difference Approximations for an Amphibian Juvenile-Adult Model. Contemporary Mathematics, 513, 1-23.
[13]	Ackleh, A.S., Deng, K., Ito, K. and Thibodeaux, J.J. (2006) A Structured Erythropoiesis Model with Nonlinear Cell Maturation Velocity and Hormone Decay Rate. Mathematical Biosciences, 204, 21-48. https://doi.org/10.1016/j.mbs.2006.08.004
[14]	Banks, H.T., Cole, C.E., Schlosser, P.M. and Tran, H.T. (2004) Modeling and Optimal Regulation of Erythropoiesis Subject to Benzene Intoxication. Mathematical Biosciences & Engineering, 1, 15-48. https://doi.org/10.3934/mbe.2004.1.15
[15]	Edelstein, L. (1982) The Propagation of Fungal Colonies: A Model for Tissue Growth. Journal of Theoretical Biology, 98, 679-701. https://doi.org/10.1016/0022-5193(82)90146-1
[16]	Ackleh, A.S. and Ito, K. (1997) An Implicit Finite Difference Scheme for the Nonlinear Size-Structured Population Model. Numerical Functional Analysis and Optimization, 18, 865-884. https://doi.org/10.1080/01630569708816798
[17]	Angulo, O. and López-Marcos, J.C. (2004) Numerical Integration of Fully Nonlinear Size-Structured Population Models. Applied Numerical Mathematics, 50, 291-327. https://doi.org/10.1016/j.apnum.2004.01.007
[18]	Smoller, J. (1994) Shock Waves and Reaction-Diffusion Equations. Springer-Verlag, New York. https://doi.org/10.1007/978-1-4612-0873-0
[19]	Crandall, M.G. and Majda, A. (1980) Monotone Difference Approximations for Scalar Conservation Laws. Mathematics of Computation, 34, 1-21. https://doi.org/10.1090/S0025-5718-1980-0551288-3
[20]	Reed, M. and Simon, B. (1980) Methods of Modern Mathematical Physis I: Functional Analysis. Academic Press, New York.
[21]	Abdulkarim, B., Sabri, S., Deutsch, E., Chagraoui, H., Maggiorella, L., Thierry, J., Eschwege, F., Vainchenker, W., Chouaïb, S. and Bourhis, J. (2002) Antiviral Agent Cidofovir Restores P53 Function and Enhances the Radiosensitivity in HPV-Associated Cancers. Oncogene, 21, 2334-2346. https://doi.org/10.1038/sj.onc.1205006
[22]	Hostetler, K.Y., Rought, S., Aldern, K.A., Trahan, J., Beadle, J.R. and Corbeil, J. (2006) Enhanced Antiproliferative Effects of Alkoxyakyl Esters of Cidofovir in Human Cervical Cancer Cells in Vitro. Molecular Cancer Therapeutics, 5, 156-159. https://doi.org/10.1158/1535-7163.MCT-05-0200
[23]	Smoke, T. and Smoking, I. (2004) IARC Mongraphs on the Evaluation of Carcinogenic Risk to Humans. IARC Press, Lyon.