Fill and Sign the Medical Certificate of Need for Medication for the Travelling Patient Form

doi 10.1098/rspb.2000.1244 A moment closure model for sexually transmitted disease transmission through a concurrent partnership network C. Bauch and D. A. Rand * Mathematics Institute, University of Warwick, Coventry CV4 7AL, UK A moment closure model of sexually transmitted disease spread through a concurrent partnership network is developed. The model employs pair approximations of higher-order correlations to derive equations of motion in terms of numbers of pairs and singletons. The model is derived from an underlying stochastic process of partnership network formation and disease transmission. The model is analysed numerically, and the ¢nal size and time evolution are considered for various levels of concurrency, as measured by the concurrency index µ3 of Kretzschmar and Morris. Additionally, a new way of calculating R0 for spatial network models is developed. It is found that concurrency signi¢cantly increases R0 and the ¢nal size of a sexually transmitted disease, with some interesting exceptions. Keywords: concurrency; moment closure; sexually transmitted diseases; sexual networks; correlation equations 1. INTRODUCTION There has recently been great interest in the role played by the pattern of sexual contacts in the spread of sexually transmitted diseases (STDs) (see Garnett 1997; Morris & Kretzschmar 1997; Watts & May 1992; Dietz & Tudor 1992). In particular, much attention has been paid to heterogeneity in the contact structure, for instance the e¡ect of the distribution of partnerships in the population. The central questions to answer concern the e¡ect of the contact structure on R0 , growth and endemicity. In this paper we use moment closure techniques to construct analytical approaches to the calculation of these quantities. The natural history of STDs mixes two dynamical processes. First, there are the partnership dynamics in which partnerships are formed and broken. This can be regarded as producing a dynamical network in which the nodes represent individuals and where edges connect individuals in a partnership. On top of this we have an infection process which is constrained by the network in that only an individual’s partners can be infected. It is the combination of these two dynamical processes that makes the mathematical analysis of such systems more di¤cult than for disease models that assume homogeneous mixing (Anderson & May 1991) or a ¢xed network (Rand 1999). A central role is played by R0 which is de¢ned as the expected number of secondary infections produced by an infected individual during his period of infectivity (Anderson & May 1991). Although this is easily calculated for mean-¢eld systems, systems in space or on networks require a more sophisticated approach. One reason is that the spread of the disease is slowed down by the screening e¡ect due to the build-up of spatial correlations between infecteds. Additionally, for such spatial systems there is no longer a simple relationship between R0 and endemicity. For STDs on dynamic partnership networks, an expression for R0 has only been calculated * Author for correspondence. Proc. R. Soc. Lond. B (2000) 267, 2019^2027 Received 10 March 2000 Accepted 7 July 2000 in the case where each individual has at most one partner and where there is no recovery (Kretzschmar et al. 1994). In more complex networks (for instance where concurrent partnerships are present) the strength of this screening e¡ect will depend upon the structure of the network and the rate at which it is changing. The analysis of this case is more formidable. Nevertheless, one of the main results of this paper is a calculation of R0 for such systems. An important aspect of the network structure is the distribution of concurrent partnerships where one of the partners has more than one partner simultaneously. Intuitively, one might expect concurrency to increase the severity of a disease because, under serial monogamy, an infectious agent must wait for the current partnership to end before spreading the infection to a new partner (Ghani 1997). This is not surprising if higher concurrency is associated with a general increase in the average number Q of partners per person. It is less clear how concurrency changes the spread if it is increased without increasing Q. To model concurrent partnerships Dietz & Tudor (1992) created deterministic models by extending conventional pair models to allow for the case where someone who is paired gains a further single partner, thus creating triples (`triangles’ in their terminology). They conclude that the introduction of concurrency in this limited way does not have a large e¡ect on the epidemic for the parameters they chose. Watts & May (1992) have constructed a deterministic model with random mixing to understand the e¡ect of concurrency on human immunode¢ciency virus (HIV) transmission. In random mixing, the number of partners an individual currently has does not in£uence the probability of gaining further partners. They ¢nd that concurrency can bring about complex dynamic patterns for HIV, including the existence of two time-scales for the spread of the epidemic through a susceptible population: the fast time-scale corresponds to the spread of infection through the existing set of partnerships and the slower time-scale corresponds to spread due to the formation of new partnerships with susceptibles. 2019 © 2000 The Royal Society 2020 C. Bauch and D. A. Rand Moment closure STD model Finally, stochastic studies of concurrency have been published by Morris & Kretzschmar (1997), who have also suggested a concurrency measure µ3 (Kretzschmar & Morris 1996). The question they speci¢cally ask is: How do di¡erent ways of distributing a given number of partnerships in£uence disease spread? When the question is put in this way it is less clear a priori whether concurrency will have a large e¡ect, and what kind of e¡ect that will be. They ¢nd concurrency increases the growth rate of the epidemic, and conclude that this is due to the growth in the average size of the largest connected components as the concurrency is increased. One strength of the model of Kretzschmar and Morris is that it allows for a range of concurrent behaviour between the extremes of serial monogamy and random mixing. We consider a system similar to that of Kretzschmar and Morris, but develop an analytical approach to it. We use moment closure and pair approximation to construct our equations. This is now a welldeveloped technique for certain sorts of spatial ecologies (see Rand (1999) for a survey and the mathematical rationale). In STDs the transmission depends on a smaller number of contacts and is largely pairwise. Thus one believes that low-order correlations will be particularly signi¢cant and probably dominate in any moment closure. To the authors’ knowledge, no previous analytical models incorporate the range of concurrent patterns which ours does, and this makes it a useful tool for studying e¡ects such as concurrency. The derivation of the model is presented in Appendix A. There we derive the equilibrium structure of the partnership network, the proportion ºk of individuals with k partners, the basic di¡erential equations and their moment closure. The variables and parameters are also introduced in tables A1 and A2. The derivation of the basic di¡erential equation begins by considering a stochastic process on a network where nodes represent individuals and edges represent sexual partnerships through which a disease may spread. The creation and removal of edges corresponds to formation and dissolution of partnerships. Individuals may recover from infection and rejoin the susceptible class. Although it is possible to treat more complex situations, in this paper we restrict ourselves to the case where there is no incubation time, no mortality and only one gender, and where space is not represented explicitly. Thus the location of nodes and the lengths of edges do not matter. The partnership dynamics are determined by the rate »=N at which any two singles form a partnership (where N is the population size) the rate »³=N at which any two individuals, at least one of whom is in a partnership, enter into a partnership (04 ³4 1, ³ controls concurrency) and the rate ¼ of partnership break-up. An infection is transmitted between an infected person and a susceptible person at rate l , and infected persons recover at rate ¸ (see Appendix A for details). 2. MOMENT CLOSURE Moment closure involves, ¢rst, the derivation of di¡erential equations for certain low-order correlations and, second, a closure scheme for truncating these equations. The latter scheme usually involves replacing higher-order Proc. R. Soc. Lond. B (2000) correlations by functions of the lower-order ones and may involve modelling combinations of high-order correlations as noise. The equations of motion (A5)^(A10) for singleton and pair quantities involve the third-order quantities Q (I jSI) and Q (I jSS) which we must approximate in terms of singletons and pairs to close the system of equations. The ¢rst of these quantities Q (I jSI) is the mean number Q xy (I) of infected partners of a susceptible x, given that the susceptible already has at least one infected partner y. Similarly Q (I jSS) is the mean number of infected partners of a susceptible, given that the susceptible has at least one susceptible partner. The fact that one has only to approximate Q (I jSI) and Q (I jSS) to close the system is general for epidemiological systems and means that the analytical treatment of such systems is not as di¤cult as is usually believed. We will be interested in calculating the value of quantities such as Q (I jSS) in the context of invasion of the disease into a purely susceptible population. In this case, instead of calculating a global value for Q (I jSS), we should restrict attention to those SS pairs that are part of the invading population (i.e. one of the susceptibles is partnered with an infective). This will be seen to produce a more accurate moment closure. We denote the appropriate value of Q (I jSS) by Q (I jSS)c. For Q (I jSI) we obviously do not need to distinguish the invasion quantity from the global one because the susceptible must be partnered with at least one infective. There are a number of ways to derive an approximation for Q (I jSI). The two conventional approaches where one assumes that the quantities Q x y (I) are Poisson or multinomially distributed are untenable for the network dynamics we have de¢ned. A multinomial distribution does not work well for networks which are not regular lattices, and a Poisson distribution does not work well unless the partnership distribution is random. Thus we develop a di¡erent, rather heuristic approach which combines simplicity, intuitive accessibility and accuracy for both equilibrium solutions and transients. Before proceeding, there are two things about this moment closure which should be noted. First, it assumes that the number of partnerships of a susceptible is not correlated with their disease statuses, and that knowledge of the disease status of one partner of the susceptible does not change the expected disease status of any of the other partners (equation (2)). Second, triangles are not likely in this partnership network, unless the population size N is small. These assumptions are clearly unrealistic; in real networks we ¢nd core groups, high-activity assortative mixing, triangles, etc., but our assumptions are a useful simpli¢cation at this preliminary level of analysis on the very particular question of the e¡ect of concurrency. In this derivation, ‰ij Š and ‰ijkŠ respectively denote the number of pairs in state ij and the number of triples in state ijk and ‰¢Šc denotes the corresponding pair or triple numbers inside the invading population. It is always possible to express quantities such as Q (ij j) in terms of pairs ‰ij Š, singletons ‰i Š, etc. For instance, Q (ij j) ˆ ‰ij Š/‰ j Š. The ¢rst step in the derivation of Q (I jSI) is to note that Q (I jSI)4 1 because the susceptible has at least one infected partner. We must also determine how many other infected partners the susceptible individual has on average. The expected number of partners of an individual with at Moment closure STD model least one partner Q 5 1 can be calculated completely aside from any consideration of disease status. We get this from equation (A3). Thus the number of extra (infected or uninfected) partners given that the individual has at least one partner is Q 5 1 ¡ 1. Also, the fraction of extra partners of the susceptible in an SI pair who are infected is (Q (I jSI) ¡ 1)=(Q (I jSI) ¡ 1 ‡ Q (SjSI)). Combining this with the expression for the total number of extra partnerships Q 5 1 ¡ 1 we obtain Q (I jSI)º 1‡ (Q 5 1 ¡1)£ Q (I jSI) ¡ 1 . Q (I jSI) ¡ 1 ‡ Q (SjSI) (1) Also, we rearrange Q (SjSI) to obtain Q (SjSI) ˆ ‰SSI Š/ ‰SI Š ˆ (‰SSI Š/‰SS Šc ) £ (‰SS Šc /‰SI Š) ˆ Q (I jSS)c ‰SS Šc /‰SI Š. For Q (I jSS)c we take Q (I jSS)c º Q (I jSI) ¡ 1. (2) This follows from our independence assumption, and we note that this expression is valid only for Q (I jSS)c and not for Q (I jSS). The upshot of this will be that our moment closure makes some progress in extending the usefulness of pair approximations to invasion analysis and not just equilibrium analysis. We need a third equation. Some thought will convince the reader that Q (I jSS)c ˆ ‰ISS Š=‰SS Šc . Also Q (I jSS) º ­ Q (I jS) where ­ is a factor which comes from the fact that the susceptible already has at least one partner. ­ is the ratio of how many partners a sexually active person has, minus one for the susceptible (i.e. Q 5 1 ¡ 1), to the number of partners anyone has (i.e. Q ). In other words ­ adjusts for what we know about the network structure in this case. We combine these facts in the following derivation: ‰ISS Š Q(I jSS)‰SS Š Q ¡1 ‰SS Š ˆ Q(I jSS)c ˆ º Q (I jS) 5 1 . ‰SS Šc ‰SS Šc Q ‰SS Šc (3) Now that we have three equations in three unknowns, we can solve equations (1), (2) and (3) for ‰SS Šc, Q (I jSI) and Q (I jSS)c : ‰SS Šc ˆ ‰SI Š ‰SSŠ , Q ‰S Š ¡ ‰SS Š (4) Q (I jSI) ˆ ‰SS Š ‡ Q 5 1 (Q ‰S Š ¡ ‰SS Š) , Q ‰S Š (5) Q (I jSS)c ˆ (Q 5 1 ¡ 1)(Q ‰S Š ¡ ‰SS Š) . Q ‰S Š (6) We can check partially if these approximations have the behaviour they should. First, consider ‰SS Šc. De¢ne Q S (respectively Q I ) as the average number of partners per susceptible (resp. infected). We expect that Q S 5 Q and Q I 4 Q , because having more partners means one is more likely to be infected. It easily follows from this and the equality ‰SI Š ‡ ‰SS Š ˆ Q S ‰S Š that ‰SSŠc 4 ‰SSŠ as we would demand. Second, we note that Q (I jSS)c is positive de¢nite. From the foregoing we have Q ‰S Š ¡ ‰SS Š4 0, and since Q 5 1 ¡ 14 0 (when ³4 0) it also follows that Q (I jSS)c 4 0. Thus Q (I jSI) ˆ 1 ‡ Q (I jSS)c 4 1 implies Proc. R. Soc. Lond. B (2000) C. Bauch and D. A. Rand 2021 Q (I jSI)4 1. Third, it can also be checked that in the case of serial monogamy, Q (I jSI) ˆ 1 and Q (I jSS)c ˆ 0. Substituting approximations (4)^(6) into the master equations (A5)^(A10) gives a closed set of ¢rst- order di¡erential equations for the pair numbers ‰SSŠ, ‰SI Š and ‰II Š, and for ‰I Š, ‰XI Š and ‰XS Š. The accuracy of this moment closure is considered in Appendix B, where a comparison between the stochastic and deterministic models is made. 3. CALCULATING R0 R0 for STD models has only been derived for cases where the processes of partnership formation and separation are not taken into account (Diekmann et al. 1990) or systems where each individual can have at most one partner at a time (Kretzschmar et al. 1994; Diekmann et al. 1991). We discuss a new method which exploits the structure of invasions in systems that are spatial and/or on a network. Thus, not only does this approach allow one to attack a seemingly much more di¤cult problem, but also, we would claim, it is much closer to reality than mean¢eld theory in the way it models the invasion. This new method relies on certain characteristics of the time evolution of the spatial structure. Let us consider the situation where an infection is introduced into a purely susceptible population by randomly placing a few infecteds among the susceptibles. Moreover, we will assume that the partnership dynamics in the purely susceptible population have had time to come to equilibrium. We then consider the invasion process. A stochastic invasion of such a system consists of three phases: (i) Inoculation phase. The ¢rst stage is largely stochastic and relatively quick: a number of infective individuals (with ‰I Š ½ N ) are randomly placed in the population and start infecting their partners. (ii) Establishment phase. If the disease does not die out, then there is a second stage which is characterized by the fact that at its beginning we still have ‰I Š ½ N, but now the individual infectives have grown into small local populations with a wellde¢ned local correlation structure that can be calculated by the techniques described below. (iii) Development phase. Here, the disease may either die out or grow until it reaches some equilibrium. If the disease grows, the characteristic local correlation structure undergoes changes as the various patches of infection start to come into contact with one another. Although ¢rst identi¢ed in lattice models, this pattern is also observed in the dynamic networks investigated in this paper, even though the network structure is totally di¡erent in the two cases. At the start of the establishment phase we expect that all the quantities of the form Q (ijI) and Q (ij jI ) (where i and j denote either infected (I ) or susceptible (S)) will have reached quasi-equilibrium values. For example, Q (SjI) will go from its initial value of Q to its quasiequilibrium value. All this can be seen very nicely in ¢gure 1, which shows the time evolution of Q (SjI ). For example, one sees the rapid decline of Q (SjI ) from its Moment closure STD model 0.3 7 0.27 5.6 0.24 4.2 log {I} Q(S|I ) 2022 C. Bauch and D. A. Rand 0.21 0.18 0.15 2.8 1.4 0 600 1200 1800 time t 2400 3000 0 0 500 1000 1500 time t 2000 2500 3000 Figure 1. Time-series of Q (SjI). » ˆ 0:01, ¼ ˆ 0:005, ³ ˆ 0:3, l ˆ 0:1, ¸ ˆ 0:02. Initial conditions: ‰SIŠ ˆ ‰IIŠ ˆ ‰IŠ ˆ ‰XI Š ˆ 1. Figure 2. Time-series of log‰I Š. » ˆ 0:01, ¼ ˆ 0:005, ³ ˆ 0:3, l ˆ 0:1, ¸ ˆ 0:02. Initial conditions: ‰SIŠ ˆ ‰IIŠ ˆ ‰IŠ ˆ ‰XI Š ˆ 1. initial value of Q (SjI) ˆ Q and the approximately zero slope of Q (SjI ) after the disease has established its local structure. Also, we see the gradual change in Q (SjI) as the number of infected individuals increases. Compare this with ¢gure 2 which shows the time-series for log‰I Š for these parameters. Note that immediately after Q (SjI) has reached quasi-equilibrium, the growth of ‰I Š is exponential and Q (SjI ) remains roughly constant during the exponential growth phase, as expected. We will exploit these growth characteristics to obtain a formula for R0. Using the equation d‰I Š/dt ˆ (¡ ¸ ‡ l Q (SjI ))‰I Š, we see that at the start of the establishment phase an individual, while infective, infects others at a rate l Q (SjI)e , where Q (SjI)e is the value of Q (SjI) in this phase (i.e. the quasi-equilibrium value). The expected length of time an individual remains infective is 1=¸. Thus we deduce that Next we substitute the expressions for Q (I jSI), Q (I jSS)c and ‰SSŠc , and also make the substitutions ‰XS Š ˆ X ¡ ‰XI Š, XI ˆ O ‰I Š, ‰II Š ˆ Q (I jI )‰I Š and ‰SI Š ˆ Q (SjI )‰I Š to produce equations in terms of Q (SjI), Q (I jI), O and ‰I Š. Finally, we make the approximation that ‰I Š=N º 0 (because the number of infecteds is initially small) to produce R0 ˆ l Q (SjI)e ¸ d Q (SjI ) ˆ (»O QXN ‡ »³QN 2 ¡ »³Q O XN dt ¡ ¼Q (SjI)N 2 Q )/((N 2 Q Q (SjI )) ‡ (¸Q (I jI )N 2 Q ‡ 4l Q (SjI )NQ 5 1 P ¡ 4l Q (SjI )NP)/N 2 Q Q (SjI)) ‡ ( ¡ l Q (SjI)N 2 Q Q 51 ¡ l Q (SjI )2 N 2 Q )/(N 2 Q Q (SjI)). . Consequently, we must estimate the quasi-equilibrium value Q (SjI )e. We use the original equations of motion (A5)^(A10) to derive equations for dQ (SjI )/dt, dQ (I jI )=dt and dO =dt where O ˆ ‰XI Š=‰I Š. These give information about the local structure evolution. These equations still include the terms Q (I jSI) and Q (SjSI), but we use the pair approximations (4)^(6) derived in the previous section to approximate them. We can then employ the quasi-equilibrium approximation dQ(SjI)/dt ˆ dQ(I jI)/dt ˆ dO /dt ˆ 0 and solve for Q (SjI )e, Q (I jI )e , and O e. By taking the derivative of Q (SjI ) ˆ ‰SI Š/‰I Š we obtain dQ (SjI)/dt ˆ (‰SI Š¡1 d‰SI Š/dt ¡ ‰I Š¡1 d‰I Š/dt)Q (SjI). Substituting the expressions for d‰I Š/dt and d‰SI Š/dt from equations (A10) and (A5) produces 1 d Q (SjI ) ˆ » ‰XS Š‰XI Š ‡ »³ ((N ¡ ‰I Š)‰I Š ¡ ‰XS Š‰XI Š) N N ‰I Š dt ¡ ¼‰SI Š ¡ ¸‰SI Š ‡ ¸‰II Š ‡ l ‰SS Šc Q (I jSS)c ¡ l ‰SI ŠQ (I jSI) ¡ (l Q (SjI) ¡ ¸)‰SI Š. Proc. R. Soc. Lond. B (2000) We can derive similar equations for dQ (I jI )=dt and dO =dt. Thus we have three nonlinear di¡erential equations in three unknowns. At the point where the disease becomes established we can set dQ (SjI )=dt ˆ dQ (I jI)/dt ˆ dO =dt ˆ 0 which is a good approximation for most parameter choices. Then we solve for the unknowns, giving us the values Q (SjI)e , Q (I jI)e and O e of Q (SjI), Q (I jI) and O at establishment. The solution for Q (SjI )e is a large third-order polynomial which is di¤cult to simplify except in special cases. Because of its length we do not write the full solution here. The complication of the solution re£ects the complexity of the model rather than any unnecessary complication introduced by our method. However, we can study it numerically (see ½ 4). In the case where partnership dynamics are much slower than infection dynamics (the orders of » and ¼ much less than the orders of l and ¸), R0 becomes R0 ˆ ¡1 ¡ w ‡ ¹ ‡ (¹ ¡ w)2 ‡ 2w, (7) where w ˆ l =¸, ¹ ˆ 1=2 ‡ l P=(N¸(1 ¡ º0 )), and º0 is as in equation (A2). Moment closure STD model 3 C. Bauch and D. A. Rand 2023 1000 800 2 600 R0 {I} 400 200 1 0 0.1 0.2 0.3 0.4 concurrency 0 0.5 Figure 3. Plot of log R0 against µ3 . » ˆ 0:02, P0 ˆ 500, l ˆ 0:05 and r ˆ 0:005. 2.0 0 BP 0.1 0.2 0.3 0.4 concurrency 0.5 0.6 Figure 5. Bifurcation diagram of the number of infecteds ‰IŠ where µ3 is the active parameter. The solid line indicates the stable branch. BP indicates the branching point where the trivial branch exchanges stability with the non-trivial branch. The trivial branch to the left of the branching point is stable, and to the right is unstable. » ˆ 0:01, P0 ˆ 442:4, l ˆ 0:02, ¸ ˆ 0:001. 1500 R0 1200 900 1.0 0.9 0.1 {I} 0.2 0.3 0.4 0.5 concurrency 0.6 0.7 Figure 4. Plot of log R0 against µ3 for the special case of slow partnership dynamics (see equation (12)). » ˆ 0:001, P0 ˆ 600, l ˆ 0:12 and ¸ ˆ 0:01. Taking the derivative of R0 with respect to º0 shows that R0 is monotone increasing in º0. Also, Appendix A, ½ (b) shows that º0 is itself an increasing function of ³ on the interval 05 ³5 1 (if P is held constant by adjusting ¼). Thus, R0 is also an increasing function of ³. As the proportion of individuals in the population with multiple partners increases, R0 also grows and the disease is able to spread more quickly through the population. Concurrency always increases R0 for slow partnership dynamics. Because R0 is monotone increasing in ³, we know that R0 is bounded above and below on the interval 05 ³5 1 according to the bounds for º0 (see Appendix A, ½ (b)). It is di¤cult to infer the shape of the R0 curve on the 05 ³5 1 interval because of the complicated dependence of º0 on ³. However, numerical analysis suggests that the dependence of R0 on the index of concurrency µ3 is roughly exponential. 4. NUMERICAL RESULTS The analysis of ½½ 2 and 3 allows us to consider the dependence of R0 and endemicity on concurrency. The di¡erential equations were analysed numerically using Proc. R. Soc. Lond. B (2000) 600 300 0 0 0.1 0.2 0.3 0.4 concurrency 0.5 0.6 Figure 6. Bifurcation diagram for ‰I Š showing the evolution of the non-trivial branch as r is varied. The topmost branch corresponds to r ˆ 0:003 and the bottom branch corresponds to r ˆ 0:006. » ˆ 0:01, P0 ˆ 300, l ˆ 0:05, ¸ ˆ 0:003, 0:004, 0:005, 0:006. CONTENT 1.5. As mentioned in }1, one wants to separate the e¡ects of increasing concurrency per se from an increase accompanied by a larger equilibrium number of partnerships P. Thus in this section, P is always held at some constant value P0 as concurrency is increased. First, there is the question of how fast R0 grows with increasing concurrency as measured by µ3. Figure 3 shows a monotonic increase in R0 with concurrency, although it is not exponential; increasing P gradually reverses the curvature. Figure 4 shows the special case of very slow partnership dynamics (as calculated in equation (7)), and indicates roughly exponential growth of R0 with concurrency; it remains exponential for other values of P. Figure 5, a bifurcation diagram, shows the most typical behaviour for realistic parameters of the number of infecteds versus µ3. As concurrency is increased, a branching 2024 C. Bauch and D. A. Rand Moment closure STD model 1500 1500 1200 1200 900 900 {I} {I} 600 600 300 300 0 0 0.1 0.2 0.3 0.4 concurrency 0.5 0.6 Figure 7. Bifurcation diagram showing the case of decreasing endemicity. » ˆ 0:01, l ˆ 0:01, ¸ ˆ 0:0006, P0 ˆ 442:4. point is reached and the trivial and non-trivial solution branches exchange stability. Beyond this threshold the endemicity increases rapidly with increasing concurrency. Figure 6, which shows the evolution of the non-trivial branch as ¸ is varied, indicates that as ¸ decreases the branching point moves closer to the origin and there is a sharper increase in endemicity after the branching point. The S-shape of the non-trivial solution branch is also notable (also apparent in ¢gure 5). In some cases the endemicity can decrease with increasing concurrency. Figure 7 shows this case. Here, the endemicity is high on account of a very low recovery rate. This decrease might be caused by a tendency of the infection to cluster in certain parts of the network. If infected individuals cluster more as concurrency is increased (perhaps on account of the increased variation in the node degree distribution), then the number of SI links is reduced, and hence overall disease transmission is also reduced. This e¡ect, in the cases of high endemicity, might be strong enough to cause a decrease in endemicity for certain values of concurrency. However, these results would only have signi¢cance for diseases with a characteristically high endemicity, and where there is recovery to a susceptible state and subsequent reinfection. Figure 8 shows three examples (for three values of ¸ ) where the partnership dynamics are very fast. In this case one would not expect network structure to be important to infection dynamics. In all three cases, the endemicity is relatively constant as concurrency is increased, although the endemicity increases or decreases somewhat according to the recovery rate. Thus, even when the partnership dynamics are occurring at such a fast rate (»=¸ º 2000), network structure has some impact. It is also interesting that the endemicity can either decrease or increase with increasing concurrency. 5. DISCUSSION We have shown that it is possible to use the ideas of moment closure to derive di¡erential equations describing interesting STDs spreading on a dynamic partnership Proc. R. Soc. Lond. B (2000) 0 0 0.1 0.2 0.3 0.4 concurrency 0.5 0.6 Figure 8. Bifurcation diagram for ‰I Š for the case where partnership dynamics are fast relative to infection dynamics (» ˆ 0:2, P0 ˆ 365:7, l ˆ 0:0005, ¸ ˆ 0:0001, 0:00015, 0:0002). The topmost branch corresponds to ¸ ˆ 0:0001 and the bottom branch corresponds to ¸ ˆ 0:0002. network. These equations allow one to analyse the e¡ects of important parameters without resorting to timeconsuming simulations. In particular, they allow one to explore much more easily issues of stability, robustness and ubiquity. In certain limiting regimes one obtains simple expressions for quantities such as R0. However, there does not seem to be a single optimal moment closure scheme for such systems. Which scheme one chooses depends upon one’s particular aim. Here we have chosen a scheme which combines simplicity, intuitive accessibility and accuracy for both equilibrium solutions and growth rates. Most important is our calculation of R0 for spatial and network systems, on grounds of its novelty and the insights into local structure evolution on which it is based (however, see Rand 1999). The concept is a very natural one for such systems. We believe that this approach will be of considerable use in a wide range of invasion problems in spatial and network systems. Finally, the results show that concurrency in sexual partnership networks greatly increases the impact of a sexually transmitted disease, both in terms of ¢nal size and the basic reproduction ratio R0. One possible elaboration of this model would be to tailor it more closely to particular STDs. This could be done just by looking at particular areas of the parameter space or by introducing elements such as disease-induced mortality, age structure or core groups. Although the sexual network in this model is not sociologically realistic, we believe that the results are robust for the other types of concurrency found in real sexual networks. One way to check this would be to introduce other types of concurrency models and compare the e¡ect they have on the disease epidemiology. D.A.R. gratefully acknowledges the support of the Engineering and Physical Sciences Research Council, the Biotechnology and Biological Sciences Research Council and the University of Warwick. C.T.B. would like to acknowledge the National Science Foundation. The authors would also like to thank two anonymous referees for their comments. Moment closure STD model C. Bauch and D. A. Rand Table A1. Dynamical variables involved in moment closure a Table A2. Parameters of the model a symbol de¢nition parameter de¢nition [XS] [XI] [SI ] [II] [SS] [I ] [S ] number of single susceptible individuals number of single infected individuals number of infected ^ susceptible partnerships 2 £ number of infected ^ infected partnerships 2 £ number of susceptible ^ susceptible partnerships total number of infected individuals total number of susceptible individuals »/N »³/N a Notice also that the variables [II ] (respectively [SS]) are determined by counting each I^I (resp. S^S) edge twice, in keeping with a convention established by earlier research on correlation equations. Also, [IS ] ˆ [SI ]. APPENDIX A. MODEL DESCRIPTION AND DERIVATION (a) Variables and parameters Tables A1 and A2 show, respectively, the dynamical variables and static parameters. The following constraints apply at equilibrium: 2P ˆ 2‰SI Š ‡ ‰II Š ‡ ‰SSŠ and X ˆ ‰XS Š ‡ ‰XI Š where X and P are as in equation (A4). (b) Network structure and the index of concurrency k 3 The index of concurrency µ3 (Kretzschmar & Morris 1996) is the number of concurrent partnerships divided by the total number of partnerships, i.e. the number of triples divided by the number of pairs. It can be shown from the ºk distribution given in equation (A3) that µ3 ˆ »³ ˆ ¿³, ¼ (A1) for ³5 0. The proportion ºk of the population with k partners can be calculated as follows. Let Fm;n be the rate at which individuals move from a state of having m partners to a state of having n. In equilibrium, ºk Fk, k‡ 1 ˆ ºk‡ 1 Fk‡ 1, k . We also require that k5 0 ºk ˆ 1. Because of the appearance of singularities at ³ ˆ 0 and ³ ˆ 1, it is convenient to calculate the equilibrium of this process separately for the three cases ³ ˆ 0 (serial monogamy), ³ ˆ 1 (Poissondistributed network), and 05 ³5 1. For º0 we obtain (¡1¡ º0 (³) ˆ ¡­ ‡ e¡¿ p 1 ‡ 4¿)=2¿ ¡­ 2 ‡ 4(­ ‡ 1)(³¡1 ¡1) 2(­ ‡ 1)(³¡1 ¡ 1) ³ˆ0 05 ³5 1, ³ˆ1 (A2) where ­ ˆ exp (¿³). For ºk, k4 1: ºk (³) ˆ 0 º0 (³) (º (³) ‡ ³(1 ¡ º0(³))) (¿³)k 0 k! ³ ¿k e¡¿ =k! ³ˆ0 05 ³5 1. ³ˆ1 (A3) º0 is the proportion of single individuals in the population. In equilibrium, the number of single individuals X is N º0 , and the number of partnerships P is given by Proc. R. Soc. Lond. B (2000) ¼ l ¸ N 2025 rate at which any two singles form a partnership rate at which any two singles, at least one of whom is in a partnership, enter into a partnership (04 ³4 1). ³ controls concurrency partnership separation rate disease transmission rate in an infected^subsceptible partnership recovery rate population size a ³ ˆ 0 corresponds to monogamy while ³ ˆ 1 corresponds to independence of partnerships from one another. By scaling time and the singleton and pair numbers, one can obtain a set of equations that depend only on the parameters ¿ ˆ »/¼, ! ˆ l /¼, ¯ ˆ ¸/¼, and ³. It follows that quantities such as R0, which are coordinate independent and do not depend on the unit of time, can only depend upon such dimensionless parameter combinations. P(¿, ³, N) ˆ N 2 k5 1 kºk ˆ ¿(³N 2 ‡ (1 ¡ ³)X 2 ) . 2N (A4) For given values of N and P there is a unique value of ¿ ˆ ¿(³) such that P(¿(³), ³, N) ˆ P0 for all ³ on the interval ‰0, 1Š. However, this value can only be determined numerically. Also, it is clear that if P is held at the constant value P0 as ³ is increased, then º0 must be monotone increasing in ³, because the proportion of individuals with multiple partners goes up as ³ increases. Because of the monotonicity of º0 in ³, º0(1) and º0 (0) are upper and lower bounds, respectively, on º0 (³). (c) Derivation of equations of motion We employ the following archetypal equation for deriving the equations of motion. If f is the time t expectation value of some function of the state of the network, then f_ ˆ r(e)¢fe , e2Events where r(e) is the rate of event e, and ¢fe is the change induced in f by event e. The main thing to understand here is that we derive the equations by summing over all individuals and considering all possible events which can a¡ect the number of pairs, singles, etc. For instance, in deriving the equation of motion for ‰SI Š, consider that the process of recovery of an infected partner can cause the transition SI ! SS; this destroys one SI pair. Because the recovery rate is ¸, SI pairs are destroyed in this way at rate ¸‰SI Š. Another example is pair formation: a susceptible and an infected can form a partnership S ‡ I ! SI which increases ‰SI Š by one. Some processes involve not only the two individuals who are part of the SI pair but also individuals connected to them (as for instance with concurrent disease transmission). In these cases higherorder correlations come into play in the form of ISI and ISS triples. These must be approximated in terms of pairs and singletons so that we do not have an in¢nite hierarchy of equations. Moment closure STD model 2026 C. Bauch and D. A. Rand Let ¯x ˆ (state description) denote the state of a node. For instance ¯x ˆ (P) indicates that node x has a non-zero degree (the individual is in a partnership). Let ¯xy ˆ (state description) denote the state of two nodes which may or may not be joined by an edge. A statement such as ¯xy ˆ (S, X; I , P) indicates that node x is susceptible and single, and node y is infected and connected to an edge, but there is no edge between them. However, a statement such as ¯xy ˆ (SI) indicates that node x is susceptible, node y is infected, and there is an edge between them. Let Q x (i) denote the number of partners (edges) of type i of node x; i can be either susceptible or infected. Q (ij j) is the population-averaged number of i partners (edges) of a j. Q (ij jk) is the average number of i partners (edges) of a j, given that the j has a k partner. Note that Q (ij j) ˆ ‰ij Š/‰ j Š and Q (ij jk) ˆ ‰ijkŠ/‰ jkŠ when i 6 ˆ k and Q (ij ji) ˆ ‰ijiŠ/‰ij Š ‡ 1. We now derive the equation for d‰SI Š=dt. We sum over all sites in the network where events can occur which a¡ect ‰SIŠ, and consider the rates at which they happen, to produce the following expression: d ‰SIŠ ˆ dt ¯ xy » ‡ N ˆ (S,X ;I ,X) ¡ ‡ ¯xy ˆ (SI ) ¼¡ ¯xy ˆ (SS)c ¯xy ˆ (S,P;I ,P)[ (S,P;I,X)[(S,X ;I ,P) ¯xy ˆ (SI) l Q x(I) ¡ ¸‡ ¯xy ¡ ‡ ¡ » ‡ N ˆ (S,X ;I ,X ) ¯xy ˆ (SI ) ¯xy ˆ (SS) ¯xy ˆ (SI ) ¼¡ ¯xy ˆ (SI) ‡ l Q x (I). ¡ ¯xy ˆ (S,X;I,X) ¯xy ˆ (SI ) ¯xy ˆ (II ) ¼£ 0 1 if y is not monogamous if y is monogamous ¼£ 0 1 if y is not monogamous if y is monogamous » »³ ‰XI Š(‰XI Š ‡ ‰XS Š) ¡ ‰XI Š(N ¡ X) ¡ r‰XI Š‡ N N ¼ £ (number of monogamous infected individuals). ¡ We approximate the number of monogamous infected individuals as ‰SI Š £ (probability I is monogamous) º ‰SI Š »³ N ¸ ¯xy ˆ (II) l fQ (I jSS)c ‡ ²x (I jSS)c g l fQ (I jSI) ‡ ²x (I jSI)g. On taking the sums, the ²x £uctuations disappear because ˆ 0. Taking means of their linearity: ¯xy ˆ ( jk) ²x (ij jk) and parameters out of the sums, and evaluating the sums, produces » »³ d ‰SI Š ˆ ‰XS Š‰XI Š ‡ f(N ¡ ‰I Š)‰I Š ¡ ‰XS Š‰XI Šg ¡ ¼‰SI Š N N dt ¡ ¸‰SI Š ‡ ¸‰II Š ‡ l ‰SS Šc Q(I jSS)c ¡ l ‰SI ŠQ(I jSI). (A5) Proc. R. Soc. Lond. B (2000) d » »³ ‰X Š ˆ ¡ ¡ ¡ ¸ N N ¯ ˆX dt I ¯xy ˆ (I,X;X) ¯xy ˆ (I,X;P) x I ¸ ¯xy ˆ (S;I) (A7) For ‰XI Š the derivation is very similar, except for one term which requires us to estimate the number of monogamous individuals: ˆ ¯xy ˆ (II) (A6) d » »³ ‰SS Š ˆ ‰XS Š2 ‡ f(N ¡ ‰I Š)2 ¡ ‰XS Š2 g¡ ¼‰SSŠ N N dt ‡ 2¸‰SI Š ¡ 2l ‰SS Šc Q (I jSS)c . ‡ ¯xy ˆ (SI ) ¸‡ » »³ d ‰II Š ˆ ‰XI Š2 ‡ (‰I Š2 ¡ ‰XI Š2 ) ¡ ¼‰II Š N N dt ¡ 2¸‰II Š ‡ 2l ‰SI ŠQ (I jSI), »³ N Next we must simplify the sums which involve Q x (I ). We do this by making the substitutions Q x (I) ˆ Q (I jSI) ‡ ²x (I jSI) and Q x (I ) ˆ Q (I jSS)c ‡ ²x (I jSS)c . Q (I jSI) and Q (I jSS)c are the population-averaged means discussed in } 2, and the ²x terms represent the £uctuation from these means at node x. Making this substitution produces d ‰SI Š ˆ dt The equations for ‰II Š and ‰SS Š are derived similarly: º1 : 1 ¡ º0 Thus, d » »³ ‰XI Š ˆ ¡ ‰XI Š(X)¡ ‰XI Š(N ¡X) ¡ ¸‰XI Š N N dt º1 ‡ ¼‰SI Š . 1 ¡ º0 (A8) Similarly, for ‰XS Š d » »³ ‰XS Š ˆ ¡ ‰XS Š(X) ¡ ‰XS Š(N ¡ X) ‡ ¸‰XI Š N N dt º1 ‡ ¼‰SI Š . 1 ¡ º0 (A9) Finally, for ‰I Š d ‰I Š ˆ dt ˆ ¯x ˆ (I ) ¸‡ l ¯xy ˆ (SI ) ¡¸‰IŠ ‡ l ‰SIŠ. (A10) To obtain a closed set of equations, Q (I jSS)c and Q (I jSI) must be expressed in terms of the given parameters and variables as discussed in ½ 2. APPENDIX B. COMPARISON OF STOCHASTIC AND DETERMINISTIC MODELS A computer simulation for the underlying stochastic model was run to compare results with the deterministic model. Comparison allows us to re¢ne the deterministic model, understand the stochastic model better, and estimate the error introduced by the moment closure. Moment closure STD model C. Bauch and D. A. Rand 2027 Table B1. Parameter values (per day) for comp arison of stochastic and deterministic predictions case 1 2 3 4 5 6 ³ l ¸ 0.08 0.1 0.15 0.1 0.1 0.1 0.10 0.10 0.10 0.20 0.08 0.15 0.006 0.006 0.006 0.012 0.006 0.006 µ3 ¸/¼ case 0.16 0.2 0.3 0.2 0.2 0.2 1.2 1.2 1.2 2.4 1.2 1.2 7 8 9 10 11 12 ³ 0.08 0.1 0.15 0.1 0.1 0.1 ¸ µ3 ¸/¼ 0.003 0.003 0.003 0.003 0.003 0.0015 0.16 0.2 0.3 0.2 0.2 0.2 0.6 0.6 0.6 0.6 0.6 0.3 l 0.05 0.05 0.05 0.08 0.04 0.025 Table B2. Comparison of stochastic and deterministic predictions measure case deterministic ‰I Š ‰SI Š 1 284 17.1 ‰I Š ‰SI Š 2 ‰I Š ‰SI Š ‰I Š ‰SI Š stochastic error (%) case deterministic stochastic error (%) 0 0 n/a n/a 7 750 45.0 674 40.4 11 11 469 28.1 275 16.4 71 33 8 811 58.7 748 44.7 8.1 31 3 718 43.1 542 32.3 32 33 9 911 54.6 880 52.8 3.5 3.4 ‰I Š ‰SI Š 4 111 6.64 0 0 n/a n/a 10 920 34.5 865 32.3 6.4 6.8 5 342 25.6 172 12.7 99 101 11 736 55.2 663 49.5 11 12 ‰I Š ‰SI Š 6 627 25.1 378 14.9 66 68 12 1002 60.1 1036 61.8 3.3 2.8 There are 12 cases, the parameters of which are shown in table B1 (in units of day¡1 ). For each case, » ˆ 0:01, ¼ ˆ 0:005 and N ˆ 1500. Table B2 compares the equilibrium of the determinstic model and the long-time average of the stochastic model (the averaging starts after the system has settled down). There is only one stochastic run for each set of parameters, but the initial number of infectives is always large so the cases where the disease dies out are still relevant. In general, the agreement increases slightly with increasing concurrency and increasing endemicity. One can see that the agreement for cases 1^6 is much poorer than for cases 7^12. This is because the infection dynamics time-scale is closer to the partnership dynamics time-scale in cases 7^12. In our derivation for Q (I jSI) and Q (I jSS)c we assumed the independence of the state of one partner from the states of other partners of an individual. Clearly this assumption works better when the infection dynamics occur as quickly as, or not much more quickly than, the time-scale of partnership dynamics. In our analysis in ½ 4 we take this weakness of the moment closure into account in our parameter choices. From this observation we suggest that the pair approximation is still useful for dynamic networks, but must be applied more carefully as in comparison with static network models. REFERENCES Diekmann, O., Heesterbeek, J. & Metz, J. 1990 On the de¢nition and the calculation of the basic reproduction ratio R 0 in models for infectious disease in heterogeneous populations. J. Math. Biol. 28, 365^382. Diekmann, O., Dietz, K. & Heesterbeek, J. 1991 The basic reproduction ratio for sexually transmitted diseases. Theoretical considerations. Math. Biosci. 107, 325^339. Dietz, K. & Tudor, D. 1992 Triangles in heterosexual HIV transmission. In AIDS ep idemiology: methodological issues (ed. N. P. Jewell, K. Dietz & V. T. Farewell), pp. 143^155. Boston, MA: Birkha«user. Garnett, G. 1997 The natural history of syphilis: implications for the transmission dynamics and control of infection. SexuallyTransmitted Diseases 24, 185^200. Ghani, A. 1997 The role of sexual partnership networks in the epidemiology of gonorrhoea. Sexually Transmitted Diseases 24, 227^238. Kretzschmar, M. & Morris, M. 1996 Measures of concurrency in networks and the spread of infectious disease. Math. Biosci. 133, 165^195. Kretzschmar, M., Reinking, D. P., Van Zessen, G., Brouwers, H. & Jager, J. C. 1994 The basic reproduction ratio R 0 for a sexually transmitted disease in a pair formation model with two types of pairs. Math. Biosci. 124, 181^205. Morris, M. & Kretzschmar, M. 1997 Concurrent partnerships and the spread of HIV. AIDS 11, 1^7. Rand, D. 1999 Correlation equations and pair approximations for spatial ecologies. In Advanced ecological theory: p rinciples and applications (ed. J. McGlade), pp. 100^142. Oxford, UK: Blackwell Science. Watts, C. H. & May, R. M. 1992 The in£uence of concurrent partnerships on the dynamics of HIV/AIDS. Math. Biosci. 108, 89^104. Anderson, R. M. & May, R. M. 1991 Infectious diseases of humans. Oxford University Press. As this paper exceeds the maximum length normally permitted, the authors have agreed to contribute to production costs. Proc. R. Soc. Lond. B (2000)