The Role of Nonlinear Relapse on Contagion Amongst Drinking Communities
Ariel Cintrón-Arias, Fabio Sánchez, Xiaohong Wang, Carlos Castillo-Chavez, Dennis M. Gorman, Paul J. Gruenewald
TThe Role of Nonlinear Relapse on Contagion AmongstDrinking Communities
Ariel Cintr´on-Arias , Fabio S´anchez , Xiaohong Wang ,Carlos Castillo-Chavez , , , , Dennis M. Gorman , Paul J. Gruenewald Center for Research in Scientific ComputationBox 8205, North Carolina State University, Raleigh, NC 27695-8205 Department of Biological Statistics and Computational BiologyCornell University, Ithaca, NY 14853-7801 Mathematical, Computational, and Modeling Sciences CenterP.O. Box 871904, Arizona State University, Tempe, AZ 85287 School of Human Evolution and Social ChangeArizona State University, Tempe, AZ 85287 School of Mathematics & StatisticsArizona State University, Tempe, AZ 85287 Santa Fe Institute1399 Hyde Park Road, Santa Fe, NM, 87501 Department of Epidemiology & Biostatistics, School of Rural Public HealthTexas A&M Health Science Center, P.O. Box 1266, College Station, TX 77843-1266 Prevention Research Center1995 University Avenue, Suite 450, Berkeley, CA 94704
December 4, 2008 a r X i v : . [ q - b i o . P E ] A p r bstract Relapse, the recurrence of a disorder following a symptomatic remission, is a fre-quent outcome in substance abuse disorders. Some of our prior results suggested thatrelapse, in the context of abusive drinking, is likely an “unbeatable” force as long asrecovered individuals continue to interact in the environments that lead to and/or re-inforce the persistence of abusive drinking behaviors. Our earlier results were obtainedvia a deterministic model that ignored differences between individuals, that is, in arather simple “social” setting. In this paper, we address the role of relapse on drink-ing dynamics but use models that incorporate the role of “chance”, or a high degreeof “social” heterogeneity, or both. Our focus is primarily on situations where relapserates are high. We first use a Markov chain model to simulate the effect of relapse ondrinking dynamics. These simulations reinforce the conclusions obtained before, withthe usual caveats that arise when the outcomes of deterministic and stochastic modelsare compared. However, the simulation results generated from stochastic realizationsof an “equivalent” drinking process in populations “living” in small world networks,parameterized via a disorder parameter p , show that there is no social structure withinthis family capable of reducing the impact of high relapse rates on drinking preva-lence, even if we drastically limit the interactions between individuals ( p ≈ Keywords: drinking behavior; deterministic model; stochastic model; small-world network;social influence; drinking dynamics.
The mechanisms responsible for observed drinking patterns within and between populationsare complex (Daido 2004; Weitzman et al. 2003; Mubayi et al. 2008; and references therein).The development of compartmental and mathematical frameworks geared towards the iden-tification of key “transition” mechanisms that increase the percentage of abusive drinkersmust factor in the impact of individuals’ socioeconomic characteristics, their propensity todrink (heavy drinking tends to run in families), changes in local environments (going to col-lege), treatment failure, ineffectiveness of educational efforts, cultural norms and communityvalues (Mubayi et al. 2008; and references therein).The term drinking (population) dynamics refers to the study and identification of “aver-age” mechanisms, at the individual level, responsible for observed drinking patterns withinthe organizational and temporal scales of interest. We model drinking dynamics at the2opulation level as the result of individuals’ social contacts in pre-specified environments(“drinking contagion”). This modeling approach has proved useful in the identification ofthe mechanisms behind social patterns that are thought to be, in part, an outcome of intenseinteractions between individuals in shared social environments. This modeling approach hasbeen applied to the study of the spread of scientific ideas and innovations (Bettencourt et al.2006 ); in studies that focus on the mechanisms behind the observed increases in prevalenceof eating disorders (Gonz´alez et al. 2003); in studies that address the impact of relapse onthe distribution of drinkers (S´anchez et al. 2007; S´anchez 2006); in studies that envisionviolence as an epidemic (Patten and Arboleda-Florez 2004); as explanation for the observedgrowth or decline of crime in cities (Gladwell 1996); and in studies that highlight the ex-plosive increases in the use of illicit drugs, such as ecstasy (Song et al. 2006; Mackintoshand Stewart 1979). Researchers are interested in studying the impact of individual drink-ing habits and preferences’ variability at multiple levels of social organization: from small“isolated” to highly connected communities; and over short or long time horizons. Modelshave been used to explore the capacity of drinking environments to support communitiesof drinkers as well as the impact of individuals’ movements between drinking venues on theoverall distribution of drinking types (Mubayi et al. 2008).The National Institute on Alcohol Abuse and Alcoholism estimates that 18 million Amer-icans suffer from alcohol abuse or dependence. Alcohol-related problems cost the UnitedStates (U.S.) nearly $185 billion annually (National Institute of Alcohol Abuse and Al-coholism 2008a) while alcohol abuse was responsible for nearly 80,000 fatalities per yearduring 2001-05, and it is now the third leading cause of death in the U.S. (Centers for Dis-ease Control and Prevention 2008a). Prevention and control efforts that include treatmentand education programs that target specific populations including children (Leadership toKeep Children Alcohol Free 2008) or adolescents (College Drinking 2008) are in need of im-provement. Among the many problems confronting these programs are the very high ratesof relapse after treatment that are observed. Up to 70% of treated alcohol abusers relapseafter treatment (reviewed in S´anchez et al. 2007). Mathematical studies can be particularlyeffective as guides to the evaluation, testing and implementation of single or multiple inter-vention strategies over short or long time scales. This is particularly true in the study ofchronic relapsing diseases such as alcohol addiction.
Social dynamics, disease transmission, and social structure
Several aspects linked to disease transmission depend strongly on a population’s social dy-namics. Disease dynamics can often be driven by factors that include heterogeneity inbehavior, frequency of use of mass transportation, travel patterns, and cultural norms andpractices. Examples where the use of mathematical models have generated useful insightsinclude studies on the role of behavior on the transmission dynamics of sexually transmitteddiseases like gonorrhea or HIV (Castillo-Chavez et al. 2003; Hethcote 2000; Hethcote andYorke 1984; Anderson and May 1991; and references therein) and studies on the intensityand frequency of travel on the spread of communicable diseases such as SARS (Chowellet al. 2003; Song et al. 2003) and influenza (Hyman et al. 2003; Chowell et al. 2006a).The most significant study of the role of heterogenous mixing on the transmission dynamicsof gonorrhea was carried out by Hethcote and Yorke (Hethcote and Yorke 1984). These3esearchers through their introduction of the concept of core group (outliers in the distribu-tion of sexually-active individuals) showed that most secondary cases of gonorrhea infectionscould be traced to the core (most connected nodes in a network of sexually-active individu-als). Furthermore, they showed that focusing surveillance and treatment on core subpopu-lations resulted in significant reductions in gonorrhea prevalence. The public health policyat that time wrongly focused on the “random” testing of women, a policy derived from datathat showed that a large percentage of gonorrhea infected women are indeed asymptomatic(Hethcote and Yorke 1984; and references therein).The systematic study of the role of heterogenous social landscapes on disease dynamics beganin direct response to efforts to stop the HIV epidemics. Efforts to compute explicit mixingmatrices (who had interactions with whom) and to study the impact of sexual preference inthe context of HIV transmission intensified (Blythe et al. 1990; Blythe et al. 1995; Blythe etal. 1991; Busenberg et al. 1989; Busenberg et al. 1991; Hsu 1993; Hsu et al. 1994; Hsu et al.1996; Castillo-Chavez et al. 1996; Hethcote 2000; Anderson and May 1991; Castillo-Chavez1989; and references therein).Most recently, efforts to explore disease dynamics in the context of heterogenous (fixed)social network structures have proved quite fruitful. The study of epidemics on networkhas increased our understanding of the role of “social” heterogeneity on disease dynamics(Newman 2003; and references therein) but the impact of the efforts of the mathematical“network” community goes beyond the study of epidemics on networks, as is evident fromthe wealth of applications found in the literature (see Watts and Strogatz 1998; Barabasiand Albert 1999; Newman 2003; Newman et al. 2006; and references therein). There is abody of research that contributes to the characterization and validation of some classes ofnetwork structures with data (Meyers et al. 2005); structures whose statistical propertiesare most often captured via power law distributions (Newman et al. 2006). The class of bestknown or more popular models of this type include small-world (Watts and Strogatz 1998)and scale-free (Barabasi and Albert 1999) networks.Social network analysis is the result (to a great degree) of major contributions by socialscientists (Wasserman and Faust 1994; and references therein). Recent contributions bymathematical scientists (Newman 2003 and references therein; Newman et al. 2006; Wattsand Strogatz 1998) have increased interactions between social and mathematical scientists.Applications that make use of specialized network structures include studies of the structureof scientific co-authorship networks (Newman 2003), the organizational structure of com-mittees in the U.S. House of representatives (Porter et al. 2005), the structure of internetnetworks (Pastor-Satorras and Vespignani 2001), the properties of contact tracing networksfor SARS (Meyers et al. 2005), and the nature of sexual partnership networks (Liljeros etal. 2001). Efforts to study stochastic epidemic and social processes on networks have alsobeen carried out in the context of homeland security (Chowell and Castillo-Chavez, 2003and references therein) and drinking (Braun et al. 2006). Our goal here is “theoretical”,that is, we focus on the study of drinking on some networks characterized by scaling laws(Newman 2003; and references therein). Specifically, the primary objective is to explore therole of network structure on the distribution of drinkers in communities (small world type)4here relapse rates are high.This manuscript is organized as follows. Section 2 revisits the results in (S´anchez et al.2007; S´anchez 2006) on the role of relapse on the distribution of drinking types. Section 3introduces the stochastic analog of the deterministic model to highlight the role of variabilityin the distribution of drinking types of Section 2. Section 4 simulates one version of thestochastic drinking dynamics in a small-world network. Finally, Section 5 discusses the roleof relapse in these settings.
In the drinking model formulation proposed in (S´anchez et al. 2007), the population is di-vided in three classes: S ( t ), moderate and occasional drinkers (Centers for Disease Controland Prevention 2008c), D ( t ), problem or heavy drinkers (Centers for Disease Control andPrevention 2008d; National Institute of Alcohol Abuse and Alcoholism 2008b), and tem-porarily recovered, R ( t ). Table 1 presents the definitions used in (S´anchez et al. 2007)where it is assumed that the population is composed of “average” individuals that interactat random with each other. The proportion of contacts of S -individuals with D -individualsper unit of time is therefore proportional to D/N where N = S + D + R , denotes the totalsize of the community. The progression rate from S to D and the relapse rate from R to D depend on frequency-dependent (random) interactions.In (S´anchez et al. 2007) the model is given by the following set of nonlinear differentialequations: dSdt = µN − βS ( t ) D ( t ) N − µS ( t ) , (1) dDdt = βS ( t ) D ( t ) N + ρR ( t ) D ( t ) N − ( µ + φ ) D ( t ) , (2) dRdt = φD ( t ) − ρR ( t ) D ( t ) N − µR ( t ) , (3) N = S ( t ) + D ( t ) + R ( t ) , (4)where β denotes the per-capita effective contact rate (transmission rate), that is, βSD/N denotes the rate of transitions from S to D , the result of the frequency-dependent interactionsbetween individuals in the classes S and D ; µ denotes the per-capita departure rate fromthe system; ρ denotes the per-capita effective relapse rate, that is, ρRD/N denotes the rateof transitions from R to D , the result of the frequency-dependent interactions between R and D ; φ denotes the per-capita recovery (treatment or education) rate; and µN denotesthe total recruitment rate into this homogeneous social mixing community. It is assumedthat all “recruits” are S -individuals. Hence, we set the S -recruitment rate equal to µN as it5uarantees constant population size. The validity of the analysis is therefore tied to a timehorizon where changes in total population size are minimal.The reproductive number under a treatment/education regime φ is given by R φ ≡ R ( φ ) = βµ + φ . (5) R φ is a dimensionless quantity (ratio or number) that can be interpreted as the number of D -individuals “generated” in a population of primarily S -individuals sharing a common envi-ronment. That is, if we start with S ≈ N individuals and introduce a “typical” D -individualthen we expect R φ secondary cases generated from the S population per D -individual, butonly at the start of the “outbreak”. Hence, R φ > D -community if N is large enough. We also expect that when R φ <
1, the introductionof D -individuals in a population where S ≈ N ( N large) will not result in the growth and(eventual) establishment of a problem-drinking community ( D -individuals). The above ob-servations are on target when the rate of relapse is linear, that is, ρR rather than ρRD/N .However, when the relapse rate is nonlinear, namely, ρRD/N , the outcome is not as “ex-pected”. The outcome depends on the ratios R ρ = ρβ [1 − R ( φ )] (6) R c = ρβ (cid:34)
11 + R − (cid:114) R − µρ (cid:35) , (7)where R ( φ ) is defined in Equation (5); R ≡ R (0) = β/µ . R ρ can be interpreted as the number of problem drinkers ( D -individuals) generated fromthe R -class as a result of the frequency-dependent interactions between the R - and D -classes( R -individuals remain in the same environment). We observe that R ρ > R ( φ ) <
1. On the other hand R c > βµ + β > (cid:114) R − µρ > . We have not been able to interpret the meaning of R c in social terms. However, the valueof R c , under some conditions, provides a sharp D -extinction threshold, that is, a thresholdthat if crossed, would lead to the eventual elimination of the D -class, independent of initialconditions ( D (0)).The distribution of drinking types, in the nonlinear relapse rate case, depends not only onthe thresholds R φ , R ρ , and R but also on the size of the initial population of problemdrinkers, D (0). In (S´anchez et al. 2007) the following results were obtained:1. If R ( φ ) > D -class becomes established.6. Whenever R c < R ( φ ) < R ρ < R ( φ ) < R c < D -classbecomes (eventually) extinct.3. Whenever R c < R ( φ ) < R ρ > D -class becomes establishedis a function of the initial size of the class of D -individuals, D (0) (see Figure 1(c), (d)).Numerical simulations —Figure 1(a), (c), (d)— illustrate the role of initial conditions ondrinking dynamics. Nonlinear relapse leads to a system that supports two socially acceptablecoexisting stable equilibria ( D ≡ D > R ( φ ) (with R ρ > ρ > φ , leads to explosivegrowth in the D -class as long as D (0) (the initial population of problem drinkers) is “largeenough” (see Figure 1(a)). The qualitative behavior displayed in Figure 1(a) is commonlycalled a “backward” bifurcation (S´anchez et al. 2007). We further observe that once thepopulation of problem drinkers becomes established ( R c < R ( φ ) <
1) their extinction canonly be carried out if φ increases to the point where R ( φ ) < R c or if ρ decreases to thepoint where R ρ <
1. Figures 1(c), (d), display D ( t ) versus t to illustrate, with a timeseries, the effects of initial conditions, D (0). We observe bistability. The size of the initialnumber of problem drinkers determines whether or not a D -community becomes establishedeven under unfavorable conditions ( R ( φ ) < ρ = φ , we observe (Figure 1(b)) that the D -class grows (gradually) with R ( φ ); multiple endemic (non-negative) stable D -equilibria will not co-exist in this case.When ρ = φ , R ( φ ) < The stochastic model of this section is built from the deterministic model given by System(1)–(4) and is used to quantify the role of variability on drinking dynamics. Here, weconcentrate on an stochastic analog to the “mean field” model given by Equations (1)–(4),the deterministic model that supports two positive equilibria ( R c < R φ < R ρ > R φ > D -class before the preselected time horizon.Setting R φ < R c < D -class in thedeterministic formulation but not always (as expected) in the stochastic formulation (Allenand van den Driessche 2006; Allen 2003).In well-established drinking communities (including college students) estimates clearly showthat R φ >
1. Thus, one may ask whether the existence of backward bifurcations (bi-stability)is just of theoretical value? If the goal is to prevent the formation of a drinking communitythen the above question “makes” sense. However, most often the goal is to reduce or eliminatethe D -class and the existence of a backward bifurcation makes this much harder.Relapse rates among problem drinkers are high (Miller et al. 2001; Daido 2004). Hence,the existence of a relapse driven backward bifurcation suggests that efforts to “eliminate”problem drinkers or reduce problem drinking may be futile as long as “ R -individuals” remainin the same social environment. Substantial reductions in the relapse parameter—withthe ultimate goal of having R φ < D -class (see bifurcation diagram in Figure 1(a)).Histograms (based on 50 stochastic realizations) of the number of problem drinkers at astoppage time T , denoted by D ( T ), are examined when R φ > R φ < R φ > D ( T ) lies in[350 , D -class (less than 10%) when R φ < ρ > φ ). Inother words, it is possible for a population of problem drinkers to become established evenif R φ < A network (graph) is a set of nodes with connections (edges) between them. Graphs providevisual representations of the contact structure of individuals in a population (Newman 2003).The fact that all social processes (including drinking) depend on contacts between distinctindividuals has, in part, motivated the study of epidemics on networks (May and Lloyd 2001;Meyers et al. 2005; Pastor-Satorras and Vespignani 2001; Grabowski and Kosinski 2005).Watts and Strogatz (Watts and Strogatz 1998) introduced a one-parameter, p , family of8etworks. As the disorder parameter p is varied in [0,1], the graph moves from a regularlattice to a random graph. The model can be formulated algorithmically as follows: theinitial network is initialized via a one-dimensional periodic ring lattice of N nodes, eachconnected to its closest (cid:104) k (cid:105) neighbors (two nodes are neighbors if there is an edge connectingthem). The network is updated by re-wiring each edge with probability p (the disorderparameter) to a randomly selected node until it reaches “fixed” statistical properties. When p → p →
1, most edges are rewired, theresulting network is a random graph (Bollobas 2001). Watts and Strogratz showed that theuse of just a few random long-range connections ( p small) drastically reduced the average distance between any pair of nodes (Watts and Strogatz 1998) —the kind of property thatenhances “transmission”, the “small-world effect”. The effect was postulated based on theresult of a series of letter-forwarding experiments carried out by S. Milgram (Milgram 1967).The statistical properties of small-world and “similar” networks have been studied (Wattsand Strogatz 1998; Newman et al. 2006; and references therein).Here we model community structure as a small-world network. The terms network andcommunity are used interchangeably, with nodes representing individuals and edges denotingthe social connections or interactions, the kind of “social mixing” that may lead to node“transition” (from the moderate drinker into the problem drinker state). Nodes can be inone of three distinct states: moderate drinker, problem drinker, and recovered drinker. Thestochastic transitions between nodes’ states are modeled as functions of time and the numberof “neighbors” in particular states (transition rates). If one starts with a community with N nodes where Node i (1 ≤ i ≤ N ) has δ ( i, t ) neighbors who, at time t , are in the state “problemdrinker”, then the probabilities that Node i changes its state given that it alters its state, ateach time step are: from moderate to problem drinker, 1 − exp( − βδ ( i, t )); from problem torecovered, 1 − exp( − φ ); and from recovered to problem drinker, 1 − exp( − ρ τ ( t ) δ ( i, t )). Thisformulation (see Table 3) defines a stochastic process on the random variables S p ( t ), D p ( t ),and R p ( t ). These random variables can also be thought of as parameterized by the disorderparameter p ∈ [0 , D p ( T ) and R p ( T ), where T denotesthe stoppage time in the simulations (see Table 4), are computed for each value of p (seeFigure 4). Figures 5 and 6 highlight the mean and variance (over 20 realizations) of D p ( T )and R p ( T ) as a function of p (Chowell and Castillo-Chavez 2003; Chowell et al. 2006b).A drinking wave is detected even as the size of the problem drinking class goes to zero forthe case ρ = 0 (no relapse) with R φ >
1. This feature agrees with deterministic (Brauer andCastillo-Chavez 2001) and stochastic “theories” (Allen 2003) on single-outbreak SIR models.Figure 5(a) shows that variations on the network structure (modeled by p ) have no effect onthe mean size of the problem drinker class D p ( T ) . However, the mean size of the recoveredclass R p ( T ) exhibits a phase transition as p → − (Figure 5(b)). Hence, in the absence ofvital dynamics (births and deaths) and relapse, we conclude that community structure does9ffect the average size of the problem drinking class during the drinking wave. Small valuesof “ p ” lead to a phase transition (Newman 2003), a “small world” effect.Figure 6 illustrates a worst case scenario in which the average relapse probability is nearone for the majority of the time. To see the impact of high, nearly stationary relapse rates,we let (cid:104) k (cid:105) denote the average number of connections per node in a one-dimensional latticewhen p = 0 and carry out simulations on this network with the average relapse probability(1 − e − ρ τ ( t ) (cid:104) k (cid:105) ) ≈
1. The relapse rate ρ τ ( t ) (defined in Table 4) is modeled as a stepwiseconstant function that drops its value at precisely t = τ . The worst case scenario herecorresponds to the case where τ = ∞ . In general, when relapse rates are high for too long,small-world structures (any value of p ) have no effect on the mean sizes of the problem andrecovered drinking classes. In fact, the size of the problem drinking community is above60% regardless of the value of p (other parameters kept fixed). Furthermore, we see thaton average D p ( T ) + R p ( T ) = N when relapse rates are high. That is, every member of thisclosed population becomes a problem drinker at least once regardless of the value of p .Reducing the relapse rate from 0.90 to 0.12 at precisely the time τ reduces the average relapseprobability from 1 − e − . (cid:104) k (cid:105) ≈ .
00 to 1 − e − . (cid:104) k (cid:105) ≈ .
50 at time τ . Figure 7 shows theimpact of increasing the values of τ = 3 , , ,
10. We do not observe a lot of differences inthe average values of D p ( T ) and R p ( T ) as a function of τ . However, these averages “improve”in the “right” direction as τ reduces its value from τ = ∞ towards τ = 0. Relapse has a significant impact on the dynamics of addictive behavior (Gonzalez et al. 2003;S´anchez et al. 2007; Song et al. 2006; and references therein). The use of a simple systemof differential equations (S´anchez et al. 2007) shows that for socially-intense processes likedrinking, the reproductive number, R φ is not always the key. Frequency dependent relapserates play a huge role. Frequency dependent relapse rates do increase the possibility of severeoutbreaks within “well-behaved” communities, but more importantly they also increase thelikelihood of failure of programs aimed at eliminating drinking. S´anchez et al. (S´anchez et al.2007) clearly delineated the possibilities from their mathematical analysis of a simple modelwhere all the mixing takes place in the same drinking environment. Mubayi et al. (Mubayiet al. 2008) recently explore the impact of individuals’ movement between heterogeneousdrinking environments. They showed that frequent movement between distinct environmentscan have a significant (negative) effect on the distribution of drinking types. Here, we onlyfocused on exploring the predictions of (S´anchez et al. 2007) in two stochastic settings. Thestochastic analog (continuous time Markov chain) of S´anchez et al.’s deterministic model wasused to highlight the role of variability. The results were consistent with those of Sanchezet al. with the usual caveats (Allen 2003). A small-world network was used to highlight thevery strong role played by relapse. 10n fact, our study of drinking in a small-world network parameterized by the disorder param-eter p leads to the following results: When there is no relapse ( ρ = 0), we recovered the wellunderstood phase transition effect previously identified from SIR simulations on small-worldnetworks (Newman 2003), as p crosses a critical value; the introduction of high relapse rates“eliminates” the role of “ p ”. In other words, the form of social connections (who interactswith whom) in populations experiencing strong patterns of relapse has no impact on theprevalence of addictive behaviors. Hence, if relapse rates are high then emphasis on pro-grams that generate substantial and sustained reductions in “mixing” will not be effective.Reducing residence times in risky environments which promote relapse, reducing recruitmentinto drinking communities and reducing movement between drinking venues are more likelyto be effective (Mubayi et al. 2008). Transitions between drinking classes involve discrete events which change the number ofindividuals in every class, one at a time. For example, when a drinking “contagion” eventoccurs, the number of moderate drinkers is decreased by one, while the number of problemdrinkers increases by one. The probability that an event takes place during an infinitesimaltime interval [ t, t + dt ] is calculated from the average rates in the deterministic model. Inthis example, the “conversion” event occurs at the rate of βS ( t ) D ( t ) /N and the probabilitythat it happens in [ t, t + dt ] is approximately ( βS ( t ) D ( t ) /N ) dt . All the events, their ratesof occurrence, and the probabilities at which they take place are listed in Table 2.It is assumed that the events are described by independent Poisson processes (Allen, 2003).The term E = µN + µS + µD + µR + βSD/N + φD + ρRD/N, denotes the rate at which an event occurs at time t . The time between events is exponentiallydistributed with mean 1 /E . The time at which the next event happens is found, for eachrealization, by sampling from an exponential distribution with mean 1 /E .To decide which event takes place (once it is known that an event occurs), we divide upthe interval (0 , E ) into subintervals that correspond to the relative occurrence probabilitiesof the various events. For example, given that an event has occurred, the probability thatit is a recruitment is µN/E , the probability of the removal of a moderate drinker is µS/E ,the probability of the removal of a problem drinker is µD/E , etc. A number U is selectedrandomly from the uniform distribution on (0 ,
1) and an event is selected if this value fallswithin the appropriate subinterval. For instance, the event is a recruitment if U satisfies0 < U < µN/E , a moderate drinker removal if U lies between µN/E and ( µN + µS ) /E , aproblem drinker removal if U lies between ( µN + µS ) /E and ( µN + µS + µD ) /E , and so on.11 cknowledgments A.C.-A. was supported in part by the Statistical and Applied Mathematical Sciences In-stitute which is funded by the National Science Foundation under Agreement No. DMS-0112069. Any opinions, findings, and conclusions or recommendations expressed in thismaterial are those of the authors and do not necessarily reflect the views of the NationalScience Foundations. A.C.-A. was also supported in part by Grant Number R01AI071915-07 from the National Institute of Allergy and Infectious Diseases. The content is solelythe responsibility of the authors and does not necessarily represent the official views of theNIAID or the NIH. X.W. and C.C.-C. were supported by NIAAA grant on “Ecosystem Mod-els of Alcohol-Related Behavior”, Contract No. HHSN2S1200410012C, ADM Contract No.No1AA410012 through Prevention Research Center, PIRE, Berkeley, the National ScienceFoundation (DMS-0502349), the National Security Agency (dod-h982300710096), the SloanFoundation, and Arizona State University. P.J.G. was supported by Grant Number R01AA06282 from the National Institute on Alcohol Abuse and Alcoholism.
References [1] Allen,L. J. S., van den Driessche, P.: Stochastic epidemic models with a backwardbifurcation. Math. Biosci. Eng. , 445-458 (2006).[2] Allen, L. J. S.: An Introduction to Stochastic Processes with Applications to Biology.Pearson, Upper Saddle River, (2003).[3] Anderson,R., May, R.: Infectious Diseases of Humans: Dynamics and Control. OxfordUniversity Press, Oxford, (1991).[4] Barabasi, A.L., Albert, R.: Emergence of scaling in random networks. Science (5439): 509-512 (1999).[5] Bettencourt, L. M. A., Cintr´on-Arias, A., Kaiser, D.I., Castillo-Ch´avez, C.: The powerof a good idea: quantitative modeling of the spread of ideas from epidemiological models.Physica A , 513-536 (2006).[6] Blythe, S., Castillo-Chavez, C.: Scaling law of sexual activity, Nature, , 202 (1990).[7] Blythe, S., Castillo-Chavez, C., Palmer, P.,Cheng, M.: Towards a unified theory ofmixing and pair formation. Math. Biosci. , 379-405 (1991).[8] Blythe, S., Busenberg, S., Castillo-Chavez, C.: Affinity and paired-event probability.Math. Biosci. , 265-284 (1995).[9] Bollob´as, B.: Random Graphs. Cambridge University Press, Cambridge, (2001).1210] Brauer, F., Castillo-Ch´avez, C.: Mathematical Models in Population Biology and Epi-demiology. Springer-Verlag, New York, (2001).[11] Braun, R. J., Wilson, R. A., Pelesko, J. A., Buchanan, J. R.: Applications of small-worldnetwork theory in alcohol epidemiology. J. Stud. Alcohol , 591-599 (2006).[12] Busenberg, S., Castillo-Chavez, C.: Interaction, pair formation and force of infectionterms in sexually transmitted diseases. In: Castillo-Chavez, C. (ed.) Mathematical andStatistical Approaches to AIDS Epidemiology. Lecture Notes Biomathematics, Vol. 83,pp. 280-300. Springer-Verlag, Berlin, (1989).[13] Busenberg, S., Castillo-Chavez, C.: A general solution of the problem of mixing ofsubpopulation, and its application to risk- and age-structured epidemic models. IMA J.Math. Appl. Med. Biol. , 1-29 (1991).[14] Castillo-Ch´avez, C. (ed.): Mathematical and Statistical Approaches to AIDS Epidemi-ology. Lecture Notes in Biomathematics, Vol. 83. Springer-Verlag, Berlin, (1989).[15] Castillo-Chavez, C., Huang, W., Li, J.: Competitive exclusion in gonorrhea models andother sexually-transmitted diseases. SIAM J. Appl. Math. , 1-8 (2003).[22] Chowell, G., Castillo-Ch´avez, C.: Worst-case scenarios and epidemics. In: Banks, H.T.,Castillo-Chavez, C. (eds.) Bioterrorism: Mathematical Modeling Applications in Home-land Security. Frontiers in Applied Mathematics, Vol. 28, pp. 35-53. Society for Indus-trial and Applied Mathematics, Philadelphia, (2003).1323] Chowell, G., Ammon, C. E., Hengartner, N. W., Hyman, J. M.: Transmission dynamicsof the great influenza pandemic of 1918 in Geneva, Switzerland: Assessing the effectsof hypothetical interventions. J. Theor. Biol. , 383-386 (2004).[27] Gladwell, M.: The tipping point. New Yorker , 32-39 (1996).[28] Gonz´alez, B., Huerta-S´anchez, E., Ortiz-Nieves, A., V´azquez-Alvarez, T., Kribs-Zaleta,C.: Am I too fat? Bulimia as an epidemic. J. Math. Psychol. , 515-526 (2003).[29] Grabowski, A., Kosinski, R. A.: The SIS model of epidemic spreading in a hierarchicalsocial network. Acta Phys. Pol. B , 1579-1593 (2005).[30] Hethcote, H., Yorke, J.: Gonorrhea Transmission Dynamics and Control. Lecture Notesin Biomathematics, Vol. 56. Springer-Verlag, Berlin (1984).[31] Hethcote, H.: The mathematics of infectious diseases. SIAM Rev. , 907-908 (2001).[38] Mackintosh, D. R., Stewart, G.T.: A mathematical model of a heroin epidemic: impli-cations for control policies. J. Epidemiol. Commun. H. , 299-304 (1979).[39] May, R. M., Lloyd, A. L.: Infection dynamics on scale-free networks. Phys. Rev. E. ,066112 (2001).[40] Meyers, L.A., Pourbohloul, B., Newman, M.E.J., Skowronski, D.M., Brunham, R.C.:Network theory and SARS: Predicting outbreak diversity. J. Theor. Biol. , 71-81(2005).[41] Miller, W. R., Walters, S. T., Bennett, M. E.: How effective is alcoholism treatment inthe United States? J. Stud. Alcohol. , 211-220 (2001).[42] Milgram, S.: The small world problem. Psychol. Today ,167-256 (2003).[47] Newman, M.E.J., Barabasi, A.L., Watts, D.J.: The Structure and Dynamics of Net-works. Princeton University Press, Princeton, (2006).[48] Orford, J., Krishnan, M., Balaam, M., Everitt, M., Van der Graaf, K.: Universitystudent drinking: the role of motivational and social factors. Drug-Educ. Prev. Polic. , 407-421 (2004).[49] Pastor-Satorras, R., Vespignani, A.: Epidemic spreading in scale-free networks. Phys.Rev. Lett. , 3200 (2001).[50] Patten, S. B., Arboleda-Florez, J. A.: Epidemic theory and group violence. Soc. Psych.Psych. Epid. , 853-856 (2004).[51] Porter, M. A., Mucha, P. J., Newman, M. E. J., Warmbrand, C. M.: A network analysisof committees in the U.S. house of representatives. P. Natl. Acad. Sci. USA, , 7057-7062 (2005). 1552] E. Renshaw, Modelling Biological Populations in Space and Time. Cambridge UniversityPress, Cambridge, (1991).[53] F. S´anchez: Studies in Epidemiology and Social Dynamics. Dissertation, Cornell Uni-versity (2006).[54] S´anchez, F., Wang, X., Castillo-Chavez, C., Gorman, D. M., Gruenewald, P. J.:Drinking as an epidemic—a simple mathematical model with recovery and relapse.In: Witkiewitz, K. A., Marlatt, G. A. (eds.) Therapist’s Guide to Evidence-Based Re-lapse Prevention: Practical Resources for the Mental Health Professional. Academic,Burlington, 353-368 (2007).[55] Song, B., Castillo-Garsow, M., Castillo-Ch´avez, C., R´ıos Soto, K., Mejran, M., Henso,L.: Raves, Clubs, and Ecstasy: The Impact of Peer Pressure. Math. Biosci. Eng. ,249-266 (2006).[56] Wasserman, S., Faust, K.: Social network analysis: methods and applications. Cam-bridge University Press, Cambridge, (1994).[57] Watts, D. J., Strogatz, S. H.: Collective dynamics of ‘small-world’ networks. Nature , 440-442 (1998).[58] Weitzman, E. R., Folkman, A., Folkman, K. L., Weschler, H.: The relationship ofalcohol outlet density to heavy and frequent drinking and drinking-related problemsamong college students at eight universities. Health Place , 1-6 (2003).16 Tables
Table 1: State variables and parameters of the contagion model in (S´anchez et al. 2007).State variable Description S ( t ) Number of occasional and moderate drinkers at time tD ( t ) Number of problem drinkers at time tR ( t ) Number of recovered individuals at time t Parameter Description β Effective transmission rate (average number of effective interactionsper occasional and problem drinker per unit of time) ρ Community-driven relapse rate (average number of effective interactionsper problem drinker and recovered individual per unit of time) φ Per-person treatment rate µ Per-person departure rate from the drinking environment N Community size (permanent population size)Table 2: Collects the transition rates and infinitestimal probabilities of occurrence of theevents linked to a single drinking model outbreak. The dependence on t is omitted, writing S , D , and R , instead of S ( t ), D ( t ), and R ( t ), respectively.Event Transition Rate at which Probability of transitionevent occurs in time interval [ t, t + dt ]Recruitment S → S + 1 µN µN dt Moderate drinker removal S → S − µS µSdt Problem drinker removal D → D − µD µDdt Sober removal R → R − µR µRdt Drinking contagion S → S − D → D + 1 βS DN βS DN dt Recovery D → D − R → R + 1 φD φDdt Relapse D → D + 1, R → R − ρR DN ρR DN dt δ ( i, t ) Number of problem drinker neighbors of node i at time tS p ( t ) Total number of moderate drinkers at time t in asmall-world community parameterized by pD p ( t ) Total number of problem drinkers at time t in asmall-world community parameterized by pR p ( t ) Total number of recovered individuals at time t in asmall-world community parameterized by p Parameter Description β Transmission rate φ Per-person treatment rate ρ τ ( t ) Time-dependent relapse rateEvent Probability of transitionNode i changes from moderate into problem drinker 1 − e − βδ ( i,t ) Node i switches from problem drinker into recovered − e − φ Node i changes from recovered into problem drinker − e − ρ τ ( t ) δ ( i,t ) Table 4: Parameter values utilized in simulations of drinking dynamics in small-world com-munities.Parameter Description Baseline value (cid:104) k (cid:105) Average connectivity per node 6 N Community size 1000 β Transmission rate 0.12 φ Per-person treatment rate 0.7 ρ τ ( t ) Time-dependent relapse rate ρ τ ( t ) = 0 .
90 whenever t < τρ τ ( t ) = 0 .
12 if t ≥ τT Stoppage time 4000 D p (0) Initial number of problem drinkers chosenuniformly at random in every community 5Number of stochastic realizations 2018 Figures ! P r ob l e m d r i n ke r s a t e qu ili b i r i u m Bifurcation diagram (a)
0 0.3 0.6 0.9 1.2 1.505001000150020002500 R eproductive number R ! P r ob l e m d r i n ke r s a t e qu ili b i r i u m Bifurcation diagram (b) ! = 0.90T im e t P r ob l e m d r i n ke r s D (t) (c) ! = 0.90T im e t P r ob l e m d r i n ke r s D (t) (d) Figure 1: Numerical simulations of drinking model in a homogeneous drinking community.Panel (a) shows a bifurcation diagram that involves the number of problem drinkers atequilibrium versus the reproductive number R φ , when φ < ρ . Panel (b) displays a bifurcationdiagram illustrating the special case when the recovery rate equals the relapse rate ( φ = ρ =0 . R φ < D ( t ) versus t under differentinitial conditions. In Panel (c) the initial conditions are S (0) = 0 . N , D (0) = 0 . N and R (0) = 0; in Panel (d) they are S (0) = 0 . N , D (0) = 0 . N and R (0) = 0. The parametervalues used are: N = 10000, µ = 0 . φ = 0 .
50 and ρ = 7 .
00, 0 . ≤ β ≤ .
50 (Panel (a)); N = 10000, µ = 0 . φ = ρ = 0 .
50, 0 . ≤ β ≤ .
50 (Panel (b)); N = 10000, µ = 0 . φ = 0 .
50 and ρ = 7 . β = 0 .
90 (Panels (c) and (d)).19 N u m b e r o f p r ob l e m d r i n ke r s Figure 2: Results from numerical simulations. 50 stochastic realizations (grey curves) andnumerical solutions of the deterministic (black curve) problem drinker class D ( t ) versus time t . For these simulations the following values of parameters were used: N = 1000, β = 1 . ρ = 7 . φ = 0 .
50 and µ = 0 .
50 with R φ = 1 .
20 and the initial number of problem drinkers D (0) = 5. φ =1.20 (a) φ =0.90 (b) Figure 3: Histograms of D ( T ), number of problem drinkers at stoppage time T = 50000,resulting from 50 stochastic realizations with R φ > R φ <
00 625 650 675 7001 2 3 D p (T) 300 325 350 375 4001 2 3 R p (T) (a ) (b) Figure 4: Histograms of the total number of problem drinkers and recovered individuals, D p ( T ) and R p ( T ), respectively, at a stoppage time T . Samples obtained from 20 stochasticrealizations in simulated communities with p = 3 . × − in community size 1000 (nodes). − − − − − − M ea n D p ( T ) − − − − M ea n R p ( T ) Network disorder parameter p (b)(a )
Figure 5: Average and variance of D p ( T ) and R p ( T ) as functions of the simulated communityarchitecture parameterized by p (logarithmic scale). The mean (circles) and mean plus andminus one standard deviation (dash curves) are computed from 20 stochastic realizations foreach fixed value of p . Panels (a) and (b) display results of simulated contagion in small-worldcommunities in the absence of relapse, ρ ≡ − − − − M ea n D p ( T ) − − − − M ea n R p ( T ) Network disorder parameter p (b)(a )
Figure 6: Dependence of the average and variance of D p ( T ) and R p ( T ) as a function ofcommunity structure p (logarithmic scale). Average (circles) and one standard deviationadded to and subtracted from the average (dash curves) are calculated from 20 stochasticrealizations for each fixed value of p . The results shown in Panels (a) and (b) assess a “worstcase scenario” of having on average every recovered node going into relapse with probabilitynearly one, in symbols 1 − e − ρ τ ( t ) (cid:104) k (cid:105) ≈
1. 22 − − − − M ean D p ( T ) ! =3 ! =5 ! =7 ! =10 ! = " − − − − M ean R p ( T ) ! =3 ! =5 ! =7 ! =10 ! = " (a) (b) Figure 7: Average D p ( T ) and R p ( T ) as functions of the community structure, p . Panels (a)and (b) display the results obtained from using a time-dependent relapse rate ρ τ ( t ). Therelapse rate jumps from 0.90 to 0.12 at time t = τ , that is, every node diminishes its proba-bility of transition from the recovered into the problem drinker state by half (probabilities gofrom 1 − e − . (cid:104) k (cid:105) ≈ − e − . (cid:104) k (cid:105) ≈ .
5) Panels (a) and (b) show the changes in averagesas a function of the timing in the jump ( τ ). The relapse reduction at times, τ = 3 (upwardtriangles), τ = 5 (diamonds), τ = 7 (right triangles), τ = 10 (circles) are highlighted. Theaverages displayed in Figure 6 are for the case τ = ∞∞