Well-Posedness and a Finite Difference Approximation for a Mathematical Model of HPV-Induced Cervical Cancer ()
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 (
), infected cell (
), precancerous cell (
), and cancerous cell (
) populations, as well as the number of HPV viral particles (
), at time t and age a. The total populations of the susceptible, infected, precancerous, and cancerous cell types at time t are denoted by
,
,
, and
, respectively, where
is the maximum lifespan of the cells. We denote their respective death rates by
, and
. The infection of healthy cells by free virus particles is modeled by the mass action term
. Infected cells are partially moved into the precancerous cell population, modeled by the term
. The rate at which precancerous cells move into the cancerous cell population is modeled dynamically with the bounded function
, although the theoretical results we provide in this work are valid for any continuously differentiable function
with
. 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
. There is a density-dependent death rate
. Since the total population of infected cells dies at a rate of
, 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
, while that of cancerous cells is
. The reproductive age range of cancer cells is shifted relative to that of the noncancerous cells [7] , and so it is assumed that
and
. Noncancerous cells are assumed to have a birth rate of
, while cancerous cells have a birth rate of
. Finally, the initial cell densities are given by
,
,
, and
, and the initial viral load is
.
These definitions and assumptions lead to the following model:
(2.1)
An uncommon aspect of this model is that the integrals in the boundary conditions are not over the entire interval
. 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
is continuous.
2) The function
is continuous and satisfies
for all x in
.
3) The parameters
are nonnegative constants.
4) The initial conditions
are nonnegative with
Now we give the definition of a weak solution to problem (2.1) as follows: A 5-tuple
is called a weak solution to the problem (2.1) if it satisfies:
for each
, every test function
. Here,
,
,
, and
. First, we divide the intervals
and
into m and K subintervals, respectively, and use the following notation throughout the paper:
, and
denote the age and time mesh sizes, respectively. The mesh points are given by:
,
and
,
. We denote by
,
,
, and
, the finite difference approximations of
,
,
, and
, respectively. Further, we denote
by
,
by
, and
by
. We approximate the total cell populations with
where
. Since the parameters
and
will not generally coincide with our mesh points, we define the following:
The total variation,
norm, and
norm are defined as
respectively. The first-order scheme that will be considered here is given by:
(3.1)
The boundary conditions are computed as follows:
(3.2)
The initial conditions for p and m are computed as follows:
Equation (3.1) is more conveniently used in the following form:
(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
.
4. Estimates for the Finite Difference Approximations
The lemma below shows that the numerical approximations are bounded in the
norm.
Lemma 4.1 There exists a positive constant
, independent of
and
, such that
for
.
Proof. Multiply the first four equations of (3.1) on both sides by
and
and sum over
(4.1)
Multiply the fifth equation of (3.1) on both sides by
(4.2)
Adding up the five equations in (4.1) and (4.2) we have
From the boundary condition, we have
Therefore,
Let
. Then
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
, independent of
and
, such that
for
.
Proof. If
, then
If
is not obtained on the boundary, then there exists a
with
such that
.
From the first equation of (3.1)
That is,
Thus,
This implies that
is bounded since
is bounded. Similarly, if
, then
If
is not obtained on the boundary, then there exists a
with
such that
.
From the second equation of (3.1)
Therefore,
Similarly, if
, then
If
is not obtained on the boundary, then there exists a
with
such that
.
From the third equation of (3.1)
Therefore,
Similarly, if
, then
If
is not obtained on the boundary, then there exists a
with
such that
.
From the third equation of (3.1)
Therefore,
In the next lemma we show that the difference approximations are of bounded total variation.
Lemma 4.3 There exists a positive constant
, independent of
and
, such that
for
.
Proof. From the first equation of (3.2) we have
Thus,
Subtracting the above two equations and summing over
:
By the first equation of (3.2):
Combining:
Here,
From the boundary condition,
From the first equation in (3.1),
Thus,
Therefore,
Similarly,
Therefore,
Similarly, by the third equation of (3.1) and (3.3),
and
Therefore,
Similarly, by the fourth equation of (3.1) and (3.3),
and
Therefore,
The next lemma shows the numerical approximations satisfy an
Lipschitz-type condition in t.
Lemma 4.4 There exists
, independent of
and
, such that for any integers
, the following estimates hold:
Proof. From the first equation of (3.1),
for some positive constant
. Hence,
Similarly, there exist positive integers
such that
and
Set
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
,
,
, and
by
and
, for
,
,
,
. Then by Lemmas 4.1 - 4.4, the each of the functions
,
,
,
is compact in the topology of
. Also,
is compact in the topology of
.
Theorem 5.1 There exists a subsequence of functions
,
,
,
, and
, which converge to functions
, and
, respectively, in the sense that for all
as
(i.e.,
). Furthermore, there exist constants
(depending on
,
,
, and
) such that the limit function satisfy
,
,
,
, and
.
Proof. The results for
,
,
, and
follow from the proof of Lemma 16.7 on P. 276 in [18] . The results for
follows from Theorem I.28 (Ascoli’s Theorem) in [20] .
We show in the next theorem that the set of limit functions
constructed by the finite difference scheme is a weak solution to problem (2.1).
Theorem 5.2 The set of limit functions
defined in Theorem 5.1 is a weak solution of problem (2.1) and satisfies
and
Proof. Let
for
. Denote the finite difference approximations
by
. Multiplying the first equation of the finite difference scheme (3.1) by
and rearranging some of the terms, we obtain
Multiplying the above equation by
, summing over
,
, and using the boundary condition we have
(5.1)
Similarly,
(5.2)
(5.3)
and
(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
. 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
,
.
(6.1)
(6.2)
It is easily verified that this results in
(6.3)
where
and
Upon inspecting the boundary conditions, we see that once
, and
are chosen, the parameters
and
must satisfy
(6.4)
The parameter values, except for
and
, were taken from [7] and are given in Table 1. The values of
and
that are implied by (6.4) are also provided there.
We ran five simulations, with the step sizes
and
being halved each time. Once simulations were done, we calculated the following error
where
is the numerical approximation of
at
with step size
. This error was also calculated for the other components of the model, i.e.,
,
,
, and V. Finally, we calculated log2 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
,
,
.
One approach to treatment would be to increase the death rate (
) 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
. Indeed, an increase in the parameter
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 (
). It was noted in [2] that the drugs 5-fluorouracil and bleomycin have shown success in this respect [23] . After experimenting with the parameter
, 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
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
.
Figure 2. Model sensitivity to the parameter
.
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
of precancerous cells to cancerous cells can be any continuously differentiable function that satisfies
. We also allow for any initial conditions that are of bounded variation on the domain
.
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 (
) and precancerous cells (
). 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] .