Measures of concurrency in networks and the spread of infectious disease

ELSEVIER

Measures of Concurrency in Networks and the Spread of Infectious Disease MIRJAM KRETZSCHMAR National Institute of Public Health and Environmental Protection (RIVM), 3720 BA BilthoL,en, The Netherlands AND

MARTINA MORRIS Department of Sociology and School of Public Health, Columbia University, New York, New York 10025

Received 5 December 1994; revised 10 July 1995

ABSTRACT An investigation is made into the impact of concurrent partnerships on epidemic spread. Starting from a definition of concurrency on the level of individuals, the authors define ways to quantify concurrency on the population level. An index of concurrency based on graph theoretical considerations is introduced, and the way in which it is related to the degree distribution of the contact graph is demonstrated. Then the spread of an infectious disease on a dynamic partnership network is investigated. The model is based on a stochastic process of pair formation and separation and a process of disease transmission within partnerships of susceptible and infected individuals. Using Monte Carlo simulation, the spread of the epidemic is compared for contact patterns ranging from serial monogamy to situations where individuals can have many partners simultaneously. It is found that for a fixed mean number of partners per individual the distribution of these partnerships over the population has a major influence on the speed of the epidemic in its initial phase and consequently in the number of individuals who are infected after a certain time period.

1.

INTRODUCTION

T h e spread of sexually t r a n s m i t t e d diseases (STDs), especially A I D S , d e p e n d s o n the p a t t e r n s of sexual contacts within a p o p u l a t i o n . H e r e the cultural a n d m o r a l t r a d i t i o n s of the p o p u l a t i o n play a n i n f l u e n t i a l role: whereas in s o m e societies serial m o n o g a m y is the n o r m a n d having m o r e t h a n o n e p a r t n e r at the same time is an exception, in o t h e r societies polygamy is the n o r m or at least widely accepted. T h e contact

MATHEMATICAL BIOSCIENCES 133:165-195 (1996) 0025-5564/96/$15.00 © Elsevier Science Inc., 1996 SSDI 0025-5564(95)00093-S 655 Avenue of the Americas, New York, NY 10010

166

MIRJAM KRETZSCHMAR AND MARTINA MORRIS

patterns formed in different populations may differ in structure even if the number of partners per person as a population average is the same. The network of sexual contacts in a population is determined by three variables: the mean number of partners per person per unit time, the duration of partnerships, and the distribution of the partnerships over the population. These factors influence the way in which an STD spreads in the population. In this article we focus on the role of simultaneous or concurrent partnerships or, more generally, on the role of the distribution of partnerships within the population in the spread of STDs. In mathematical modeling of STDs the focus of attention has been on the modeling of heterogeneity in the population with respect to sexual activity and other attributes and the mixing between different subgroups [1-12]. An important ingredient in these models is the so-called mixing matrix, which describes the distribution of contacts within and among subgroups. Apart from random or proportional mixing, many different types of mixing matrices have been proposed. In [3, 9, 11, 13], ways of separating random and selective effects in the mixing matrix are discussed. With respect to the description of the contact process, models can be classified according to their partnership dynamics. At one end of the spectrum are mixing models as mentioned above in which duration of contacts is neglected, contacts are assumed to take place instantaneously, and every contact is with a new partner. In that case the notion of concurrent partnerships does not make sense. At the other end of the spectrum are models that assume a static structure of partnerships and disease transmission within these partnerships (e.g., lattice models [14, 15] and static network models). Between these two extremes are pair-formation models as introduced into epidemiological modeling by Dietz and Hadeler [16, 17]. In pair-formation models, duration of partnerships is taken into account by modeling a dynamic process of pair formation and separation. In most cases these models are based on the assumption that the population is monogamous in the sense that an individual cannot have more than one partnership simultaneously [18, 19]. A stochastic pair formation model has been discussed in [20]. It describes the process of pair formation and separation on the population level, that is, the numbers of pairs and singles are counted, but the partnerships themselves are not modeled explicitly. In contrast, in dynamic network models as introduced in [21, 22], partnerships are modeled on the individual level, resulting in a complete description of the contact network, and on the population level the dynamics are described by a stochastic pair-formation model.

C O N C U R R E N C Y AND EPIDEMICS

167

To model concurrent partnerships, Dietz and Tudor [23] extended a pair-formation model to include "triangles," by allowing an existing pair to form a relationship with a single. More generality, however, quickly leads to large systems of equations. In a paper by Watts and May [24] the overlapping of partnerships in time is modeled indirectly by incorporating a time delay between sexual contact and the infection of the next partner. A model in which partnerships form independently of each other for every individual (i.e., where the number of partners per individual is approximately Poisson distributed) is discussed by Altmann [25]. Furthermore, simulation models have been constructed [22, 26] that include partnership dynamics and multiple partnerships, but there the problem is the limited analytic tractability. In this article we study two aspects of concurrent partnerships. First we discuss how to quantify the "amount" of concurrency in a network; then we use a stochastic simulation model to study the spread of infection in different network structures. 2.

CONCURRENCY IN NETWORKS

By a sexual network we mean the total population together with all the partnerships at a given moment. In graph theoretical terms, individuals are referred to as vertices or nodes and a partnership as an edge between two vertices. The network is then the set of all vertices together with the set of all edges. So the sexual network may also contain isolated nodes, that is, individuals who are single. Concurrent partnerships are described by edges that emanate from the same vertex. In this section we do not consider dynamic aspects of a network. Rather, we assume that we have a population of size N with a given fixed partnership constellation. 2.1.

THE FRACTION OF SINGLES

In a completely monogamous population the number of persons involved in a partnership equals twice the number of partnerships. For every concurrent partnership the number of persons involved in partnerships decreases by 1. So for a given number of partnerships the fraction of the population who are single gives an indication of the number of concurrent partnerships. We denote by N P C m

the the the the

total population size, total number of partnerships, number of persons involved in a partnership, mean number of partners per person.

168

MIRJAM KRETZSCHMAR AND MARTINA MORRIS

So N - C is the number of singles. The mean number of partners per person is given as m = 2 P / N . Now we define K~ as the ratio of the number of persons involved in a partnership to twice the number of partnerships: C K1 := 2---fi"

(1)

For a monogamous population we have K1 1; for a population with concurrent partnerships, K1 < 1. We can write =

CN K1 = -2-ff-N

C raN"

(2)

