Numerical and bifurcation analyses for a population model of HIV chemotherapy

Abstract

A competitive implicit finite-difference method will be developed and used for the solution of a non-linear mathematical model associated with the administration of highly-active chemotherapy to an HIV-infected population aimed at delaying progression to disease. The model, which assumes a non-constant transmission probability, exhibits two steady states; a trivial steady state (HIV-infection-free population) and a non-trivial steady state (population with HIV infection). Detailed stability and bifurcation analyses will reveal that whilst the trivial steady state only undergoes a static bifurcation (single zero singularity), the non-trivial steady state can not only exhibit static and dynamic (Hopf) bifurcations, but also a combination of two types of bifurcation (a double zero singularity).

Although the Gauss–Seidel-type method to be developed in this paper is implicit by construction, it enables the various sub-populations of the model to be monitored explicitly as time t tends to infinity. Furthermore, the method will be seen to be more competitive (in terms of numerical stability) than some well-known methods in the literature. The method is used to determine the impact of the chemotherapy treatment by comparing the population sizes at equilibrium of the treated and untreated infecteds.

Introduction

The Human Immunodefiency Virus (HIV) is the causative agent of Acquired Immunodefiency Syndrome in humans (AIDS). HIV is largely transmitted sexually, by IV drug use or through mother-to-child transmission. Over the past decade, substantial progress has been made in terms of HIV drugs development (see, for instance, [7]). Thanks to such advances, HIV replication and most opportunistic infections can be controlled without eliminating the virus using highly-active anti-retroviral drugs (HAART) [5].

Because of the ever-increasing numbers of reported cases of HIV infection and AIDS, much collaborative research is being conducted by mathematicians, biologists and physicians with the hope of gaining better insight into the transmission dynamics of HIV and designing effective methods of control (see [3], [4], [5], [6], [7], [8], [9], [10], [12], [14]).

Anderson et al. [3], Kiessling et al. [8] and Knox [10] developed epidemiological models which break down the total population into sub-groups according to their sexual preferences and rate of partner change. Dietz et al. [4] generalized the classical Susceptible-Infective-Susceptible (SIS) models (that is, models without immunity) by introducing heterogeneous contact and pairing rates. Ho et al. [7], Kirschner and Webb [9] and Wei et al. [14] studied the dynamics of HIV infection in patients receiving drug treatment. McLean et al. [12] provided a mathematical explanation of the co-existence between drug-resistant and drug-sensitive virus in treated HIV patients.

Most of these studies, however, appear to give scant or no detail of the numerical methods used, together with their associated stability analyses, to solve the resulting initial value problems (IVP’s).

This paper focusses on the bifurcation and stability of the equilibrium states, and the development of a Gauss–Seidel-like, implicit, finite-difference scheme (see [11], [13]), for a modified homogeneous model for the administration of anti-viral therapy developed by Gupta et al. [6] with a non-constant probability of transmission. The model was obtained by sub-dividing the infected and uninfected populations into treated and untreated categories. The type of treatment used was the administration of chemotherapy to infectious persons (which acts to delay progression to disease). Chemotherapy here refers to any drug treatment that is effective in the control of HIV replication, most importantly, HAART-combination therapy. The model enables the populations of susceptibles, together with treated and untreated infecteds, to be monitored separately as t→∞.

It is shown that the model exhibits two steady states, one of which can only have a static bifurcation, jumping to the other equilibrium solution and, consequently, the two states exchange their stability properties on the same static stability boundary. The second equilibrium solution, on the other hand, can have both a single zero singularity and a Hopf singularity, and may also have a co-dimension two singularity (that is, a double zero bifurcation).

It is worth mentioning that discretizing the ordinary differential equations (ODE’s) of the model by explicit Euler finite-difference schemes, can result in contrived chaos whenever the discretization parameters exceed certain values [11], [13]. The same applies to the use of explicit Runge–Kutta methods. Although chaos can often be avoided, even for Euler methods, by using small time steps, the extra computing costs incurred when examining the long-term behaviour of a dynamical system may be substantial. It is, therefore, essential to use a numerical method which allows the largest possible time steps that are consistent with stability and accuracy.

To avoid contrived chaos, whilst retaining accuracy and numerical stability, it may be necessary to forego the ease-of-implementation of inexpensive, explicit, numerical methods in favour of implicit methods (which are, at times, computationally-intensive). The novelty of the implicit finite-difference scheme to be developed in this paper is that it is easy to implement, fast-convergent (expected of Gauss–Seidel-type methods) and has very competitive stability properties.