Eur. J. Epidemiol.0392-2990

Vot. 8, No. 4

EUROPEAN JOURNAL

July 1992,p. 585-593

OF EPIDEMIOLOGY

AN AIDS MODEL WITH DISTRIBUTED INCUBATION AND VARIABLE INFECTIOUSNESS: APPLICATIONS TO IV DRUG USERS IN LATIUM, ITALY M. IANNELL1.1, R. LORO*, F. MILNER**, A. PUGLIESE* and G. RABBIOLO** *Dipartimento di M a t e m a t i c a - Universit6 di Trento - 38050 - P o v o ( T N ) - Italy. **Dipartimento di M a t e m a t i c a - 11 Universitft di R o m a - Italy, a n d D e p a r t m e n t o f M a t h e m a t i c s P u r d u e University - W e s t L a f a y e t t e - I N 4 7 9 0 7 - U.S.A..

Key words: HIV/AIDS - Mathematical model - IV drug users An AIDS model with distributed incubation and variable infectiousness is considered and simulated via a second-order numerical method. The method is applied to the HIV epidemic among IV drug users in the Latium region of Italy, using available data on the lenght of the incubation period before the onset of AIDS, on the infectivity of infected individuals during that period, and on the demography of drug users. The contact rate is adjusted to match the actual number of AIDS cases. The sensitivity of the model to uncertainties in the parameters is finally investigated, by performing several simulations.

INTRODUCTION

Since HIV infection is characterized by a very long incubation period and apparently by variable infectiousness, the mathematical modelling of this disease has dealt with these two features from early on (1). As a matter of fact, models which rest upon the use of ordinary differential equations implicitly assume that the incubation period is exponentially distributed; experimental evidence shows that this is not the case. It seems, therefore, necessary to introduce a new variable 0, the time elapsed since the infection, and to formulate models in terms of partial differential equations. Moreover, this formulation also allows the modelling of variable infectiousness, which can be assumed to depend on 0. An approximation to this approach has been developed (13) by subdividing the incubation period into several stages, each exponentially distributed, so that the resulting model still rests on ordinary differential equations. 1 Corresponding author. 585

Here we are concerned with the spread of infection in a single homogenous population (as, for instance, homosexuals or IV drug users) which does not reproduce itself but has an external recruitment. This situation is described by the scheme in figure 1, where: S(t) = i(O,t) =

number of susceptible individuals at time t; density function (at time t) of individuals who were infected 0 time before;

A(t) = number of individuals in the AIDS class; A = recruitment rate of new susceptibles; tl

= exit rate, i.e. rate at which individuals die of causes other than AIDS, or leave the active compartment;

v V

= additional death rate due to AIDS; = removal rate, i.e. probability that an individual of infection age 0 develops AIDS in the time unit;

Iannelli M. et aL

Eur. J. Epidemiol.

K(t) = infection rate, i.e. rate at which susceptible individuals get infected;

Appendix 1 The model described in Section 1 translates into the following system of equations:

As for this latter variable, we assume the following constitutive form:

K(t)

=

Oi(O,t) . Oi(O,t) ........ t O0 = - ' r t v ) * W , ~ ) - ~i(o,

c(S(t) I(t)) fo °° ¢(0)i(0, t) dO S(t) ++ ±(t)

t) (i)

d A(t) = G(t) - (U + v)A(t)

with:

i(O, t) = K ( t ) S ( t )

0(0) = probability of contracting the infection (if susceptible) in a contact with an infected individual of infection age 0; c(T) = contact rate per individual, when total population size is T;

i(o, o) = io(O) >_ o .

I(t)

S

s ( o ) = So > o

with

= f o i(O, t) dO = total number of infectives;

_} io=o. , , q

£(OX)

c(t) =

I

no = c(hh,)

A

f0oo¢.O.e)

o~l(a)da'dO

(2)

in fact: i) if Ro < 1 the disease free equilibrium (S = A/la, I = A = 0) is globally stable: that is, the infection goes to extinction whatever the initial distribution and size of infectives. ii) ifRo > 1 the disease free state is unstable, i.e., the infection is sustained. Moreover, in this case an endemic state exists (S* I*), S * > O , I* > O. Some information is available about the endemic state; in particular, when the function c(T) is constant, is can be explicitly expressed as:

~+v 1 Figure 1. - Block scheme of the model.

The previous scheme translates into a system of differential equations which has been considered by several authors (1, 4, 6, 26). As in simpler epidemic models, a crucial role in the dynamics is played by the reproduction number Ro: if Ro < l the disease is not sustained while if Ro > 1 epidemics will develop and will generally approach an endemic level (see Appedix 1 for the exact formula of Ro and for a short account of the actual results). We performed several simulations of this system of equations, choosing parameter values obtained from available information on the demography and on the spread of HIV infection among IV drug users in Latium, Italy. Our simulations should give some idea about the way HIV infection has spread, possible outcomes and perhaps measures of control.