With c .'= C / N the fraction of "non-singles" in the population, we get K 1 = c/m. In graph theoretical terms the number of partners of an individual is described by the degree of the vertex. The degree distribution in the network is a discrete distribution on t~ 0. We denote by r i the fraction of the population with i partners. If m is the mean of that distribution, we get 1 -- F 0

K~

m

'

(3)

where r 0 is the fraction of singles in the population. If the degree distribution can be described by a Poisson distribution with mean m, we get K1

l_e-m m

(4)

For a negative binomial with mean m and clumping p a r a m e t e r k we get

m[l

=1

k

k

(5)

These are monotone decreasing functions of the mean. The quantity K1 describes what fraction of the population is responsible for all existing partnerships. For a given mean, a small value of K1 means that these partnerships are concentrated on a small part of the population and therefore many partnerships are concurrent. However, the measure Ka cannot distinguish between various distributions of these partnerships within the active fraction c. Consider, for example,

CONCURRENCY AND EPIDEMICS

169

the graphs shown in Figures l b and c. T h e value of K1 is the same for both graphs, although their structure differs in a way that might be important for the spread o f a disease. T o define a measure that is sensitive to these differences we have to take m o r e details o f the network structure into account. 2.2. THE CONTACT GRAPH W e want to formalize the graph theoretical description o f the contact graph. T h e individuals of the population form a set o f vertices V = {v 1..... UN}. If a partnership exists between two vertices u i and vj, then the pair (vi, v j) is an element of the set of edges E. W e will use parentheses to indicate pairs, and we assume that (v i, vj)= (Vy, Vi) and

G H I

L(G) •

H H (a)

a

L(G)

ID

(b)

c

L(G)

H (c) FIG. 1. The graphs in (a), (b), and (c) have the same number of nodes and the same number of edges but differ in the distribution of the edges. The index of concurrency is K3 = 0 for (a), K3 = 1.5 for (b), and /£'3 0.75 for (c). =

170

MIRJAM KRETZSCHMAR AND MARTINA MORRIS

v i --/: vj, that is, (vi, vj) is a two-element subset of V.

More precisely, if the symbol ~ denotes "has a partnership with," then we define a function X v on V × V as follows: {~

)(v( Ui'Uj)

=

if v i v j, if v i ~ vj. ~

(6)

The function X v defines a matrix that is called the adjacency matrix in the literature [27, 28]. Let N = IVI denote the number of vertices in V. Then E is defined as

e = {(v,, O

v x v: x v ( v , , O

= 1}.

(7)

L e t e i denote the elements of E, and let P = IEI denote the number of

edges of E. Then E = {e 1..... e e l The graph G is given by G = (V, E). We can now describe the structure of the sexual contact network in terms of graph theory. The degree of a vertex v c V is defined as d e g ( v ) := E

(8)

Xv(Ui'U) "

Ui~U

The degree distribution of the graph describes the distribution of the number of partners of individuals of the population. Now we can define concurrency of partnerships in the following way. First we define a function XE as if there exists a v ~ V with v ~ e i and v ~ ej,

(9)

otherwise. Two edges e i and ej are concurrent in G if and only if x e ( e i , e ) = 1. We will also make use of the concept of the line graph of G, which we denote by L ( G ) . The line graph is also called an interchange graph in the literature (see, e.g., [27-29]). The vertices of L ( G ) are the edges of G, that is, the elements of E. The function Xe defines the edges of L ( G ) in the sense that (ei, e ) ~ E × E is an edge of L ( G ) if and only if x e ( e i , e ) = 1. We denote the set of vertices of L ( G ) by W and its set of edges by F; that is, L ( G ) = ( W , F ) . By definition, W = E and F = {(ei,ej) ~ E X E: xE(ei,ej)

= 1}.

(10)

CONCURRENCY AND EPIDEMICS

171

The number of concurrent partnerships in the population determines the structure of the line graph; for a completely monogamous population the line graph consists only of isolated vertices. The mean degree of the line graph increases with increasing numbers of concurrent partnerships. We want to study the relationship between the degree distributions of the graph and its line graph. 2.3. MEASURING CONCURRENCY The number of edges in the line graph is a measure of the number of concurrent partnerships in the original graph. However, this should be related to the number of vertices in the line graph. Let us first assume that N and P are given and fixed. In other words, the population size and the number of partnerships in the population are constant. The number of vertices of the line graph is P. A measure of concurrency is the number of edges of L ( G ) normalized by dividing by the maximal number of edges of L(G). We define a concurrency measure K 2 a s

IFI

(11)

This is just the density of the graph L(G). The measure K2 takes values between 0 and 1 and describes the fraction of the maximal density taken by the graph. This fraction is computed within the class of graphs with a given number of nodes IWI of the line graph, that is, a given number of edges in the contact network. In other words, the values of K2 are comparable only for graphs with a fixed given number of partnerships. This measure is discussed to some extent in [30]. In this article we focus on the following measure K3, which we call the index of concurrency: K3 = [FI/IWI.

(12)

In other words, K3 is the mean number of concurrent partnerships per partnership in G. It differs from the degree mean of the line graph by a factor 1/2. One can easily check that for a graph G with N nodes, K3 decreases if one adds edges that connect two previously isolated vertices; furthermore, if the graph G is a chain and one adds one link to the chain, then K3 increases. If one adds isolated vertices to G, K3 does not change.

MIRJAM K R E T Z S C H M A R AND MARTINA MORRIS

172

In the following theorem we use the mean m and the variance S z of the degree distribution, which are given by 1

m=~

N

N

E deg(vi)

and

$ 2 = E [ d e g ( v i ) - m ] z.

i=1

(13)

i=1

THEOREM 1 Let G be a graph with N nodes and a degree distribution (ri)i= 1..... N 1, where ri denotes the fraction of nodes of degree i. Let m be the mean and S 2 the variance of the degree distribution. Then S2

K3 = -~- + m - 1. Proof.

(14)

It can be shown that the following equation holds (see [29]): 1

N

PF]= 7 •

N

]~ X v ( v , , v j ) [ d e g ( v j ) - l ] .

(15)

i=lj=l

This is equivalent to 1

IFI = 5

N

~ i(i-1)Nri'

(16)

i=1

and so using IWI = m N / 2 we get 1

N

K3= ~- ~ i ( i - 1 ) r / .

(17)

i=1

We define the generating function for the frequency distribution ri as N

g ( z ) = E ri Zi.

(18)

i=1

Then m=g'(1)

and

S 2 =g"(1)+g'(1)-[g'(1)]

2.

(19)

We can express % as tt

K3 = ~-g (1),

(20)

CONCURRENCY AND EPIDEMICS

173

and so we get 1 K3 = ~ - ( S 2 - m + m

S2

2) = -~-- + m - 1.

(21)

As N --+ ~, the degree of the nodes of a random graph can be viewed as independent realizations of a random variable X with values in N 0. Then P [ X = i ] - - P i is the probability that a node has degree i. We denote by /* the expectation E ( X ) and by o 2 = E ( ( X - / , ) 2 ) the variance of X. The line graph of an infinite graph can be defined locally, but the numbers of nodes and edges are not necessarily finite. So in defining the index of concurrency for an infinite graph we make use of the following COROLLARY

Let G be an infinite graph with a degree distribution described by a random variable X with values in N 0 with mean E ( X ) =/* and variance 0 -2 = E((X - /*)2). Let G(N) be a graph with N nodes, and let ri(N) be the fraction of nodes of degree i. Assume that lim ri(N ) =Pi = P [ X = i ]

for all i ~ N 0 .

N - - ~ ~a

Define the index of concurrency of G as = lim K 3 ( N ) = lim IF(N)I IW(N)I "

Then 0-2

K3 = - - i f - + / * - 1 .

Proof

(22)

Let g(z) be the generating function of X,

g ( z ) = E Pi zi.

(23)

i=0

Then /*=g'(1)

and

0 - 2 = g " ( 1 ) + g ' ( 1 ) - [ g ' ( 1 ) ] 2.

(24)

174

MIRJAM KRETZSCHMAR

AND MARTINA MORRIS

Using Equation (17) we get zc

IF(N) I

1 • i(i-1)ri(N ) . K3(N) = IW(N)I = m ( N ) i = 0

(25)

For N ~ ~ we have m ( N ) --*/z, and so :~

g,,

~ 3 = 1 y, i ( i _ l ) p i /zi= °

(1) /x

(26)

This gives O-2 = ~ K 3 + ~ -- ~ 2 ,

(27)

and so O-2 =

+

-

1.

(28)

By replacing tz and or 2 in Equation (28) by the mean /2 and the variance 9 2 found in a sample of N nodes from G, we define an estimator k3 for ~3, 9 2

k3=--~-+/2-1,

(29)

with d-2 defined as the sample variance

6.2 : N--------T N E ( i - / 2 ) 2 r i,

(30)

i

where r i is now the fraction of nodes of degree i in the sample. Note that k3 is different from K3 for a graph with N nodes. If the degrees of two adjacent nodes in G are independent conditional on their adjacency, then we can derive an expression for the generating function of the degree distribution of the line graph L(G). THEOREM 2

Let G and L( G ) be infinite graphs with degree distributions described by the random variables X and Y, respectively. Let g( z ) denote the generating function of X and h(z) the generating function Y. Let tz = g'(1). If the degrees of two adjacent vertices in G are independent conditional on their

CONCURRENCY AND EPIDEMICS

175

adjacency, then h(z)=(~-)

2.

(31)

Proof. Let o¢

~c

g(z) = E rezi

and

i=0

h(z)

= E Si Ze.

(32)

i=0

The probability that an individual of degree k has a partnership with an individual of degree j is

krk(~) lz analogous to a proportional mixing assumption. For all i we can then write s i as i+1

s i = Y'=1 ~ 7kr k -i-+- 2- T - k- - r +2_k.

(33)

So .

1

i+1

SiZ,=___~ y l=k r k z k - l ( i + 2 _ k h _ Iri+2-kZ i+l-k

(34)

and sS =

(35)

i=O It is now easy to see that

K3 =

h'(1)/2.

•

The question arises whether higher moments of Y also have an interpretati__on in terms of properties of G. If the degrees of adjacent nodes in G are not independent (conditional on adjacency), then their correlation influences the higher moments of Y. One can expect an increase in the variance of the degree distribution of the line graph with increasing correlation of the degrees of adjacent nodes in G. This will be discussed in more detail elsewhere. We conclude that the index of concurrency K3 of a (finite or infinite) graph G is characterized by the mean and the variance of the degree distribution of G. As K3 is insensitive to the number of isolated vertices in G, a more complete picture is obtained by using both K3 and K1.

176 3.

MIRJAM KRETZSCHMAR AND MARTINA MORRIS C O N C U R R E N C Y AND P A R T N E R S H I P DYNAMICS

We want to investigate the impact of concurrent partnerships on the spread of an infectious disease. Therefore it is necessary to consider the formation and dissolution of partnerships and the transmission of disease within partnerships as dynamic processes. We present simulation results that have been obtained with a stochastic pair-formation model developed by Kretzschmar et al. [22]. The basic idea of the model is a generalization of a pair-formation process as described in a deterministic framework by Dietz and Hadeler [17]. The generalization is in two respects: (1) The pair-formation and disease transmission processes are fully stochastic; (2) concurrent partnerships are possible. The first feature can be described on a population level; for the second we have to consider individuals in the graph theoretical sense. In the following we first describe the model dynamics on the population level (Section 3.1), then we show how different contact structures can be realized within the population level framework (Section 3.2). 3.1. POPULATIONLEVEL DYNAMICS We first consider the process of pair formation and separation without disease. We consider a population of N individuals who can be either single or paired. Let X be the number of singles and P the number of pairs at time t; so N = X + 2P. For simplicity we assume that N is even. Possible transitions are (X,P) ~ (X-2,

P +1)

with probability pXAt

P-1)

with probability crPA t

if a new pair is formed and

(X,P)~(X+2,

if an existing pair splits up. For this process the stationary distribution of the numbers of pairs and singles can be computed (see Appendix). The expected number of pairs at any point in time is

P* = pN / 2 ( p + o').

(36)

If we include the spread of the disease, the system can be in many more different states. The state space is then eight-dimensional, and we have to describe transitions of the form ( Yo ,

, Zo,

, Poo,

, Pxo , e . )

( Y/, , r ; , z'o ,

, p(,o,

, e o, P h ) ,

CONCURRENCY AND EPIDEMICS

177

where the index 0 refers to susceptibles, the index 1 to infecteds. Y~ is the number of female singles, Z i the number of male singles, and Pij the numbers of pairs of different types. So X = Y0 + I11 + Z0 + ZI and P = Ei,j Pij. Transmission of the disease, for example, is described by (Yo, Y1, Zo, Z1, Poo, Po,, Plo, Pll) ~(Yo,Y1,Zo,Z1,Poo,Pol-

1,Plo, Pll +1)

and (Yo, II1, Zo, Zl, Poo, Pol, Pao, PI~ ) -~ (Yo,Y1,Zo,Z,,Poo,Po~,P,o - 1 , P ~ +1). These transitions have probability APolAt and APloAt, respectively. That is, we assume that within a partnership sexual intercourse takes place with a certain frequency and there is a constant probability per contact of transmitting the infection. As we have formulated it here, the model is symmetric with respect to males and females, but that can easily be generalized to asymmetric situations such as the case that transmission probabilities differ between male to female and vice versa or the case that the sex ratio is not 1:1. In this article, however, we restrict ourselves to the simplest situation. The stochastic pair-formation model with disease transmission can be approximated by a deterministic pair-formation model. We denote with X 0 and X 1 the numbers of susceptible and infected singles, and with P0o, Pop and P~ the numbers of pairs consisting of two susceptibles, one susceptible and one infected, and two infecteds, respectively. Then a deterministic model is given by the equations d -diXo = - p X o + o(2Poo + Pol ),

(37a)

d --d-iX1 = - p X 1 + o(2P11 + Pol),

(37b)

d

Xo

-diPoo = ~ P Xo + X 1 d pol

d

X° X1 = P-Xo + X1

1

x?

-d-iP11 = 2 P Xo + X~

°'Poo, o-Pol - APol , O'Pll + APol.

(37c) (37d) (37e)

MIRJAM KRETZSCHMARAND MARTINA MORRIS

178

We denote with P = P0o +/°01 d- PI1 the total number of pairs and with X = X 0 + X~ the total number of singles. Then the total population size is N = X + 2P. One can easily check that in steady state the numbers of pairs P* and singles X* are given by

P*

pN

2(p+o')

and

X*-

o'N

p+o'"

(38)

The epidemic in a population where the pair-formation process is at equilibrium can be described by a three-dimensional system of equations: d

-diX1 = -pX1

dpol dP11

= P

+

(X* -- X 1 ) X 1

1 X~

= ~p~-~-

(39a)

°'(2Pll + POl),

X*

- o'P11 +

- °'P°l -

AP0~.

hP°l'

(39b)

(39c)

As we have included neither demographic flow through the population nor recovery from the infectious state, the basic reproduction ratio R 0 is always greater than 1 and the epidemic will spread until everybody is infected. For a model with demographic flow and recovery, the basic reproduction ratio can be analyzed as in [18, 31]. The dynamic behavior of general models of this type has been analyzed in [17, 32]. 3.2. I N D I V I D U A L L E V E L D Y N A M I C S To formulate the possible transitions in the simulation model we have to describe the processes of pair formation, pair separation, and disease transmission on an individual level. Again we consider a population of N individuals, N even, with a male-to-female ratio of 1:1. An individual of the population at time t is characterized by three variables: (1) disease status (susceptible or infected); (2) degree (number of partners at time t); and (3) a vector containing identities of the partners. With this information both the contact graph and the state of the epidemic at time t are uniquely determined. An individual can undergo three types of transitions: forming a new partnership, transition from

CONCURRENCY AND EPIDEMICS

179

susceptible to infected, and separation of an existing partnership. Consequently, in every time step three types of events occur: (1) Two randomly chosen singles (one male, one female) form a pair with probability p. This is repeated n times. (2) In every pair consisting of a susceptible and an infected, the disease is transmitted with probability A. (3) Every pair splits up with probability tr. The number n that determines how many new pairs are formed per time step will depend on the number of singles or on the number of partnerships already present at time t. In a strictly monogamous population we will choose n = X(t)/2, which means that every single individual has a probability p per time step to form a new partnership. The maximum number of partnerships in this situation is N/2. This corresponds to the population level process described in the previous section. The situation changes when we abandon the assumption of monogamy. Then the process of formation of a partnership is not restricted to singles; individuals who have partners already also participate in the process. Consequently many more partnerships are possible at one time. Although in the case of serial monogamy the maximum number of partnerships is N/2, if we allow for, say, dma x c o n c u r r e n t partnerships per individual the maximum is Ndma x / 2 . However, as we want to compare serial monogamy and situations with concurrent partnerships independently of the mean number of partners per individual, we will restrict the maximum number of partnerships in all cases to N/2. This will be done by setting n = N / 2 - P(t) in all simulations. For a monogamous population, N / 2 - P(t) = X(t) as before. For a nonmonogamous population, we consider the same number of partnerships as in the monogamous case, but they can be distributed differently. This means that in general the fraction of singles will increase compared to the monogamous case. This is a somewhat artificial construction but serves to investigate the following question. Given a certain number of partnerships in a population of N individuals, how does the spread of an epidemic depend on the distribution of these partnerships among individuals? Apart from a totally random distribution where the relational status of an individual has no influence at all on his propensity to form a new partnership, there are intermediate cases where the presence of a partnership decreases the probability that this individual will participate in another partnership but will not decrease it to 0. This will be expressed with a mixing function ~b. We use the term mixing here in the sense that it describes the probability of two individuals forming a partnership depending on their degrees. On an encounter between a

180

MIRJAM KRETZSCHMAR AND MARTINA MORRIS

random female x and a random male y a partnership will be formed with probability c~(x,y) that depends on the degrees of x and y. We distinguish the following situations.

Serial monogamy

~b(x,y) =

1 0

if deg(x) = deg(y) = 0, otherwise.

(40)

Random mixing

c~(x,y) =

1 1 - ~c

if deg(x) = deg(y) = O, otherwise.

(41)

Here the parameter ~ ~ [0,1] describes a continuous transition between serial monogamy ( ~ = 1) and a pure random distribution of partnerships over the population ( ~: = 0). In the literature there has been much discussion about the influence of different types of mixing on the spread of an infectious disease [9, 10, 33, 34], and in particular assortative and disassortative mixing have been investigated (see, e.g., [1]). In assortative mixing, individuals with many partners mix with individuals who also have many partners and lowactivity individuals mix among themselves. In disassortative mixing, high-activity individuals mix with low-activity individuals and vice versa.

Assortative mixing da(x,y)=(1-~)+~deg'x'deg'y'( ~

dZa×

t ~

(42)

Here the parameter ~: ~ [0,1) describes a continuous transition between a pure random distribution of partnerships over the population ( ~ = 0) and situations where the probability of forming a partnership increases with the number of partners individuals already have.

Disassortative mixing ~b(x,y) = ( 1 - ~:)+ ~:(deg(x) - deg(y) )2 dmax

(43)

Here the parameter ~: ~ [0,1) describes a continuous transition between a pure random distribution of partnerships over the population (~ = 0) and situations where the probability of forming a partnership increases with the difference in the numbers of partners individuals already have.

CONCURRENCY AND EPIDEMICS

181

Note that in Equations (42) and (43) for sc > 0 two singles will have a probability of forming a partnership that is below 1, in contrast to the mixing function in (41). The function (41) leads to underdispersed distributions (i.e., variance-to-mean ratio below 1), while (42) and (43) lead to overdispersed distributions. The mixing functions given above should not be interpreted as behavioral and mechanistic descriptions of pair-formation processes. They are only a mathematically convenient way to produce different contact patterns in the model population. If one wants to describe mixing on a mechanistic basis, one either has to incorporate a sociologically motivated mixing function on the level of individuals or one has to use a phenomenological description that reconstructs mixing matrices from observed distributions of numbers of partners as discussed by among others, Altmann and Morris [13], Blythe et al. [3], Hsu Schmitz and Castillo-Chavez [35], and Morris [11]. Including the mixing functions 05 we can now formulate the transition rules per time step that were used for the simulations presented in the next section: (1) (a) A new partnership is formed with probability p. The individuals participating in the pair are chosen according to the mixing function 05: (i) Draw two random individuals, one male and one female; (ii) decide whether they form a pair according to 05(x,y); (iii) if yes, done; else go to (i). (b) Repeat (a) n = N / 2 - P times. (2) In every pair consisting of a susceptible and an infected, the disease is transmitted with probability A. (3) Every pair splits up with probability tr. The simulation strategy is as follows. First a contact network is built up following the rules for pair formation and separation. Once the pair-formation process is stationary in the sense of Section 3.1, the infectious disease is introduced by one index case; by definition this is at time t = 0. The simulation is continued until t = T. As we have excluded demographic processes and possible recovery, our strategy will be to investigate the spread of the epidemic over a relatively short period of time. 4.

C O N C U R R E N C Y AND EPIDEMICS

We now discuss the relationship between the structure of the partnership network and the resulting spread of the epidemic. The simulation results presented in this section were obtained with the set of parameters given in Table 1. The simulation time step was 1 day. The parameter values are not chosen on an empirical basis but are chosen such that the qualitative

182

MIRJAM KRETZSCHMAR AND MARTINA MORRIS TABLE 1 Parameter Values Used in Simulations Population size Pair-formation probability Pair separation probability Probability of transmission Simulation time

N = 2000 p = 0.01 per day tr = 0.005 per day h = 0.05 per day T = 1825 days

features of the dynamics of the system can be demonstrated. However, the orders of magnitude are within a sociologically and biologically meaningful range. Average partnership duration is 6 - 7 months, and the interpartnership period about 3 months. If we assume that the frequency of sexual contact is 1 per day, then the transmission probability per contact is 5%. The number of partnerships in the population is on average around 666 in all simulations. In the case of serial m o n o g a m y this means that about 67% of the population are paired at any point in time. The simulation time T is the time the simulation runs after introduction of one infected individual into the population. As a reference point we first want to look at the results for a population with serial monogamy. This means that we set the p a r a m e t e r in (41) to 1 and we have K3 0. We performed 100 runs to get an idea of the stochastic variation. At every point in time we evaluated the minimum; the 25%, 50%, and 75% percentiles; the maximum; and the mean of the 100 runs. These statistics are plotted in Figure 2a. The epidemic grows exponentially, and the mean epidemic agrees perfectly with the epidemic as predicted by the deterministic model (39) with the above p a r a m e t e r values. The number of infecteds in the deterministic model is given by X(t)+ Pol(t)+2Pll(t). This is shown in Figure 3 on a logarithmic scale. In the next step we cover the interval 0 ~< ~ ~< 1 in the mixing function defined by Equation (41) with a series of simulations. Trials 1-10 of Table 2 correspond to these simulations, where each trial consists of 100 runs. The maximum number of simultaneous partnerships dmax was set to 10, but, as Table 2 shows, this maximum was never attained. Trial 1 is the case of serial m o n o g a m y discussed above; trial 10 corresponds to ~ = 0. In Table 2 the index of concurrency K3 and the degree distributions are shown for the different trials. Every n u m b e r in the table is an average over time and the 100 runs. We observe how the degree distribution shifts from a distribution concentrated on 0 and 1 to approximately a Poisson distribution. -----

CONCURRENCY AND EPIDEMICS

183

I000

I00

I0 a .......

~1

500

,,,,,,1,,,,',

. . . . . . . . . . .

,,,,,,,

1000 time (a)

2000

1500

1000

10

I

1,1

o

°

,

J

500

'

I000 time (b)

: :=

= ~

2::~:i:

-

.

.

.

1500

.

i,?,i

=

2000

FIG. 2. For 100 runs with the same parameter values we computed at every time step the minimum; 25%, 50%, 75% percentiles; the maximum; and the mean of the n u m b e r of infected individuals. These statistics are plotted on a logarithmic scale for (a) a population with serial monogamy and (b) a population with a random distribution of partnerships. ( O minimum; [] 25%; ~ median; /x 75%; <~ maximum; - mean.)

184

MIRJAM K R E T Z S C H M A R AND MARTINA MORRIS 101)0

100

10

I

I

I

i

500

1000 time

1500

2000

FIG. 3. The epidemic in the deterministic model compared to the mean epidemic in the stochastic simulation with serial monogamy. For the deterministic model we plot y = X(t)+ Pm(t)+2Plx(t). ( 0 deterministic; + stochastic mean.)

In Figure 2b the statistics of the epidemic at every point in time are shown for trial 10. Compared with the same picture for trial 1, we see that now the maximum epidemic has reached almost the total population while also the mean epidemic has reached more than half the population. On the other hand, the minimum epidemic is now smaller - - i n fact, the epidemic never takes off--whereas in the serial monogamy TABLE 2 Index of Concurrency and Degree Distributions Averaged over 100 Runs for 10 Different Mixing Scenarios Degree Trial 1 2 3 4 5 6 7 8 9 10

K3

0

1

2

3

4

5-9

10+

0 0.05 0.09 0.17 0.26 0.37 0.44 0.50 0.54 0.66

33.37 35.11 36.32 38.76 41.32 44.35 46.14 47.52 48.52 51.32

66.63 63.26 60.88 56.10 51.36 45.96 42.82 40.49 38.87 34.27

0 1.60 2.71 4.84 6.70 8.54 9.48 10.09 10.45 11.42

0 0.03 0.08 0.28 0.57 1.05 1.39 1.68 1.88 2.53

0 0 0 0.01 0.04 0.10 0.15 0.21 0.25 0.41

0 0 0 0 0 0.01 0.01 0.02 0.03 0.06

0 0 0 0 0 0 0 0 0 0

CONCURRENCY AND EPIDEMICS

185

trial there were always a few secondary cases. This also has the consequence that the mean epidemic in trial 10 differs in shape from the other statistics. Given the bimodal outcomes under concurrency, the mean is no longer a sufficient summary measure of the expected size of the epidemic. This becomes clearer when we look at the distribution of the number of infecteds at t = T for trials 1-10 in Figure 4. Here we see clearly how with increasing K3 the peak of the distribution shifts to the right and variability increases. As Figure 5 shows, the log of the average number of infecteds at t = T depends linearly on K3; that is, the number of infecteds increases exponentially with increasing concurrency index. Of course, this applies only during the exponential growth phase of the epidemic; in later phases saturation effects will prevail. In Figure 6a a rough estimate of the exponential growth rate of the mean epidemic in its initial phase is plotted against K3. To estimate the growth rate we estimated the slope of the log of the mean number of infecteds as a function of time over 100 runs as long as the number of infecteds was below Nv/-N ~ (compare, e.g., [36, 37]). This was also done for a number of trials (each consisting of 100 runs) with the assortative and disassortative mixing functions. Again dmax = 10, which is here needed as a parameter of the mixing functions (42) and (43). We see that for the same value of K3 the estimated exponential growth rates differ for assortative and disassortative mixing. In the assortative case the epidemic grows faster in its initial phase than in the disassortative case. An explanation for this can be found in the differences in network structure for the two situations. In Figure 6b we have plotted the average size of the largest network component against K3. We see that for assortative mixing the number of nodes in the largest network component is greater at the same value of K3 as for disassortative mixing. Also, if we look at the structure of the largest component (Figure 7) we observe that in the disassortative case the degree variance within the largest network component increases with increasing component size, while in the assortative case the degree variance stays approximately constant. The degree mean within the largest network component is in both cases around 2 independent of the component size. So with assortative mixing for a given value of K3 there is a relatively large and homogeneous network component through which the epidemic can spread very fast, while in the disassortative case the largest network component is much more heterogeneous with respect to degree distribution, which slows the epidemic spread. The differences in the epidemics for assortative and disassortative mixing are analyzed statistically in more detail in [30].

186

MIRJAM KRETZSCHMAR AND MARTINA MORRIS 6

50

0

400

800

1200

1600

0

400

800

1200

1600

50

400

800

1200

1600

5oI

2000

o ~ 0

2

400

3]50 0

2000

400

400

2000

0

400

800

1200

1600

50 ~

0

0

1200

1600

2OOO

80O

1200

1600

2000

81 800

1200

1600

2OO0

~

9

1600

2000

41 o l__t__~ , . ~

o~

800

2000

5

,400

8(]0

1200

1600

2000

o

400

T 800

q

~

1200

50

0

10

0

400

800

1200

1600

2000

FIG. 4. The distribution of the number of infecteds after 5 years in 100 runs for trials 1-10 of Table 1.

CONCURRENCY AND EPIDEMICS

187

8.0

7.0

6.0 _Q

5.0

4"00.0

0.2

0.4 kappa 3

06

0.8

FIG. 5. The log of the mean number of infecteds after 5 years plotted against the

index of concurrency K3. (0 simulation; -- regression.)

We conclude that the concurrency index alone cannot capture effects of degree correlation, which was to be expected from the definition of K3. On the other hand, it does seem to describe the effects of degree distribution on epidemic spread for situations without assortative or disassortative bias. 5.

DISCUSSION

In the foregoing sections we have investigated the effect of concurrent partnerships on the transmission dynamics of a sexually transmitted disease. More specifically, we have studied the following question. Given a certain number of partnerships in a constant population, how does the distribution of partnerships among the individuals of the population influence the dynamics of the epidemic? Concurrent partnerships, although possibly leading to complex network structure, are described by the degree distribution of the contact graph G, that is, by local network information. This means that a considerable amount of detail about the structure of the network can be inferred from local data if they include information about the number of partners at one point in time rather than only the number of partners accumulated over a longer time period. This raises the question of how

188

MIRJAM K R E T Z S C H M A R A N D M A R T I N A MORRIS

0.10

o

0

0.05 0

0

0

0

0

OO

0:5

0.0

1:0

'

115

20

kappa 3 (a) 100

o o 50

o o

o

o

tP" * o oo o o

o

o

o

i

0.0

015

110

1.5

kappa 3 (b) FIG. 6. (a) The estimated rate of exponential growth in the initial phase o f the epidemic plotted against the index o f concurrency K3. (b) The average size o f the largest network c o m p o n e n t plotted against the index o f concurrency K3. T h e average is taken over time and over 100 runs. ( © random; [] assortative; @ disassortative.)

CONCURRENCY AND EPIDEMICS

189

4.0

2.0

0.0

(9ooOO

o

o~Kla

Q

[] O

oo

I

0

100

50

#nodes (a) 4.0

0

"~> 2.0 0

d~° ooO o 0.0

0

O

OnI~ °

0

0

0

50

100

#nodes (b) FIG. 7. The degree mean (a) and degree variance (b) within the largest network component is plotted against the ~ize of the component. (© random; [] assortative; disassortative.)

a sexual partnership is defined other than by the actual sexual act. In retrospect one can use the definition that it constitutes the time between the first and the last sexual contact; if one does not want to go into the details of the frequency distribution of sexual acts within a partnership, o n e would then assume a certain constant frequency o f sexual contacts between the first and last contact. E x t r e m e cases like two persons having two sexual contacts years apart and no contacts in

190

MIRJAM KRETZSCHMAR AND MARTINA MORRIS

between should then be described as two independent casual contacts. Given the dates of first and last sexual contact, one can then compute the time that two or more partnerships were overlapping. To measure the amount of concurrency in a population we introduced the index of concurrency K3 based on a graph theoretical approach. The index of concurrency can be computed from the mean and variance of the degree distribution of the population. This means that it can be computed from local or "egocentric" network data without any knowledge of the further network structure. It is interesting to note the similarity of the index of concurrency K3 with a quantity introduced by May and Anderson [38] as "effective contact number." These authors considered the distribution of the number of new partners per unit time. The effective contact rate is the sum of the mean and the variance-to-mean ratio of this distribution. The effective contact rate then determines the basic reproduction number R 0. In contrast, we consider the distribution of numbers of partners at one moment in time rather than accumulated over some time period. We then introduced a stochastic simulation model with dynamic partnerships. On the population level the model is a stochastic version of a model first introduced by Dietz and Hadeler [17], and on the individual level it provides full network information at any point in time. To investigate the effects of concurrent partnerships we compared a series of simulations that differ only in the distribution of partnerships over the population. In all simulations the duration of partnerships and the mean number of partners per individual were kept constant, but the degree distribution of the contact graph was varied. Outcome variables were studied as a function of the concurrency index. Although the constant mean degree across simulations allowed a comparison of the simulation results based only on higher moments of the distribution, it should be noted that it introduced the artificial effect of increasing the fraction of singles in the population. The effects of concurrency turn out to be striking. With increasing concurrency the epidemic grows faster in its initial phase and the average number of infecteds after a given time increases. The variability of the distribution of the number of infecteds also increases with more concurrency. The reasons can be found in the changing network structure. The average size of the largest network component that can be observed at any point in time increases from 2 in a monogamous population to larger connected components in a purely random distribution of partnerships. With assortative or disassortative mixing, the size of the largest component increases even more. At the same time

CONCURRENCY AND EPIDEMICS

191

the degree mean and variance within the largest network component increase. These changes in network structure imply a change in the distribution of infection risk in the population, but this heterogeneity is not based on individual attributes. With increasing concurrency, the overall variance in the number of partners per person increases, making the exposure to infectious material also more unevenly distributed. In the case of assortative mixing this is more pronounced, because risk is aggregated within a closely linked subgroup. In contrast to the conventional definition of a core group, however, this high-risk group is not a well-defined group of individuals that stays constant over time in its composition. Instead, heterogeneity in our model is a momentary concept that is not based on individual characteristics that are constant over time. The variability in individual risk in the long run is determined only by the timing of an individual's high-risk behavior within the time course of the epidemic. Our results show clearly that epidemic spread depends not only on the number of partners an individual accumulates over time but also on whether these partnerships exist simultaneously or sequentially in time. More generally, we conclude that the influence of behavioral patterns such as concurrency of partnerships on the spread of HIV/AIDS might have been underestimated up to now (see also [39]). In this respect much needs to be done in modeling culturally specific sexual network structures and their consequences for the spread of STDs. We are currently extending our framework to include more complex situations. One interesting case is the gender asymmetric situation, where one gender is monogamous and the other is not. Another important extension is to include different types of partnerships such as steady and casual, and to assume differences in behavior in those types of partnerships. It will then be interesting to study epidemic spread with realistic estimates of the parameters for different populations. Furthermore, we are extending the graph theoretical approach for describing important structural characteristics of the contact network in terms of line graph properties to more general situations like assortative and disassortative mixing. We are also investigating how to derive deterministic models that can deal with concurrent partnerships in order to obtain analytic results about the relationship between network structure and the spread of infectious disease. APPENDIX We compute the stationary state for the stochastic pair-formation process described in Section 3.1.

192

MIRJAM K R E T Z S C H M A R AND MARTINA MORRIS

We denote by qt[P] the probability that there are P pairs at time t. Then we can write

qt+~t[P] = q , [ P ] ( 1 - p~AtX -o-PAt) + qt[P-1]p( X+ I)At+ qt[P+ I]O-( P+ I)At. We now replace

X/2

by

N / 2 - P. We get

qt+at[P]-qt[P] = _qt[p][p (_2_N p)+o.p] At + qt[ P-1]p[ N -( p-1)] + qt[ P + l]o'( P + l). We now define the generation function

N/2

g(t,z) := E q,[P] ze.

(44)

P=0

For At ~ 0 we then get N

d

~ig( t,z) = - p ~ g ( t,z) + ( p - o-) z ~ g ( t,z)

N d + pTzg(t,z) - pz 2d-~g(t,z)+ o---d~g(t,z), and after manipulation, d N d g ( t , z ) = ( 1 - z)(( pz + o-)--d-~g(t,z)p--~g(t,z)).

(45)

We can now compute the stationary distribution ~(z) as solution of ~z ~ ( z )

2(

pN pz + o-) ~,(z).

(46)

Integration results in

g(z)=~

[ pz + o'tN/2 p+o- j

.

(47)

CONCURRENCY AND EPIDEMICS

193

This is the g e n e r a t i n g function o f a b i n o m i a l d i s t r i b u t i o n with p a r a m e t e r p / ( p + ~ ) . T h e e x p e c t e d n u m b e r o f p a i r s in t h e s t a t i o n a r y state is given by t h e derivative ~ ' ( z ) e v a l u a t e d at z = 1: g'(1)

pN 2( p + o-) "

(48)

T h e v a r i a n c e can b e c o m p u t e d as g"(1)+g'(1)-[g'(1)]2=-~

(0+o.)2

•

(49)

REFERENCES 1 R. M. Anderson, S. Gupta, and W. Ng, The significance of sexual partner contact networks for the transmission dynamics of HIV, J. AIDS 3:417-429 (1990). 2 R.M. Anderson, G. F. Medley, R. M. May, and A. M. Johnson, A preliminary study of the transmission dynamics of the human immunodeficiency virus (HIV), the causative agent of AIDS, IMA J. Math. Appl. Med. Biol. 3:229-263 (1986). 3 S.P. Blythe, C. Castillo-Chavez, J. S. Palmer, and M. Cheng, Toward a unified theory of sexual mixing and pair formation, Math. Biosci. 107:379-405 (1991). 4 S. Busenberg and C. Castillo-Chavez, A general solution of the problem of mixing in subpopulations and its application to risk- and age-structured epidemic models, IMA J. Math. Appl. Med. Biol. 8:1-29 (1991). 5 C. Castillo-Chavez, S. Busenberg, and K. Gerow, Pair formation in structured populations, in Differential Equations with Applications in Biology, Physics, and Engineering (Lect. Notes Pure Appl. Math. 133), J. A. Goldstein, F. Kappel, and W. Schappacher, Eds., Marcel Dekker, New York, 1991, pp. 47-65. 6 S. Haraldsdottir, S. Gupta, and R. M. Anderson, Preliminary studies of sexual networks in a male homosexual community in Iceland, J. AIDS 5:374-381 (1992). 7 H.W. Hethcote and J. A. Yorke, Gonorrhea Transmission Dynamics and Control (Lect. Notes Biomath. 56), Springer-Verlag, New York, 1984. 8 J.M. Hyman and E. A. Stanley, Using mathematical models to understand the AIDS epidemic, Math. Biosci. 90:415-473 (1988). 9 J.A. Jacquez, C. P. Simon, J. Koopman, L. Sattenspiel, and T. Perry, Modeling and analyzing HIV transmission: the effect of contact patterns, Math. Biosci. 92:119-199 (1988). 10 J . A . Jacquez, C. P. Simon, and J. Koopman, Structured mixing: heterogeneous mixing by the definition of activity groups, in Mathematical and Statistical Approaches to AIDS Epidemiology (Lect. Notes Biomath. 83), C. Castillo-Chavez, Ed., Springer-Verlag, New York, 1989, pp. 301-315. 11 M. Morris, A log-linear modeling framework for selective mixing, Math. Biosci. 107:349-377 (1991). 12 L. Sattenspiel, Population structure and the spread of disease, Hum. Biol. 59:411-438 (1987). 13 M. Altmann and M. Morris, A clarification of the 4~-mixingmodel, Math. Biosci. 124:1-7 (1994).

194

MIRJAM KRETZSCHMAR AND MARTINA MORRIS

14 P. Grassberger, On the critical behavior of the general epidemic process and dynamical percolation, Math. Biosci. 63:157-172 (1983). 15 K. Wu and R. M. Bradley, A percolation model for venereal epidemics. I. Mean field theory, J. Phys. A.: Math. Gen. 24:2569-2580 (1991). 16 K. Dietz, On the transmission dynamics of HIV, Math. Biosci. 90:397-414 (1988). 17 K. Dietz and K. P. Hadeler, Epidemiological models for sexually transmitted diseases, J. Math. Biol. 26:1-25 (1988). 18 M. Kretzschmar, J. C. Jager, D. P. Reinking, G. van Zessen, and H. Brouwers, The basic reproduction ratio R o for a sexually transmitted disease in a pair formation model with two types of pairs, Math. Biosci. 124:181-205 (1994). 19 R. Waldst~itter, Pair formation in sexually transmitted diseases, in Mathematical and Statistical Approaches to AIDS Epidemiology, C. Castillo-Chavez, Ed., Springer-Verlag, Berlin, 1989, pp. 260-274. 20 X. Luo and C. Castillo-Chavez, Limit behaviour of pair formation for a large dissolution rate, J. Math. Syst. Estimation, Control 3:247-264 (1993). 21 M. Kretzschmar, Deterministic and stochastic pair formation models for the spread of sexually transmitted diseases, Presented at the 2nd European Conference on Mathematics Applied to Biology and Medicine, Lyon 1993. 22 M. Kretzschmar, D. P. Reinking, H. Brouwers, G. van Zessen, and J. C. Jager, Network models: from paradigm to mathematical tool, in Modeling the AIDS Epidemic: Planning, Policy, and Prediction, E. H. Kaplan and M. Brandeau, Eds., Raven Press, New York, 1994, pp. 561-583. 23 K. Dietz and D. Tudor, Triangles in heterosexual HIV transmission, in AIDS Epidemiology: Methodological Issues, N. P. Jewell, K. Dietz, and V. Farewell, Eds., Birkh~iuser, Boston, 1992, pp. 143-155. 24 C.H. Watts and R. M. May, The influence of concurrent partnerships on the dynamics of HIV/AIDS, Math. Biosci. 108:89-104 (1992). 25 M. Altmann, SIR epidemic models with dynamic partnerships, J. Math. Biol. 33:661-675 (1995). 26 P. Blanchard, G. F. Bolz, and T. Kriiger, Modeling AIDS epidemics or any venereal disease on random graphs, Lect. Notes Biomath. 86:104-117 (1990). 27 C.W. Marshall, Applied Graph Theory, Wiley-Interscience, New York, 1971. 28 O. Ore, Theory of Graphs, AMS Colloquium Publ. Vol. 38, Am. Math. Soc., Providence, RI, 1962. 29 F. Buckley and F. Harary, Distance in Graphs, Addison-Wesley, Reading, MA, 1990. 30 M. Morris and M. Kretzschmar, Concurrent partnerships and transmission dynamics in networks, Social Networks (in press). 31 O. Diekmann, K. Dietz, and J. A. P. Heesterbeek, The basic reproduction ratio for sexually transmitted diseases. Theoretical considerations, Math. Biosci. 107:325-339 (1991). 32 R. Waldst~itter, Models for pair formation with applications to demography and epidemiology, Ph.D. Thesis, University of Tiibingen, 1990. 33 S.P. Blythe and C. Castillo-Chavez, Like-with-like preference and sexual mixing models, Math. Biosci. 96:221-238 (1989). 34 J. M. Hyman and E. A. Stanley, The effects of social mixing patterns on the spread of AIDS, in Mathematical Approaches to Resource Management and

CONCURRENCY AND EPIDEMICS

35

36 37 38 39

195

Epidemiology (Lect. Notes Biomath. 81), C. Castillo-Chavez, S. A. Levin, and C. Shoemaker, eds., Springer-Verlag, New York, 1989, pp. 190-219. S.-F. Hsu Schmitz and C. Castillo-Chavez, Parameter estimation in non-closed social networks related to the dynamics of sexually transmitted diseases, in Modeling the AIDS Epidemic: Planning, Policy, and Prediction, E. H. Kaplan and M. L. Brandeau, Eds., Raven Press, New York, 1994, pp. 533-559. F. Ball and P. Donnelly, Branching process approximation of epidemic models, Proc. 2nd World Congress of the Bernoulli Society, Uppsala, 1990. J. Jacquez and P. O'Neill, Reproduction numbers and thresholds in stochastic epidemic models. I. Homogeneous population, Math. Biosci. 107:161-186 (1991). R. M. May and R. M. Anderson, Transmission dynamics of HIV infection, Nature 326:137-142 (1987). C.P. Hudson, Concurrent partnerships could cause AIDS epidemics, Int. J. STD AIDS 4:249-253 (1993).