We have studied a numerical scheme to solve the system of equations in Appendix 1. The numerical

7(O)i(O,t)dO

Note that since the equations for S(t) and i(O,t) do not depend on the value of A(t), the model can be analyzed without the equation for A(t); the latter is kept only for comparison with available data. The qualitative behaviour of system (1) has been analyzed by Castillo-Chavez et aL (5) and by Thieme and Castillo-Chavez (16); here we give an account of the results obtained. A crucial role is played by the "reproduction number".

iI

METHODS

f0 °

F = A f~o e-Ua-f: -r(o)a,,d0[R° _ 11 " f 2 e-"O- f° "(~)a"dO + Ro - 1

s'=

(3)

a lg /~ +/',°° e-~a-J0 , , "t(°)~d0 ,

Ro - 1

Concerning its stability, the following is known: if o (0) is constant, the endemic equilibrium (when it exists) is always locally asymptotically stable (5); in this case it is generally believed that the stability is global, and our numerical evidence seems to confirm that. When 0(0) is not constant, several sufficient conditions for the (local) stability are given by Thieme and Castillo-Chavez (16); however it is also possible that the endemic equilibrium loses stability through a pair of roots of the characteristic equation passing through the imaginary axis (16), probably giving rise to a stable periodic solution. We note that the actual parameters of the model, as discussed in section 2, give a value Ro ~ 10; moreover, as it can actually be seen from equation (2), 1" tends to a limiting value when c (and consequently Ro) tends to 0% and this limit value is already approached for Ro ~ 10. This explains why varying c has a barely noticeable effect on the endemic level of the infection, as shown in figure 6.

586

VoL 8, 1992

An AIDS model and applications

If h is small enough, for each T > 0 there exists Q > such that the following estimate holds for all j and n such that nh < T:

scheme we used is a second order one, so that it allows fast results and good accuracy; some details on the method are given in Appendix 2.

I "1 + Io;'l + ICI _< Qh 2 That is, the method is second order accurate, as stated in Section 2.

APPENDIX 2

Here is a brief description of the numerical scheme. Its mathematical analysis will be performed elsewhere. Choosing the same step h for both variables t and o, we set: =-yOh),

=

Cj = COh)

and use as main variables S", t)", and An to approximate, respectively, S(nh), iOh, nh) and A(nh). Our scheme is then the following: Sn - S~-2 - A - ~ + Kn-l) Sn + S"-2 2h 2 A" - A n-2 v ) A" + A "-2 2h = - (1~ + 2 ~- G " - 1

.n

i;' + iy:

The exit rate la.

2h •

The exit rate la used in the model includes real deaths and other types of exit from the active population, such as old age or giving up the intravenous use of drugs. Here, for the sake of simplicity, we assume that individuals start using intravenous drugs when 15 years old and quit when 45, if they are still alive. Mortality rate among IV drug users is estimated to be about 10 times that of individuals in the same age-groups (19). Therefore we obtain, for IV drug users, a true yearly mortality rate m = 0.012. Adding to the mortality rate the proportion of individuals that are between 44 and 45 year old, and therefore will exit from the population within one year, we obtain

.n--1

h In

1.

oo

.~=1

n

j=l

[[1

Sn + In

-h

j=l

hi

la = m +

~o.,, = s n R n 1

(1 - m)2~m 1 - (1 - m) '°

- 0.04

T h e r e c r u i t m e n t rate A .

.~

oo

j=l

The first formulae in (5) are obtained approximating the derivatives at t = (n - 1)h with centered finite differences. The subsequent formulae approximate the integrals defining l(t), K(t) and G(t) according to the trapezoidal rule. Note that, although the summations are written forj up to 0% they are actually finite sums, if io(O) is positive only on a bounded interval. In fact, the interval on which in is positive increases by just one index at each time step; therefore, if io(o) = 0 for o > hNo, the sums are actually up to j = No + n. Formulae (5) are clearly only applicable for n _>2. For n = 1 a mixed implicit-explicit formula is used; this is not discussed here. As for the convergence of this method, setting

¢" = ~

The parameters used in the simulations of the model were chosen on the basis of information on the demography and the spread of HIV infection among intravenous drug users of the Latium region. To reflect uncertainties of parameters and/or conceivable measures of control, we also performed simulations with different parameter values. Information about the population of IV drug users, as well as the data on AIDS cases in that population, have been provided by the "Osservatorio Epidemiologico della Regione Lazio" (20). Some details on the choice of the parameters follow:

-

As for the recruitment rate, the estimate of the population of IV drug users we use is 30,000 people (20), although more recent estimates are lower than that (19). Since the model, in absence of disease, yields a stationary population of A/la (see Appendix I), we estimated A = 30,000 x la = 1,200. In order to study the sensitivity of the model, we performed one simulation with A = 600, corresponding to a population o f 15,000 IV drug users; this is probably a lower estimate for the true size o f the population, while 30,000 should be an upper estimate. Certainly better, direct estimates on the exit rate (not due to AIDS) from the active population, and on the recruitment rate, would be helpful.

s n

,n

= i ( j h , n h ) - ~j

~n = A ( n h ) - A n,

T h e i n c u b a t i o n curve y.

Discussions on the curve y are based on data on the incubation periods of selected cohorts of infected individuals. These data pose difficult statistical 587

Iannelli M. et aL

Eur. J. Epidemiol.

problems, discussed by several authors (7, 14, 15, 22, 25). We used the curve y (see Fig. 2) obtained by Bacchetti and Moss (3) analyzing, in a non-parametric way, data on seropositivity in selected cohorts and on AIDS cases in San Francisco.

Infectivity Curve

Phi(Thetal

0.1 9. ooooooe-o2 8. ooooooe-o2 7, ooooooe-o2 0.00000oe-02 5. DOOOOOe-02

Incubation

Curve

Gamma(Theta)

.0.14

0.12 •

4. 000000~-02

~~

3,000000e-02

2. O00000e-02

0.1

I. O 0 0 0 0 0 e - 0 2

8. oO0000e-02 0

1

2

3

4

5 Theta

6

7

8

9

i0

6, ooooooe-02

Figure 3. - T h e infectivity curve o (0).

4 . oo0oooe-02

2, ooooooe-o2

o

[

I

I

I

I

I

I

I

[

I

2

3

4

5 '~heta

6

7

O

9

Figure 2. - T h e i n c u b a t i o n curve y(0)

The infective curve o .

The curve o we used (Fig. 3) comes from Blythe and Anderson (4), using data on H1V antigen and antibody titers by Pedersen et al. (18). The curve reflects the common belief that HIV virus has two peaks of infectivity, one shortly after the infection and the second when the infected individual begins to progress towards AIDS. Obviously, actual data on the infectivity with varying infection age do not exist, and one has to be aware that alternative patterns of HIV infection have been found (12). However, there seems to be enough evidence for a two-peaked infectivity to justify a theoretical study of the epidemiological implications of such a fact. Several studies (16, 17, 21) report a very strong variability in infectivity per partnership; it seems likely that infection age plays a strong role in that pattern. Note that the absolute numbers chosen for o are not important, since o enters the model only multiplied times c(T) and c is estimated from data (see below). We make no claim that the function o of figure 3 really represents the infectivity per contact or per partnership, nor that c represents the rate of partner exchange. Only the product co will reflect reality. The contact number

correct cumulative number of AIDS cases up to December 1989, as available in the spring of 1990. Due to later notifications, the model does not fit exactly the number of AIDS cases reported now and used in figure 4. However, the initial value for the infective population is also unknown. Therefore, we made four estimates about the initial distribution, namely that there were 1 or 30 individuals just infected in 1981 or 1979. These are based on data on seropositivity among IV drug users in several Italian cities (10). Certainly we could perform a joint estimate of initial values and c, using some statistical assumptions. Given the uncertainty in all parameters, we preferred a simple method to obtain an estimate of the order of magnitude of c. Using this method, we obtained four estimated for c in Table 1. The simulation obtained, assuming a single infected individual in 1981, best reproduces the overall spread up to 1989, as can be seen in figure 4. There we show, in linear and logarithmic scales, the number of AIDS cases predicted by the model up to 1993, compared with data up to 1989. One may note that the first part of the curve is essentially exponential, while it bends down later.

TABLE 1. - Values of c that fit the cumulative number of AIDS cases with different initial values. Values for c

c(T).

C(T) has been chosen to be a constant, at least for Taround 30,000, the only part of the curve that may be empirically determined. This constant has been chosen trying to fit the data on AIDS cases. We used the simulator and chosen the value of c that gave the

1 infectives

30 infectives

in 1979

45.92

37.96

in 1981

50.82

41.75

588

Vol. 8, 1992

A n AIDS model and applications

RESULTS 800

-

-

Mo0et

•*"Q""

Dala

-

Mo~el

200

3

6

9

Yoar - 1981

100

-

• --..v- - *