- Email: [email protected]
Static behavioral effects on gonorrhea transmission dynamics in a MSM population
Journal of Theoretical Biology 267 (2010) 35–40
Contents lists available at ScienceDirect
Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi
Static behavioral effects on gonorrhea transmission dynamics in a MSM population Benjamin R. Morin a,,1, Liana Medina-Rios b, Erika T. Camacho c, Carlos Castillo-Chavez a,1,2 a
School of Human Evolution and Social Change, SHESC 233, P.O. Box 872402, Tempe, AZ 85287-2402, USA Department of Mathematics and Statistics, University of Maryland, Baltimore County, 1000 Hilltop Circle, Baltimore, MD 21250, USA c Division of Mathematical and Natural Sciences, Arizona State University, 4701 W. Thunderbird Rd, Glendale, AZ 85306, USA b
a r t i c l e in fo
abstract
Article history: Received 22 September 2009 Received in revised form 28 June 2010 Accepted 21 July 2010 Available online 27 July 2010
An SIS/SAS model of gonorrhea transmission in a population of highly active men-having-sex-with-men (MSM) is presented in this paper to study the impact of safe behavior on the dynamics of gonorrhea prevalence. Safe behaviors may fall into two categories—prevention and self-awareness. Prevention will be modeled via consistent condom use and self-awareness via STD testing frequency. Stability conditions for the disease free equilibrium and endemic equilibrium are determined along with a complete analysis of global dynamics. The control reproductive number is used as a means for measuring the effect of changes to model parameters on the prevalence of the disease. We also find that appropriate intervention would be in the form of a multifaceted approach at overall risk reduction rather than tackling one specific control individually. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Mathematical epidemiology Gonorrhea
1. Introduction The task of disease eradication and prevention, while of interest to a society as a whole, is in the hands of ‘‘active’’ individuals. Active here can be taken to mean any behavior or circumstance that puts that individual in the process of the disease spreading. Of interest here is the spread of sexually transmitted diseases (STDs), more specifically gonorrhea, and thus active is taken to mean highly sexually active. The study of so-called core groups has been used, to some success, in modeling the prevalence of an STD (Hethcote and Van Ark, 1992; Hethcote and Yorke, 1984; Dietz, 1988; Dietz and Hadeler, 1988; CastilloChavez, 1989). In classical model controls are viewed as amalgamations of several affects that are treated via a single parameter (i.e. b as a ‘‘force of infection’’) (Hethcote and Van Ark, 1992; Hadeler, 1989; Bailey, 1975; Anderson, 1982, 1991). The behavior of individuals in a population can limit the spread of disease and we propose to place them into two categories—prevention and self-awareness. Prevention is specific to the characteristics of each disease. With STDs such as gonorrhea, prevention education advocates safe sex practices, condom use for a larger proportion of the time and reduction in the number of sexual partners. Self-awareness is Corresponding author.
E-mail address: [email protected] (B.R. Morin). Mathematical and Computational Modeling Sciences Center, Arizona State University, P.O. Box 871904, Tempe, AZ 85287, USA. 2 Santa Fe Institute, SFI, 1399 Hyde Park Road, Santa Fe, NM 87501, USA. 1
0022-5193/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2010.07.027
simply possessing a knowledge of symptoms, the presence of asymptomatic dynamics, and getting tested and treated. Self-awareness leads to frequent STD screening as a disease control method since individuals may be infected but have no knowledge of their own status. In particular, infected individuals may have asymptomatic gonorrhea but still be equally infectious to the symptomatic (Hethcote and Yorke, 1984; CDC, 2007). If an infected individual is asymptomatic, but relatively uneducated about the disease, they will not seek medical attention and thus their average infectious period is potentially much greater than that of a symptomatically infected individual. Typically, symptomatically infected individuals recover 14 days after the start of treatment which, due to the pain associated with the symptoms, usually begins once symptoms show a few days after gonorrhea is contracted (CDC, 2007). Seminal mathematical work was done by Hethcote and Yorke (1984) on heterosexual gonorrhea transmission. There, a two-sex model of symptomatically infected, asymptomatically infected and susceptible populations with different activity levels was used. The focus was the effect that contact tracing and increased STD testing would have on the dynamics of the disease. It was found that contact tracing of infectees was less effective than that of infectors. Here the infectors were identified as a core group, or highly sexually active subpopulation. Hethcote concluded that focusing on the core group’s activities was key to controlling the spread of gonorrhea. Li et al. specifically modeled STDs like gonorrhea using an SIS model with a multiple strains and varied reaction to infection (Li et al., 2003). It was found that sufficient heterogeneity of the
36
B.R. Morin et al. / Journal of Theoretical Biology 267 (2010) 35–40
female in the form of contact structure, immune response or activity level was necessary in order to have coexistence of the multiple strains. The conditions for the existence and stability of an endemic equilibrium with two strains were found. Coexistence demanded that one strain should be better at infecting one subpopulation while the other be better at infecting another. It was concluded that each strain creates reservoirs in the population that it is less able to infect. Although not specifically looking at gonorrhea, previous mathematical epidemiological studies by Kemper (1978) have developed a general model for curable diseases with symptomatic or asymptomatic infection. There an SIS/SAS model was developed to consider the impact of asymptomatic attacks. However, this model does not treat the different recovery time of asymptomatic infected individuals, nor does it account for the different proportion of contacts with symptomatically infected individuals that lead to symptomatic versus asymptomatic infection and vice versa. As a case study in the effects of behavior change in the spread of a curable disease that confers no permanent immunity, we present an SIS/SAS model of gonorrhea transmission in a men-having-sex-with-men (MSM) population that incorporates the effects of education previously mentioned. The choice of the MSM population was because it is a population that is typically high risk and very active (Hethcote and Yorke, 1984) and can make occasional contacts with the highly susceptible female population. Thus focusing on the eradication of STD in this high risk and active population can shed light on STD control measures in the population as a whole. In this work we assume that the occasional male–female contacts in the MSM population is negligible from the point of view of spread within the MSM population but has the potential to cause gonorrhea to spread in the larger population not modeled here. The influence of safe sex will be modeled exogenously with a weighted average on the effective contact rate that accounts for changing behavior with respect to condom use, both of which are fixed parameter values. Self-awareness will be modeled via an increase in infectious period for asymptomatic versus symptomatic individuals with the asymptomatic individuals exhibiting either high or low levels of awareness (frequent or rare testing). Analysis is done describing the disease free and endemic dynamics as well as quantifying the changes necessary to eradicate gonorrhea in this population.
2. Basic homosexual gonorrhea transmission model 2.1. Mathematical model The population modeled is single sex with three homogeneous compartmental states available: susceptible individuals, S, symptomatically infected, I, and asymptomatically infected individuals, A. The model equations are dS IþA ¼ mðNSÞðl1 þ el2 ÞbS þ aI þ ðrl3 þ sl4 ÞaA, dt N dI I A ¼ ðl1 þ el2 ÞbS q1 þ ð1q2 Þ ðm þ aÞI, dt N N dA I A ¼ ðl1 þ el2 ÞbS ð1q1 Þ þ q2 m þ ðrl3 þ sl4 aÞA: dt N N
l2 . The proportion of time safe sex prevents the transmission of gonorrhea is ð1eÞ. Thus for the proportion of the time an individual participates in safe sex, l2 is multiplied by the reduction factor, e. These parameters combine to become a total reduction factor on the force of infection, i.e. ðl1 þ el2 Þb. Individuals are recruited into the I class from the S class at a rate ðl1 þ el2 ÞbSðq1 I=N þ ð1q2 ÞA=NÞ, where q1 and (1 q2) are the proportion of individuals that become symptomatically infected from contact with symptomatically infected individuals and asymptomatically infected individuals, respectively. Similarly individuals are recruited into the A class from the S class at a rate ðl1 þ el2 ÞbSðð1q1 ÞI=N þ q2 A=NÞ, similarly where q2 and (1 q1) are the proportion of individuals that become asymptomatically infected from contact with asymptomatically infected individuals and symptomatically infected individuals, respectively. Individuals from the I class reenter the susceptible class due to treatment at a rate a. Since individuals in the A class do not know they have gonorrhea the average duration of infection of an asymptomatic person is assumed to be longer than that of a symptomatic individual. We represent this increase in infectious period via rl3 þ sl4 A ½0,1. The proportion of people who engage in low levels of self-awareness, do not get tested often, is denoted by l3 while the population of those who have a high level of selfawareness are represented by l4 . The parameters r and s are extensions to the symptomatic treatment rate a caused by low and high levels of self-awareness, respectively, and by the nature of their meanings r o s. Thus reentry into the susceptible population from the asymptomatic class occurs at the rate ðrl3 þ sl4 ÞaA. The model assumes a constant population size, N, with constant recruitment and removal, m. This limits our ability to extend the time scale of this work but allows for analysis to be carried out that would not be possible. For this reason we are choosing college bound (a time of high-sexual activity and educational influence) as a further targeting of the population. The reasonable time scale for the model is thus over 4 years, where it is assumed that model parameters cannot change and the population will hold relatively constant. We also do not have memory in our population. That is, once diagnosed with asymptomatic gonorrhea an individual is not more likely to get frequent testing (Hethcote and Yorke, 1984). Here, new individuals will enter the system as old ones leave making the absence of memory an acceptable assumption. We present control here as a meaningful separation of the more classical parameters, but these are still averaged measures of behavior. Thus, we cannot make any tracking of how often a specific individual has been diagnosed asymptomatic or has become infected at all. We also do not look at dynamic changes in sexual mixing through the use of evolving contact structures. A few proposed generalizations will be made in the end of the paper as veins of future work but the thought governing this work is to start with as little detail as seems relevant.
2.2. Parameter estimation
ð1Þ
The rate at which the susceptible population is lost to infection is bSðI þ AÞ=N, where b is the effective contact rate. We introduce control into the system via l1 , l2 , and e. The proportion of the population participating in non-safe sex, or equivalently the proportion of time an individual participates in non-safe sex, is l1 and the proportion of the population participating in safe sex is
The system we wished to model was a college attending (18–22 years old), MSM population thus m is taken to equal 14 years. This group was chosen due to the prevalence of both highly sexually active individuals and the ability for policy changes to potentially reach the entire population. The behavioral parameters desired for the model were not found from a single source, or over a single demographic across many sources, and thus we have put together several sources’ worth of parameter estimates as an approximation to the potential behavior of the study group. The effective contact rate is the product of the number of risky
B.R. Morin et al. / Journal of Theoretical Biology 267 (2010) 35–40
contacts per year and the proportion of these that lead to infection. Based on the information given by NIH (2000), roughly 50% of MSM contacts with an infected individual leads to a new infection in the susceptible partner. Making an estimate of 100 risky contacts per year we arrive at b ¼ 50 infectious contacts per year. To determine the amount of time spent practicing safe versus risky behavior, according to a study of an MSM population by Shlay et al. (2004), 25.6% of participants report consistent condom usage, thus l2 ¼ 0:256 and l1 ¼ 1l2 ¼ 0:744. To determine e we take into account that although condoms are 97% effective at preventing gonorrhea infection with perfect use, many uses are imperfect due to slippage, breaking, etc. According to Shlay et al. and Stone et al. (1989), there is a 16.6–17.3% usage failure of condoms in a MSM population. We chose to let e 0:173. According to NIH (2000) among men, symptoms 4 6 typically appear within 4–6 days, 365 2365 years, following 7 infection while treatment duration is 7 days, 365 years, thus we take a 365 12 years. A study by Fortenberry et al. showed that roughly 50% of men got tested for STDs, including gonorrhea, annually (Fortenberry et al., 2002). This study was done without differentiating between hetero- or homo-sexual males and the age range was larger than that of your typical college bound of 18–22, with about 50% of their sample in our age range. This therefore serves as a poor estimate, but an estimate none the less, and is probably an overestimate of safe behavior. Using this information we let l3 ¼ l4 ¼ 0:5. For the values of r and s we will use testing once every year and once every 3 months, the CDC recommendation for highly active MSM (CDC, 2007). Thus we 365 12 1 have r a ¼ 365 372 years and s a ¼ 99 years, or r ¼ 372 ¼ 31 and 12 4 s ¼ 99 ¼ 33. Making r represent a year as asymptomatic is probably an underestimate for this group, but it is our hope that this models some form of ‘‘college’’ pressure: student health efforts, peer pressure, partner desire. The probabilities q1 and q2 are two difficult parameters to estimate. If data were collected that found the probability of exhibiting a certain type of infection given a successfully infectious contact with an infectious individual (not necessarily of the same type) then we would immediately have these parameters. However, what we generally have is some presence proportion (i.e. a percentage of the infected population who is of a particular type). Thus, the probabilities would solve the necessary relationship in assumed endemic populations. According to Bozicevic et al. (2006) 29.9% of infected individuals are asymptomatic. Thus we would have to solve 0:299I1 ¼ 0:701A1 for q1 and q2. This is daunting since the expression for the endemic equilibrium is very large and the result would not necessarily be an invertible function (i.e. more than one pair of the probabilities would give the same relation). A numerical investigation, using the nominal values for every other parameter and varying q1 and q2, provides a few values for the respective probabilities that are near the 29.9% estimate. Choosingq1 ¼0.965 and q2 ¼ 0.034 we get an endemic percentage of asymptomatic individuals of 29.85%.
3. Analysis 3.1. Global stability analysis With the model constructed we wish to do a full analysis on the qualitative behavior of the system. To that end the main focus of this section is to prove Theorem 1. The discussion proving this will involve a treatment of the stability of the disease free equilibrium, conditions on the number of endemic equilibrium that may exist, and a preclusion of closed orbits which will make all stability arguments global.
37
Theorem 1 (Global stability). System (1) has two fixed points: a disease free and an endemic equilibrium. The disease free equilibrium is globally stable when the control reproductive number is less than one and unstable when the control reproductive number is greater than one. The endemic equilibrium does not exist biologically when the control reproductive number is less than one and is globally stable when the control reproductive number is greater than one. To start we wish to formulate System (1) into a slightly easier form. Since S(t)+ I(t)+ A(t)¼N we may eliminate one state variable for the purpose of analysis. We may also rescale in both state and time to reduce the overall system. Using xðtÞ ¼ I=N, yðtÞ ¼ A=N, and t ¼ t obð1xyÞ, where o ¼ l1 þ el2 and p ¼ rl3 þ sl4 , we arrive at dx q1 x , ¼ q1 x þ ð1q2 Þy dt RII ð1xyÞ dy q2 y , ¼ ð1q1 Þx þ q2 y dt RAA ð1xyÞ
ð2Þ
where RII ¼ boq1 =ðm þ aÞ and RAA ¼ boq2 =ðm þ paÞ. Since each state variable S, I and A are positive we have that x þ y r1. Thus the rescaling is positive invariant. There is a disease free equilibrium that always exists, DFE :¼(x*, y*)¼ (0,0), implying S* ¼N, I* ¼0 and A* ¼0. To determine the local stability of the DFE the system is linearized about (x*, y*) resulting in the following: 0 B Jðx ,y Þ ¼ B @
q1 q1 RII ð1q1 Þ
ð1q2 Þ
1
C C q2 A: q2 RAA
ð3Þ
The characteristic polynomial of the above Jacobian is q q q q l2 q1 þq2 1 2 l þ q1 1 q2 2 ð1q1 Þð1q2 Þ, RII RAA RII RAA 2
which is in the form l bl þ c. It may be shown that conditions for the determinant of the Jacobian to be positive, c 4 0, are identical to RE, the basic control number, being less than one. First, we use the next generation operator to compute the number of secondary infections a typical infectious individual creates in a completely susceptible population (van den Driessche and James, 2002). This method essentially creates two vectors F and V. The vector F contains any rates in the infected classes’ ODEs that constitute recruitment from a non-infected class. The vector V contains the negative of all other rates except for the rates in noninfected classes that stand for recruitment into the infected classes. This is more easily done with Eqs. (1) than the rescaled system in Eqs. (2). This gives us 0
1 0 B I A C B obS q1 þ ð1q2 Þ C B N N C F ¼B , C B I A C @ A obS ð1q1 Þ þq2 N N 0 B V ¼@
mðNSÞaIpaA
1
ðm þ aÞI
C A:
ðm þpaÞA
ð4Þ
ð5Þ
We then define F and V as the Jacobians of F and V, respectively. Specifically these are matrices where the (i, j)th element is the partial derivative of the ith term of the vector with respect to the jth state variable. The spectral radius, or dominant eigenvalue, of FV 1 is then our basic reproductive number. The numerator of
38
B.R. Morin et al. / Journal of Theoretical Biology 267 (2010) 35–40
each term in Eq. (6) is from F and each denominator is from V 1 : 2 3 boq1 boð1q2 Þ " # 6 aþm RII RAI m þpa 7 6 7 1 FV ¼ 6 ¼ : ð6Þ boq2 7 RIA RAA 4 boð1q1 Þ 5 aþm m þpa Interestingly, the terms in Eq. (6) are each control reproductive numbers in and by themselves. The term Rij is the average number of individuals recruited into class j from a typical member in class i per unit time. Thus the control reproductive number of the entire system is a function of each of these individual recruitment thresholds. The spectral radius of Eq. (6) is qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi RII þ RAA þ ðRII RAA Þ2 þ 4RIA RAI RE ¼ : ð7Þ 2
Corollary 1. The condition RE o1 is identical to c 40 for the 2 characteristic polynomial l bl þ c for System (2). Proof. We begin with the condition for RE o1: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi RII þ RAA þ ðRII RAA Þ2 þ 4RIA RAI o 1, 2
Corollary 2. If q1 þ q2 o1 then RE o 13C o 0 and RE 413C 40. If q1 þq2 41 then RE o 13C 4 0 and RE 413C o0. Proof. If q1 þ q2 o 1 then C and c have differing sign and thus if c 40 then C o 0 and their relationships to RE are the opposite of one another. If q1 þ q2 41 then C and c have the same sign and their relationships to RE are identical. & The equation f(z) is a quadratic and thus may have 0, 1, or 2 roots in (0,1). The relative signs of f(0) and f(1) will allow us to determine under what conditions f(z) has a particular number of roots in the unit interval. Consider f(0) ¼C and f(1)¼1+ C B. We already have that the sign of C may be viewed as being dependent on the magnitude of RE. We also have that 1þ CB ¼
q1 q2 1 , 1q1 q2 RII RAA
whose sign depends on 1 q1 q2. Thus in order to study the zeros of f(z)¼ z2 Bz+ C we must consider all four combinations of the sum of q1 and q2 with the magnitude of the control reproductive number.
ðRII RAA Þ2 þ4RIA RAI o ð2ðRII þ RAA ÞÞ2 ,
1. RE o 1 and q1 þ q2 o1: In this situation f ð0Þ ¼ C o0 and f ð1Þ ¼ 1 þCB o0. Thus there are no zeros of f ðzÞ A ð0,1Þ. 2. RE 4 1 and q1 þ q2 o1: Here we have that f ð0Þ 4 0 and f ð1Þ o0. Thus there is a single root for f ðzÞ A ð0,1Þ. 3. RE o 1 and q1 þ q2 4 1: In the most difficult situation we have that f ð0Þ 4 0 and f ð1Þ 4 0. The quadratic having exactly two roots occurs when (a) B=2A ð0,1Þ, and (b) B2 4C 4 0.
R2II 2RII RAA þR2AA þ4RIA RAI o 44ðRII þ RAA Þ þðRII þ RAA Þ2 , 4RII RAA þ 4RIA RAI þ4ðRII þRAA Þ o 4, RII þ RAA þ RII RAA
where c is from the characteristic polynomial. If z¼0 then x¼ y¼0, the DFE. The interest thus lies in where f(z)¼0. A relationship between C and RE can be made using the existing relationship found in Corollary 1.
1q1 q2 o1: q1 q2
If we then consider the condition for c 40: q1 q2 q1 q2 ð1q1 Þð1q2 Þ 40, RII RAA
In order for B 4 0 we require that q1 q2 ð1RII RAA ÞRII RAA ð1q1 q2 Þ 4 0, 1RII RAA 4 RII RAA
RII þ RAA 2ðq1 þ q2 1Þ o : q1 q2 RII RAA
1q1 q2 , q1 q2
The second condition yields
1q1 q2 RII þ RAA þ RII RAA o1, q1 q2 then we see Eq. (8) identical to the condition for RE o 1.
ð8Þ &
Thus we have a somewhat easier condition for stability and may define 1q1 q2 R^ E ¼ RII þ RAA þRII RAA : q1 q2 This quantity is somewhat easier to understand and its relationship to 1 is identical to RE’s. If RE 4 1 then the DFE is unstable, but we have not discussed how many equilibrium may exist. Consider z¼x+y and the fact that if dz=dt ¼ 0 and dy=dt ¼ 0 then we would be at a fixed point for System (2). Solving dz=dt ¼ 0 we get an expression for y in terms of z by noting x¼z y. Plugging this expression into dy=dt and solving the new expression for zero we get zf(z)¼0, where f(z)¼z2 Bz+C, 1 1 q1 q2 þ RAA RII , B ¼ 2þ 1q1 q2 and 1 1 ð1q1 Þð1q2 Þq1 q2 1 1 c RII RAA ¼ , C¼ 1q1 q1 1q1 q2
ðRII þRAA Þ2 4ðq1 þq2 1Þ 4 : q1 q2 RII RAA
ð9Þ
These two conditions are contradictory. To illustrate the contradiction we invoke Eq. (7) to get that RE o 1-RII þ RAA o2. Rearranging Eq. (9) results in RII þ RAA 4 ðq1 þq2 1Þ 4 , RII þ RAA q1 q2 RII RAA
ð10Þ
but 4=ðRII þ RAA Þ 4 2 which results in the contradiction. Thus there are no zeros of f ðzÞ A ð0,1Þ. 4. RE 4 1 and q1 þ q2 4 1: In this situation f ð0Þ 4 0 and f ð1Þ o0. Thus there is a single root for f ðzÞ A ð0,1Þ. Combining the arguments gives us that if RE o1 then the only solution that is biologically relevant for our system is the DFE. When RE 4 1 the two solutions, which are biologically relevant, the unstable DFE and a single endemic equilibrium, EE. The entire above argument is valid only if 1q1 q2 a 0. When 1q1 q2 ¼ 0, the system exhibits a single EE, ! RII ðRII þRAA 1Þ RAA ðRII þ RAA 1Þ ðx,yÞ ¼ , , ðRII þ RAA Þ2 ðRII þ RAA Þ2 which is only valid if RII þ RAA 1 Z 0. We wish to make assertions about the stability of the EE without having to do a linearization
B.R. Morin et al. / Journal of Theoretical Biology 267 (2010) 35–40
around the fixed point which is very term intensive. We may now disprove the existence of closed orbits in the plane and thus assert that when the EE exists, RE 4 1, it is globally stable. Corollary 3. System (1) has no closed orbits. Proof. By Dulac’s criterion, if there exists a function j A C 1 such that _ _ @ðjxÞ @ðjyÞ þ a 0, @x @y _ y_ has no closed orbits. Let j ¼ 1=xy. then the planar system x, Then _ @ðjxÞ @ q1 ð1q2 Þ q1 ¼ þ @x y x @x RII ð1xyÞy ð1q2 Þ q1 ¼ , x2 RII ð1xyÞ2 y _ @ðjyÞ @ ð1q1 Þ q2 q2 þ ¼ @y y @y x RAA ð1xyÞx ð1q1 Þ q2 ¼ : y2 RAA ð1xyÞ2 x _ _ þ @ðjyÞ=@y is Since q1 ,q2 A ½0,1 and RAA ,RII 4 0 the sum @ðjxÞ=@x always negative. Thus there are no closed orbits for System (2). Since the dynamics are identical for System (1) we have precluded limit cycles in it and have shown what was intended. & 3.2. Sensitivity of gonorrhea transmission model In a perfect world, public health initiatives would be simple, multifaceted, and have great effects on the dynamics of a disease. However, this is not always the case and moreover, economic costs have to be considered in determining which interventions to support. Another important question is whether the size of the population is important to the dynamics of the disease and intervention. In order to address these concerns, in this section we examine the sensitivity of the system to changes in the parameter values. We take two approaches, first the time-dependent sensitivity of the original system to changes in parameter values is calculated via elasticity and second we compute what changes to the parameters would be necessary to bring the control reproductive number to a value less than 1.
39
difference quotient to arrive at i @X Xðt; y^ ÞXðt; y^ Þ ðtÞ : @yi Di
ð11Þ
We may use this approximation in our computation of the elasticity of J(t) with respect to each parameter. The sensitivity of the total infected population with respect to the nominal parameters is presented in Fig. 1 with Di ¼ 0:001yi used in each of the approximations. Not surprisingly b acts as a large control on the infected population. Small changes in b will result in the largest changes to J(t) and increase the number of individuals infected. The parameter modeling non-condom use, l1 , intuitively is also very large. The remaining four control parameters, l3 , r, s, and e are all near zero and thus each will require large changes to impact the infected population’s size.
3.2.2. Reducing RE Holding all but one parameter constant we varied the free parameter until RE o 1. Each parameter was varied by 0.1% in the correct direction to reduce the infected population, analogous to reducing RE. Of these only two single effort campaigns resulted in a control reproductive number less than 1, l1 and b. Reducing l1 to roughly 0.392, condom use increases to 61% of the time, or b to 28.239, 56.5 risky contacts per year, will result in RE o 1. Each represents a fairly significant change in behavior. Even with the entire population following the CDC recommended testing frequency of once every 3 months the disease will remain prevalent. Instead if one considers a multifaceted campaign to disease reduction we may vary each parameter at once, by 1/10%, until RE o1. This will represent smaller changes to each parameter. However, we fixed s at the CDC recommendation to measure its effectiveness. If the control parameter vector reads ðl1 , l3 , b, r, eÞ ð0:5605,0:3767,37:67,0:0428,0:1303Þ then gonorrhea will no longer exist within the population. This translates to 43.95% of the population using a condom, up from 25.6%, 62.3% of the population using the CDC’s testing frequency recommendation, reduce the number of risky a contacts per year from 100 to 75.34, even if you do not follow the CDC recommendation for testing get tested once every 9 months and finally reduce condom failure rate from 17.3% to 13%.
3.2.1. Sensitivity of I(t)+ A(t) In order to discuss the importance of individual parameters one must investigate how their value affects the state variables over time. We will consider the concept of elasticity. Formally, one may define elasticity of a function Xðt; yÞ with respect to a parameter yi A y as Eyi ðtÞ ¼
yi Xðt; yÞ
@Xðt; yÞ : @yi
This measures the ratio of a percent change in a parameter to that of the function. This gives us a unit-less and scaled method with which to compare each parameter’s affect on the solution J(t):¼I(t)+ A(t). However, we do not have a closed form for J(t) and thus we must make an approximation to both its value and its partial derivative. In general, consider dX=dt ¼ f ðt,X; yÞ where y is a parameter vector. Now consider the vector of nominal parameter values, y^ , and a small perturbation, Di , of the ith element, y^ i , and call this i i new parameter vector y^ . If what we are interested in is @X=@y near the nominal value then we could do the following. i Numerically find the solutions Xðt; y^ Þ and Xðt; y^ Þ then take the
Fig. 1. Sensitivity of the infected population, regardless of symptom expression, with respect to each model parameter using the nominal values ðm, l1 , l3 , b, a, r, s,q1 ,q2 , eÞ ¼ ð1=4,0:744,0:5,50,365=12,1=31,4=33,0:965,0:034,0:173Þ and N¼ 10 000 with S(0) ¼ 99 999, I(0) ¼0, and A(0)¼ 1.
40
B.R. Morin et al. / Journal of Theoretical Biology 267 (2010) 35–40
4. Conclusion and discussion There are a few shortcomings in this model that need to be addressed. The removal of a constant population size could potentially remove the monotonicity of the system and thus would remove the hope for any kind of qualitative global stability analysis via criteria as supplied herein. It would fall upon the modeler to do a more complex stability analysis where cycles may develop. It is unclear to the authors here if gonorrhea prevalence in this type of population merits the need for cyclic dynamics. In this paper we have considered the average, safe and risky behavior of the population. While this has given insights into how education affects disease transmission, it does not explicitly treat the heterogeneity present in such populations. This is a step forward from more classical models where the individual affects were not treated explicitly, but more work may be done for, if not analytical then numerical, results to suggest policy. Several treatments could simply increase the number of compartments while still not dealing with individual level concerns. Specifically, we may take each of the S, I, and A classes and divide them into those who are ‘‘safe’’ and those that are ‘‘risky’’ with respect to each behavior resulting in a system of 12 equations (i.e. susceptible population who are risky with respect to protection and self-awareness, risky with only one behavior or risky with neither) with flow between like-state groups of different ‘‘risk’’ types. Work done in Castillo-Chavez (1989), Blythe et al. (1991), Blythe and Castillo-Chavez (1989) and Morin et al. (accepted for publication) has also shown the importance of affinity-based mixing while others have turned to simulation to look at prevention strategies on networks (see Kemper, 1978). Here we have proportionate mixing which should not be accurate given the pain of symptoms and the mixing assumed between risky and not risky individuals. Furthermore, there is this intuitive desire to model stubbornness in a system such as this. Two difficulties with educational measures is how broad of outreach one should invest in, how many people get to see the information, and how intense the education should be. The trade offs between should we remind everyone a little or remind a few constantly can be addressed with more game theoretic or individual based approaches (Shim and Galvani, 2009; Shim et al., 2010). In this paper, a simple model of homosexual gonorrhea transmission between men is analyzed and the parameter sensitivities are determined. A caveat to the sensitivity portion of analysis is that the results are based off of a set of imperfect parameters. The populations studied in Shlay et al. (2004), Stone et al. (1989) and Fortenberry et al. (2002) are non-identical in age, region or sexuality and thus the parameters chosen are done as a starting point. The Fortenberry study covered men from ages 14 to 59, with mean age 24 with no distinction made for sexuality, from seven different US cities. The Shlay study covered MSM subjects with mean age 33.4 from a single clinic in Denver. The Stone paper includes a sample of homosexual men where half of which are between the ages of 18 and 30 participating in a multicity HIV study. Work could be put to good use if it sampled all the necessary data from several groups of men to determine if the behaviors are similar across region for the ages 18–22. In Section 3.2, it is shown that disease education is useful in disease control, however, education should be multifaceted in order to reduce the need of changing any one parameter by too much. In detail, the system is most sensitive to b (the number of risky contacts per year), and this parameter, as well as l1 , may be changed all on its own enough to eradicate gonorrhea. These changes, however, appear to represent too drastic a change to be realistic. The sensitivity analysis has also revealed that the two methods of prevention, frequent condom use and the number of partners taken, have not only the strongest effect on the infectious population but the methods of self awareness, the proportion who participate in frequent testing and the frequency of this testing,
have a sensitivity more than 15 times smaller. This supports a ‘‘prevention first’’ intuition to reducing the prevalence of the disease. We calculated time-dependent sensitivities, and although gonorrhea is best viewed as an endemic disease there are fouryear residence systems, like a high school, where the introduction of a single infected individual is feasible. In such situations we have shown that within the first 6 months prevention measures serve the goal of disease eradication best since after a few years the effect is approximately one seventh what it was at its highest.
Acknowledgements This project has been partially supported by grants from the National Science Foundation (NSF - Grant DMPS-0838704), the National Security Agency (NSA - Grant H98230-09-1-0104), the Alfred P. Sloan Foundation and the Office of the Provost of Arizona State University. The authors would like to thank the suggestions and effort put forth by Carlos Castillo-Chavez, Christopher Kribs-Zaleta and Baojun Song in order to accomplish the analysis that has been done to the model. References Anderson, R.M., 1982. Population Dynamics of Infectious Diseases: Theory and Applications. Chapman and Hall, London, New York. Anderson, R.M., 1991. Infectious Diseases of Humans. Oxford Science Publications, Great Britain. Bailey, N.T.J., 1975. The Mathematical Theory of Infectious Diseases and its Applications. Griffin, London. Blythe, S.P., Castillo-Chavez, C., 1989. Like-with-like preference and sexual mixing models. Math. Biosci. 96, 221–238. Blythe, S.P., Castillo-Chavez, C., Palmer, J., Cheng, M., 1991. Towards a unified theory of mixing and pair formation. Math. Biosci. 107, 379–405. Bozicevic, I., et al., 2006. Epidemiological correlates of asymptomatic gonorrhea. Sexually Transmitted Dis. 33 (5), 289–295. Castillo-Chavez, C. (Ed.), 1989. Mathematical and Statistical Approaches to AIDS Epidemiology, Lecture Notes in Biomathematics, vol. 83, Springer-Verlag, New York. Centers for Disease Control, CDC Fact Sheet: Gonorrhea. December 2007. Dietz, K., 1988. On the transmission dynamics of HIV. Math. Biosci. 90, 397–414. Dietz, K., Hadeler, K.P., 1988. Epidemiological models for sexually transmitted diseases. J. Math. Biol. 26, 1–25. Fortenberry, et al., 2002. Relationships of stigma and shame to gonorrhea and HIV screening. Am. J. Public Health 92 (3). Hadeler, K.P., 1989. Modeling AIDS in structured populations. In: 47th Session of the International Statistical Institute, Paris, August/September, Conference Proceedings C1-2, pp. 83–99. Hethcote, H.W., Yorke, J.A., 1984. Gonorrhea Transmission Dynamics and Control, Lecture Notes in Biomathematics, vol. 56. Springer-Verlag, New York. Hethcote, H.W., Van Ark, J.W., 1992. Modeling HIV Transmission and AIDS in the United States, Lecture Notes in Biomathematics, vol. 95. Springer-Verlag, New York. Kemper, J., 1978. The effect of asymptomatic attacks on the spread of infectious disease: a deterministic model? Bull. Math. Biol. 4 (I) 707–718. Li, J., Ma, Castillo-Chavez, C., 2003. Coexistence of pathogens in sexuallytransmitted disease models. J. Math. Biol. 47, 547–568. doi:10.1007/s00285003-0235-5. Morin, B., et al., accepted for publication. Notes from the heterogeneous: a few observations on the implications and necessity of affinity. J. Biol. Dyn. National Institute of Allergy and Infectious Diseases, National Institutes of Health, Department of Health and Human Services. Workshop Summary: Scientific Evidence on Condom Effectiveness for Sexually Transmitted Disease (STD) Prevention, June 2000. Shim, E., Galvani, A., 2009. Impact of transmission dynamics on the costeffectiveness of rotavirus vaccination. Vaccine 27, 4025–4030. Shim, E., Chapman, G., Galvani, A., 2010. Decision making with regard to antiviral intervention during an influenza pandemic. Medical Decision Making 30 (4), E64–E81. Shlay, J., et al., 2004. Comparison of sexually transmitted disease prevalence by reported condom use: errors among consistent condom users seen at an urban sexually transmitted disease clinic. Sexually Transmitted Dis. 31 (9), 526–532. Stone, E., et al., 1989. Correlates of condom failure in a sexually active cohort of men who have sex with men. J. Acquired Immune Defic. Syndr. Human Retrovirology 20 (5), 495–501. van den Driessche, P., Watmough, J., 2002. Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 180, 29–48.
Contents lists available at ScienceDirect
Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi
Static behavioral effects on gonorrhea transmission dynamics in a MSM population Benjamin R. Morin a,,1, Liana Medina-Rios b, Erika T. Camacho c, Carlos Castillo-Chavez a,1,2 a
School of Human Evolution and Social Change, SHESC 233, P.O. Box 872402, Tempe, AZ 85287-2402, USA Department of Mathematics and Statistics, University of Maryland, Baltimore County, 1000 Hilltop Circle, Baltimore, MD 21250, USA c Division of Mathematical and Natural Sciences, Arizona State University, 4701 W. Thunderbird Rd, Glendale, AZ 85306, USA b
a r t i c l e in fo
abstract
Article history: Received 22 September 2009 Received in revised form 28 June 2010 Accepted 21 July 2010 Available online 27 July 2010
An SIS/SAS model of gonorrhea transmission in a population of highly active men-having-sex-with-men (MSM) is presented in this paper to study the impact of safe behavior on the dynamics of gonorrhea prevalence. Safe behaviors may fall into two categories—prevention and self-awareness. Prevention will be modeled via consistent condom use and self-awareness via STD testing frequency. Stability conditions for the disease free equilibrium and endemic equilibrium are determined along with a complete analysis of global dynamics. The control reproductive number is used as a means for measuring the effect of changes to model parameters on the prevalence of the disease. We also find that appropriate intervention would be in the form of a multifaceted approach at overall risk reduction rather than tackling one specific control individually. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Mathematical epidemiology Gonorrhea
1. Introduction The task of disease eradication and prevention, while of interest to a society as a whole, is in the hands of ‘‘active’’ individuals. Active here can be taken to mean any behavior or circumstance that puts that individual in the process of the disease spreading. Of interest here is the spread of sexually transmitted diseases (STDs), more specifically gonorrhea, and thus active is taken to mean highly sexually active. The study of so-called core groups has been used, to some success, in modeling the prevalence of an STD (Hethcote and Van Ark, 1992; Hethcote and Yorke, 1984; Dietz, 1988; Dietz and Hadeler, 1988; CastilloChavez, 1989). In classical model controls are viewed as amalgamations of several affects that are treated via a single parameter (i.e. b as a ‘‘force of infection’’) (Hethcote and Van Ark, 1992; Hadeler, 1989; Bailey, 1975; Anderson, 1982, 1991). The behavior of individuals in a population can limit the spread of disease and we propose to place them into two categories—prevention and self-awareness. Prevention is specific to the characteristics of each disease. With STDs such as gonorrhea, prevention education advocates safe sex practices, condom use for a larger proportion of the time and reduction in the number of sexual partners. Self-awareness is Corresponding author.
E-mail address: [email protected] (B.R. Morin). Mathematical and Computational Modeling Sciences Center, Arizona State University, P.O. Box 871904, Tempe, AZ 85287, USA. 2 Santa Fe Institute, SFI, 1399 Hyde Park Road, Santa Fe, NM 87501, USA. 1
0022-5193/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2010.07.027
simply possessing a knowledge of symptoms, the presence of asymptomatic dynamics, and getting tested and treated. Self-awareness leads to frequent STD screening as a disease control method since individuals may be infected but have no knowledge of their own status. In particular, infected individuals may have asymptomatic gonorrhea but still be equally infectious to the symptomatic (Hethcote and Yorke, 1984; CDC, 2007). If an infected individual is asymptomatic, but relatively uneducated about the disease, they will not seek medical attention and thus their average infectious period is potentially much greater than that of a symptomatically infected individual. Typically, symptomatically infected individuals recover 14 days after the start of treatment which, due to the pain associated with the symptoms, usually begins once symptoms show a few days after gonorrhea is contracted (CDC, 2007). Seminal mathematical work was done by Hethcote and Yorke (1984) on heterosexual gonorrhea transmission. There, a two-sex model of symptomatically infected, asymptomatically infected and susceptible populations with different activity levels was used. The focus was the effect that contact tracing and increased STD testing would have on the dynamics of the disease. It was found that contact tracing of infectees was less effective than that of infectors. Here the infectors were identified as a core group, or highly sexually active subpopulation. Hethcote concluded that focusing on the core group’s activities was key to controlling the spread of gonorrhea. Li et al. specifically modeled STDs like gonorrhea using an SIS model with a multiple strains and varied reaction to infection (Li et al., 2003). It was found that sufficient heterogeneity of the
36
B.R. Morin et al. / Journal of Theoretical Biology 267 (2010) 35–40
female in the form of contact structure, immune response or activity level was necessary in order to have coexistence of the multiple strains. The conditions for the existence and stability of an endemic equilibrium with two strains were found. Coexistence demanded that one strain should be better at infecting one subpopulation while the other be better at infecting another. It was concluded that each strain creates reservoirs in the population that it is less able to infect. Although not specifically looking at gonorrhea, previous mathematical epidemiological studies by Kemper (1978) have developed a general model for curable diseases with symptomatic or asymptomatic infection. There an SIS/SAS model was developed to consider the impact of asymptomatic attacks. However, this model does not treat the different recovery time of asymptomatic infected individuals, nor does it account for the different proportion of contacts with symptomatically infected individuals that lead to symptomatic versus asymptomatic infection and vice versa. As a case study in the effects of behavior change in the spread of a curable disease that confers no permanent immunity, we present an SIS/SAS model of gonorrhea transmission in a men-having-sex-with-men (MSM) population that incorporates the effects of education previously mentioned. The choice of the MSM population was because it is a population that is typically high risk and very active (Hethcote and Yorke, 1984) and can make occasional contacts with the highly susceptible female population. Thus focusing on the eradication of STD in this high risk and active population can shed light on STD control measures in the population as a whole. In this work we assume that the occasional male–female contacts in the MSM population is negligible from the point of view of spread within the MSM population but has the potential to cause gonorrhea to spread in the larger population not modeled here. The influence of safe sex will be modeled exogenously with a weighted average on the effective contact rate that accounts for changing behavior with respect to condom use, both of which are fixed parameter values. Self-awareness will be modeled via an increase in infectious period for asymptomatic versus symptomatic individuals with the asymptomatic individuals exhibiting either high or low levels of awareness (frequent or rare testing). Analysis is done describing the disease free and endemic dynamics as well as quantifying the changes necessary to eradicate gonorrhea in this population.
2. Basic homosexual gonorrhea transmission model 2.1. Mathematical model The population modeled is single sex with three homogeneous compartmental states available: susceptible individuals, S, symptomatically infected, I, and asymptomatically infected individuals, A. The model equations are dS IþA ¼ mðNSÞðl1 þ el2 ÞbS þ aI þ ðrl3 þ sl4 ÞaA, dt N dI I A ¼ ðl1 þ el2 ÞbS q1 þ ð1q2 Þ ðm þ aÞI, dt N N dA I A ¼ ðl1 þ el2 ÞbS ð1q1 Þ þ q2 m þ ðrl3 þ sl4 aÞA: dt N N
l2 . The proportion of time safe sex prevents the transmission of gonorrhea is ð1eÞ. Thus for the proportion of the time an individual participates in safe sex, l2 is multiplied by the reduction factor, e. These parameters combine to become a total reduction factor on the force of infection, i.e. ðl1 þ el2 Þb. Individuals are recruited into the I class from the S class at a rate ðl1 þ el2 ÞbSðq1 I=N þ ð1q2 ÞA=NÞ, where q1 and (1 q2) are the proportion of individuals that become symptomatically infected from contact with symptomatically infected individuals and asymptomatically infected individuals, respectively. Similarly individuals are recruited into the A class from the S class at a rate ðl1 þ el2 ÞbSðð1q1 ÞI=N þ q2 A=NÞ, similarly where q2 and (1 q1) are the proportion of individuals that become asymptomatically infected from contact with asymptomatically infected individuals and symptomatically infected individuals, respectively. Individuals from the I class reenter the susceptible class due to treatment at a rate a. Since individuals in the A class do not know they have gonorrhea the average duration of infection of an asymptomatic person is assumed to be longer than that of a symptomatic individual. We represent this increase in infectious period via rl3 þ sl4 A ½0,1. The proportion of people who engage in low levels of self-awareness, do not get tested often, is denoted by l3 while the population of those who have a high level of selfawareness are represented by l4 . The parameters r and s are extensions to the symptomatic treatment rate a caused by low and high levels of self-awareness, respectively, and by the nature of their meanings r o s. Thus reentry into the susceptible population from the asymptomatic class occurs at the rate ðrl3 þ sl4 ÞaA. The model assumes a constant population size, N, with constant recruitment and removal, m. This limits our ability to extend the time scale of this work but allows for analysis to be carried out that would not be possible. For this reason we are choosing college bound (a time of high-sexual activity and educational influence) as a further targeting of the population. The reasonable time scale for the model is thus over 4 years, where it is assumed that model parameters cannot change and the population will hold relatively constant. We also do not have memory in our population. That is, once diagnosed with asymptomatic gonorrhea an individual is not more likely to get frequent testing (Hethcote and Yorke, 1984). Here, new individuals will enter the system as old ones leave making the absence of memory an acceptable assumption. We present control here as a meaningful separation of the more classical parameters, but these are still averaged measures of behavior. Thus, we cannot make any tracking of how often a specific individual has been diagnosed asymptomatic or has become infected at all. We also do not look at dynamic changes in sexual mixing through the use of evolving contact structures. A few proposed generalizations will be made in the end of the paper as veins of future work but the thought governing this work is to start with as little detail as seems relevant.
2.2. Parameter estimation
ð1Þ
The rate at which the susceptible population is lost to infection is bSðI þ AÞ=N, where b is the effective contact rate. We introduce control into the system via l1 , l2 , and e. The proportion of the population participating in non-safe sex, or equivalently the proportion of time an individual participates in non-safe sex, is l1 and the proportion of the population participating in safe sex is
The system we wished to model was a college attending (18–22 years old), MSM population thus m is taken to equal 14 years. This group was chosen due to the prevalence of both highly sexually active individuals and the ability for policy changes to potentially reach the entire population. The behavioral parameters desired for the model were not found from a single source, or over a single demographic across many sources, and thus we have put together several sources’ worth of parameter estimates as an approximation to the potential behavior of the study group. The effective contact rate is the product of the number of risky
B.R. Morin et al. / Journal of Theoretical Biology 267 (2010) 35–40
contacts per year and the proportion of these that lead to infection. Based on the information given by NIH (2000), roughly 50% of MSM contacts with an infected individual leads to a new infection in the susceptible partner. Making an estimate of 100 risky contacts per year we arrive at b ¼ 50 infectious contacts per year. To determine the amount of time spent practicing safe versus risky behavior, according to a study of an MSM population by Shlay et al. (2004), 25.6% of participants report consistent condom usage, thus l2 ¼ 0:256 and l1 ¼ 1l2 ¼ 0:744. To determine e we take into account that although condoms are 97% effective at preventing gonorrhea infection with perfect use, many uses are imperfect due to slippage, breaking, etc. According to Shlay et al. and Stone et al. (1989), there is a 16.6–17.3% usage failure of condoms in a MSM population. We chose to let e 0:173. According to NIH (2000) among men, symptoms 4 6 typically appear within 4–6 days, 365 2365 years, following 7 infection while treatment duration is 7 days, 365 years, thus we take a 365 12 years. A study by Fortenberry et al. showed that roughly 50% of men got tested for STDs, including gonorrhea, annually (Fortenberry et al., 2002). This study was done without differentiating between hetero- or homo-sexual males and the age range was larger than that of your typical college bound of 18–22, with about 50% of their sample in our age range. This therefore serves as a poor estimate, but an estimate none the less, and is probably an overestimate of safe behavior. Using this information we let l3 ¼ l4 ¼ 0:5. For the values of r and s we will use testing once every year and once every 3 months, the CDC recommendation for highly active MSM (CDC, 2007). Thus we 365 12 1 have r a ¼ 365 372 years and s a ¼ 99 years, or r ¼ 372 ¼ 31 and 12 4 s ¼ 99 ¼ 33. Making r represent a year as asymptomatic is probably an underestimate for this group, but it is our hope that this models some form of ‘‘college’’ pressure: student health efforts, peer pressure, partner desire. The probabilities q1 and q2 are two difficult parameters to estimate. If data were collected that found the probability of exhibiting a certain type of infection given a successfully infectious contact with an infectious individual (not necessarily of the same type) then we would immediately have these parameters. However, what we generally have is some presence proportion (i.e. a percentage of the infected population who is of a particular type). Thus, the probabilities would solve the necessary relationship in assumed endemic populations. According to Bozicevic et al. (2006) 29.9% of infected individuals are asymptomatic. Thus we would have to solve 0:299I1 ¼ 0:701A1 for q1 and q2. This is daunting since the expression for the endemic equilibrium is very large and the result would not necessarily be an invertible function (i.e. more than one pair of the probabilities would give the same relation). A numerical investigation, using the nominal values for every other parameter and varying q1 and q2, provides a few values for the respective probabilities that are near the 29.9% estimate. Choosingq1 ¼0.965 and q2 ¼ 0.034 we get an endemic percentage of asymptomatic individuals of 29.85%.
3. Analysis 3.1. Global stability analysis With the model constructed we wish to do a full analysis on the qualitative behavior of the system. To that end the main focus of this section is to prove Theorem 1. The discussion proving this will involve a treatment of the stability of the disease free equilibrium, conditions on the number of endemic equilibrium that may exist, and a preclusion of closed orbits which will make all stability arguments global.
37
Theorem 1 (Global stability). System (1) has two fixed points: a disease free and an endemic equilibrium. The disease free equilibrium is globally stable when the control reproductive number is less than one and unstable when the control reproductive number is greater than one. The endemic equilibrium does not exist biologically when the control reproductive number is less than one and is globally stable when the control reproductive number is greater than one. To start we wish to formulate System (1) into a slightly easier form. Since S(t)+ I(t)+ A(t)¼N we may eliminate one state variable for the purpose of analysis. We may also rescale in both state and time to reduce the overall system. Using xðtÞ ¼ I=N, yðtÞ ¼ A=N, and t ¼ t obð1xyÞ, where o ¼ l1 þ el2 and p ¼ rl3 þ sl4 , we arrive at dx q1 x , ¼ q1 x þ ð1q2 Þy dt RII ð1xyÞ dy q2 y , ¼ ð1q1 Þx þ q2 y dt RAA ð1xyÞ
ð2Þ
where RII ¼ boq1 =ðm þ aÞ and RAA ¼ boq2 =ðm þ paÞ. Since each state variable S, I and A are positive we have that x þ y r1. Thus the rescaling is positive invariant. There is a disease free equilibrium that always exists, DFE :¼(x*, y*)¼ (0,0), implying S* ¼N, I* ¼0 and A* ¼0. To determine the local stability of the DFE the system is linearized about (x*, y*) resulting in the following: 0 B Jðx ,y Þ ¼ B @
q1 q1 RII ð1q1 Þ
ð1q2 Þ
1
C C q2 A: q2 RAA
ð3Þ
The characteristic polynomial of the above Jacobian is q q q q l2 q1 þq2 1 2 l þ q1 1 q2 2 ð1q1 Þð1q2 Þ, RII RAA RII RAA 2
which is in the form l bl þ c. It may be shown that conditions for the determinant of the Jacobian to be positive, c 4 0, are identical to RE, the basic control number, being less than one. First, we use the next generation operator to compute the number of secondary infections a typical infectious individual creates in a completely susceptible population (van den Driessche and James, 2002). This method essentially creates two vectors F and V. The vector F contains any rates in the infected classes’ ODEs that constitute recruitment from a non-infected class. The vector V contains the negative of all other rates except for the rates in noninfected classes that stand for recruitment into the infected classes. This is more easily done with Eqs. (1) than the rescaled system in Eqs. (2). This gives us 0
1 0 B I A C B obS q1 þ ð1q2 Þ C B N N C F ¼B , C B I A C @ A obS ð1q1 Þ þq2 N N 0 B V ¼@
mðNSÞaIpaA
1
ðm þ aÞI
C A:
ðm þpaÞA
ð4Þ
ð5Þ
We then define F and V as the Jacobians of F and V, respectively. Specifically these are matrices where the (i, j)th element is the partial derivative of the ith term of the vector with respect to the jth state variable. The spectral radius, or dominant eigenvalue, of FV 1 is then our basic reproductive number. The numerator of
38
B.R. Morin et al. / Journal of Theoretical Biology 267 (2010) 35–40
each term in Eq. (6) is from F and each denominator is from V 1 : 2 3 boq1 boð1q2 Þ " # 6 aþm RII RAI m þpa 7 6 7 1 FV ¼ 6 ¼ : ð6Þ boq2 7 RIA RAA 4 boð1q1 Þ 5 aþm m þpa Interestingly, the terms in Eq. (6) are each control reproductive numbers in and by themselves. The term Rij is the average number of individuals recruited into class j from a typical member in class i per unit time. Thus the control reproductive number of the entire system is a function of each of these individual recruitment thresholds. The spectral radius of Eq. (6) is qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi RII þ RAA þ ðRII RAA Þ2 þ 4RIA RAI RE ¼ : ð7Þ 2
Corollary 1. The condition RE o1 is identical to c 40 for the 2 characteristic polynomial l bl þ c for System (2). Proof. We begin with the condition for RE o1: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi RII þ RAA þ ðRII RAA Þ2 þ 4RIA RAI o 1, 2
Corollary 2. If q1 þ q2 o1 then RE o 13C o 0 and RE 413C 40. If q1 þq2 41 then RE o 13C 4 0 and RE 413C o0. Proof. If q1 þ q2 o 1 then C and c have differing sign and thus if c 40 then C o 0 and their relationships to RE are the opposite of one another. If q1 þ q2 41 then C and c have the same sign and their relationships to RE are identical. & The equation f(z) is a quadratic and thus may have 0, 1, or 2 roots in (0,1). The relative signs of f(0) and f(1) will allow us to determine under what conditions f(z) has a particular number of roots in the unit interval. Consider f(0) ¼C and f(1)¼1+ C B. We already have that the sign of C may be viewed as being dependent on the magnitude of RE. We also have that 1þ CB ¼
q1 q2 1 , 1q1 q2 RII RAA
whose sign depends on 1 q1 q2. Thus in order to study the zeros of f(z)¼ z2 Bz+ C we must consider all four combinations of the sum of q1 and q2 with the magnitude of the control reproductive number.
ðRII RAA Þ2 þ4RIA RAI o ð2ðRII þ RAA ÞÞ2 ,
1. RE o 1 and q1 þ q2 o1: In this situation f ð0Þ ¼ C o0 and f ð1Þ ¼ 1 þCB o0. Thus there are no zeros of f ðzÞ A ð0,1Þ. 2. RE 4 1 and q1 þ q2 o1: Here we have that f ð0Þ 4 0 and f ð1Þ o0. Thus there is a single root for f ðzÞ A ð0,1Þ. 3. RE o 1 and q1 þ q2 4 1: In the most difficult situation we have that f ð0Þ 4 0 and f ð1Þ 4 0. The quadratic having exactly two roots occurs when (a) B=2A ð0,1Þ, and (b) B2 4C 4 0.
R2II 2RII RAA þR2AA þ4RIA RAI o 44ðRII þ RAA Þ þðRII þ RAA Þ2 , 4RII RAA þ 4RIA RAI þ4ðRII þRAA Þ o 4, RII þ RAA þ RII RAA
where c is from the characteristic polynomial. If z¼0 then x¼ y¼0, the DFE. The interest thus lies in where f(z)¼0. A relationship between C and RE can be made using the existing relationship found in Corollary 1.
1q1 q2 o1: q1 q2
If we then consider the condition for c 40: q1 q2 q1 q2 ð1q1 Þð1q2 Þ 40, RII RAA
In order for B 4 0 we require that q1 q2 ð1RII RAA ÞRII RAA ð1q1 q2 Þ 4 0, 1RII RAA 4 RII RAA
RII þ RAA 2ðq1 þ q2 1Þ o : q1 q2 RII RAA
1q1 q2 , q1 q2
The second condition yields
1q1 q2 RII þ RAA þ RII RAA o1, q1 q2 then we see Eq. (8) identical to the condition for RE o 1.
ð8Þ &
Thus we have a somewhat easier condition for stability and may define 1q1 q2 R^ E ¼ RII þ RAA þRII RAA : q1 q2 This quantity is somewhat easier to understand and its relationship to 1 is identical to RE’s. If RE 4 1 then the DFE is unstable, but we have not discussed how many equilibrium may exist. Consider z¼x+y and the fact that if dz=dt ¼ 0 and dy=dt ¼ 0 then we would be at a fixed point for System (2). Solving dz=dt ¼ 0 we get an expression for y in terms of z by noting x¼z y. Plugging this expression into dy=dt and solving the new expression for zero we get zf(z)¼0, where f(z)¼z2 Bz+C, 1 1 q1 q2 þ RAA RII , B ¼ 2þ 1q1 q2 and 1 1 ð1q1 Þð1q2 Þq1 q2 1 1 c RII RAA ¼ , C¼ 1q1 q1 1q1 q2
ðRII þRAA Þ2 4ðq1 þq2 1Þ 4 : q1 q2 RII RAA
ð9Þ
These two conditions are contradictory. To illustrate the contradiction we invoke Eq. (7) to get that RE o 1-RII þ RAA o2. Rearranging Eq. (9) results in RII þ RAA 4 ðq1 þq2 1Þ 4 , RII þ RAA q1 q2 RII RAA
ð10Þ
but 4=ðRII þ RAA Þ 4 2 which results in the contradiction. Thus there are no zeros of f ðzÞ A ð0,1Þ. 4. RE 4 1 and q1 þ q2 4 1: In this situation f ð0Þ 4 0 and f ð1Þ o0. Thus there is a single root for f ðzÞ A ð0,1Þ. Combining the arguments gives us that if RE o1 then the only solution that is biologically relevant for our system is the DFE. When RE 4 1 the two solutions, which are biologically relevant, the unstable DFE and a single endemic equilibrium, EE. The entire above argument is valid only if 1q1 q2 a 0. When 1q1 q2 ¼ 0, the system exhibits a single EE, ! RII ðRII þRAA 1Þ RAA ðRII þ RAA 1Þ ðx,yÞ ¼ , , ðRII þ RAA Þ2 ðRII þ RAA Þ2 which is only valid if RII þ RAA 1 Z 0. We wish to make assertions about the stability of the EE without having to do a linearization
B.R. Morin et al. / Journal of Theoretical Biology 267 (2010) 35–40
around the fixed point which is very term intensive. We may now disprove the existence of closed orbits in the plane and thus assert that when the EE exists, RE 4 1, it is globally stable. Corollary 3. System (1) has no closed orbits. Proof. By Dulac’s criterion, if there exists a function j A C 1 such that _ _ @ðjxÞ @ðjyÞ þ a 0, @x @y _ y_ has no closed orbits. Let j ¼ 1=xy. then the planar system x, Then _ @ðjxÞ @ q1 ð1q2 Þ q1 ¼ þ @x y x @x RII ð1xyÞy ð1q2 Þ q1 ¼ , x2 RII ð1xyÞ2 y _ @ðjyÞ @ ð1q1 Þ q2 q2 þ ¼ @y y @y x RAA ð1xyÞx ð1q1 Þ q2 ¼ : y2 RAA ð1xyÞ2 x _ _ þ @ðjyÞ=@y is Since q1 ,q2 A ½0,1 and RAA ,RII 4 0 the sum @ðjxÞ=@x always negative. Thus there are no closed orbits for System (2). Since the dynamics are identical for System (1) we have precluded limit cycles in it and have shown what was intended. & 3.2. Sensitivity of gonorrhea transmission model In a perfect world, public health initiatives would be simple, multifaceted, and have great effects on the dynamics of a disease. However, this is not always the case and moreover, economic costs have to be considered in determining which interventions to support. Another important question is whether the size of the population is important to the dynamics of the disease and intervention. In order to address these concerns, in this section we examine the sensitivity of the system to changes in the parameter values. We take two approaches, first the time-dependent sensitivity of the original system to changes in parameter values is calculated via elasticity and second we compute what changes to the parameters would be necessary to bring the control reproductive number to a value less than 1.
39
difference quotient to arrive at i @X Xðt; y^ ÞXðt; y^ Þ ðtÞ : @yi Di
ð11Þ
We may use this approximation in our computation of the elasticity of J(t) with respect to each parameter. The sensitivity of the total infected population with respect to the nominal parameters is presented in Fig. 1 with Di ¼ 0:001yi used in each of the approximations. Not surprisingly b acts as a large control on the infected population. Small changes in b will result in the largest changes to J(t) and increase the number of individuals infected. The parameter modeling non-condom use, l1 , intuitively is also very large. The remaining four control parameters, l3 , r, s, and e are all near zero and thus each will require large changes to impact the infected population’s size.
3.2.2. Reducing RE Holding all but one parameter constant we varied the free parameter until RE o 1. Each parameter was varied by 0.1% in the correct direction to reduce the infected population, analogous to reducing RE. Of these only two single effort campaigns resulted in a control reproductive number less than 1, l1 and b. Reducing l1 to roughly 0.392, condom use increases to 61% of the time, or b to 28.239, 56.5 risky contacts per year, will result in RE o 1. Each represents a fairly significant change in behavior. Even with the entire population following the CDC recommended testing frequency of once every 3 months the disease will remain prevalent. Instead if one considers a multifaceted campaign to disease reduction we may vary each parameter at once, by 1/10%, until RE o1. This will represent smaller changes to each parameter. However, we fixed s at the CDC recommendation to measure its effectiveness. If the control parameter vector reads ðl1 , l3 , b, r, eÞ ð0:5605,0:3767,37:67,0:0428,0:1303Þ then gonorrhea will no longer exist within the population. This translates to 43.95% of the population using a condom, up from 25.6%, 62.3% of the population using the CDC’s testing frequency recommendation, reduce the number of risky a contacts per year from 100 to 75.34, even if you do not follow the CDC recommendation for testing get tested once every 9 months and finally reduce condom failure rate from 17.3% to 13%.
3.2.1. Sensitivity of I(t)+ A(t) In order to discuss the importance of individual parameters one must investigate how their value affects the state variables over time. We will consider the concept of elasticity. Formally, one may define elasticity of a function Xðt; yÞ with respect to a parameter yi A y as Eyi ðtÞ ¼
yi Xðt; yÞ
@Xðt; yÞ : @yi
This measures the ratio of a percent change in a parameter to that of the function. This gives us a unit-less and scaled method with which to compare each parameter’s affect on the solution J(t):¼I(t)+ A(t). However, we do not have a closed form for J(t) and thus we must make an approximation to both its value and its partial derivative. In general, consider dX=dt ¼ f ðt,X; yÞ where y is a parameter vector. Now consider the vector of nominal parameter values, y^ , and a small perturbation, Di , of the ith element, y^ i , and call this i i new parameter vector y^ . If what we are interested in is @X=@y near the nominal value then we could do the following. i Numerically find the solutions Xðt; y^ Þ and Xðt; y^ Þ then take the
Fig. 1. Sensitivity of the infected population, regardless of symptom expression, with respect to each model parameter using the nominal values ðm, l1 , l3 , b, a, r, s,q1 ,q2 , eÞ ¼ ð1=4,0:744,0:5,50,365=12,1=31,4=33,0:965,0:034,0:173Þ and N¼ 10 000 with S(0) ¼ 99 999, I(0) ¼0, and A(0)¼ 1.
40
B.R. Morin et al. / Journal of Theoretical Biology 267 (2010) 35–40
4. Conclusion and discussion There are a few shortcomings in this model that need to be addressed. The removal of a constant population size could potentially remove the monotonicity of the system and thus would remove the hope for any kind of qualitative global stability analysis via criteria as supplied herein. It would fall upon the modeler to do a more complex stability analysis where cycles may develop. It is unclear to the authors here if gonorrhea prevalence in this type of population merits the need for cyclic dynamics. In this paper we have considered the average, safe and risky behavior of the population. While this has given insights into how education affects disease transmission, it does not explicitly treat the heterogeneity present in such populations. This is a step forward from more classical models where the individual affects were not treated explicitly, but more work may be done for, if not analytical then numerical, results to suggest policy. Several treatments could simply increase the number of compartments while still not dealing with individual level concerns. Specifically, we may take each of the S, I, and A classes and divide them into those who are ‘‘safe’’ and those that are ‘‘risky’’ with respect to each behavior resulting in a system of 12 equations (i.e. susceptible population who are risky with respect to protection and self-awareness, risky with only one behavior or risky with neither) with flow between like-state groups of different ‘‘risk’’ types. Work done in Castillo-Chavez (1989), Blythe et al. (1991), Blythe and Castillo-Chavez (1989) and Morin et al. (accepted for publication) has also shown the importance of affinity-based mixing while others have turned to simulation to look at prevention strategies on networks (see Kemper, 1978). Here we have proportionate mixing which should not be accurate given the pain of symptoms and the mixing assumed between risky and not risky individuals. Furthermore, there is this intuitive desire to model stubbornness in a system such as this. Two difficulties with educational measures is how broad of outreach one should invest in, how many people get to see the information, and how intense the education should be. The trade offs between should we remind everyone a little or remind a few constantly can be addressed with more game theoretic or individual based approaches (Shim and Galvani, 2009; Shim et al., 2010). In this paper, a simple model of homosexual gonorrhea transmission between men is analyzed and the parameter sensitivities are determined. A caveat to the sensitivity portion of analysis is that the results are based off of a set of imperfect parameters. The populations studied in Shlay et al. (2004), Stone et al. (1989) and Fortenberry et al. (2002) are non-identical in age, region or sexuality and thus the parameters chosen are done as a starting point. The Fortenberry study covered men from ages 14 to 59, with mean age 24 with no distinction made for sexuality, from seven different US cities. The Shlay study covered MSM subjects with mean age 33.4 from a single clinic in Denver. The Stone paper includes a sample of homosexual men where half of which are between the ages of 18 and 30 participating in a multicity HIV study. Work could be put to good use if it sampled all the necessary data from several groups of men to determine if the behaviors are similar across region for the ages 18–22. In Section 3.2, it is shown that disease education is useful in disease control, however, education should be multifaceted in order to reduce the need of changing any one parameter by too much. In detail, the system is most sensitive to b (the number of risky contacts per year), and this parameter, as well as l1 , may be changed all on its own enough to eradicate gonorrhea. These changes, however, appear to represent too drastic a change to be realistic. The sensitivity analysis has also revealed that the two methods of prevention, frequent condom use and the number of partners taken, have not only the strongest effect on the infectious population but the methods of self awareness, the proportion who participate in frequent testing and the frequency of this testing,
have a sensitivity more than 15 times smaller. This supports a ‘‘prevention first’’ intuition to reducing the prevalence of the disease. We calculated time-dependent sensitivities, and although gonorrhea is best viewed as an endemic disease there are fouryear residence systems, like a high school, where the introduction of a single infected individual is feasible. In such situations we have shown that within the first 6 months prevention measures serve the goal of disease eradication best since after a few years the effect is approximately one seventh what it was at its highest.
Acknowledgements This project has been partially supported by grants from the National Science Foundation (NSF - Grant DMPS-0838704), the National Security Agency (NSA - Grant H98230-09-1-0104), the Alfred P. Sloan Foundation and the Office of the Provost of Arizona State University. The authors would like to thank the suggestions and effort put forth by Carlos Castillo-Chavez, Christopher Kribs-Zaleta and Baojun Song in order to accomplish the analysis that has been done to the model. References Anderson, R.M., 1982. Population Dynamics of Infectious Diseases: Theory and Applications. Chapman and Hall, London, New York. Anderson, R.M., 1991. Infectious Diseases of Humans. Oxford Science Publications, Great Britain. Bailey, N.T.J., 1975. The Mathematical Theory of Infectious Diseases and its Applications. Griffin, London. Blythe, S.P., Castillo-Chavez, C., 1989. Like-with-like preference and sexual mixing models. Math. Biosci. 96, 221–238. Blythe, S.P., Castillo-Chavez, C., Palmer, J., Cheng, M., 1991. Towards a unified theory of mixing and pair formation. Math. Biosci. 107, 379–405. Bozicevic, I., et al., 2006. Epidemiological correlates of asymptomatic gonorrhea. Sexually Transmitted Dis. 33 (5), 289–295. Castillo-Chavez, C. (Ed.), 1989. Mathematical and Statistical Approaches to AIDS Epidemiology, Lecture Notes in Biomathematics, vol. 83, Springer-Verlag, New York. Centers for Disease Control, CDC Fact Sheet: Gonorrhea. December 2007. Dietz, K., 1988. On the transmission dynamics of HIV. Math. Biosci. 90, 397–414. Dietz, K., Hadeler, K.P., 1988. Epidemiological models for sexually transmitted diseases. J. Math. Biol. 26, 1–25. Fortenberry, et al., 2002. Relationships of stigma and shame to gonorrhea and HIV screening. Am. J. Public Health 92 (3). Hadeler, K.P., 1989. Modeling AIDS in structured populations. In: 47th Session of the International Statistical Institute, Paris, August/September, Conference Proceedings C1-2, pp. 83–99. Hethcote, H.W., Yorke, J.A., 1984. Gonorrhea Transmission Dynamics and Control, Lecture Notes in Biomathematics, vol. 56. Springer-Verlag, New York. Hethcote, H.W., Van Ark, J.W., 1992. Modeling HIV Transmission and AIDS in the United States, Lecture Notes in Biomathematics, vol. 95. Springer-Verlag, New York. Kemper, J., 1978. The effect of asymptomatic attacks on the spread of infectious disease: a deterministic model? Bull. Math. Biol. 4 (I) 707–718. Li, J., Ma, Castillo-Chavez, C., 2003. Coexistence of pathogens in sexuallytransmitted disease models. J. Math. Biol. 47, 547–568. doi:10.1007/s00285003-0235-5. Morin, B., et al., accepted for publication. Notes from the heterogeneous: a few observations on the implications and necessity of affinity. J. Biol. Dyn. National Institute of Allergy and Infectious Diseases, National Institutes of Health, Department of Health and Human Services. Workshop Summary: Scientific Evidence on Condom Effectiveness for Sexually Transmitted Disease (STD) Prevention, June 2000. Shim, E., Galvani, A., 2009. Impact of transmission dynamics on the costeffectiveness of rotavirus vaccination. Vaccine 27, 4025–4030. Shim, E., Chapman, G., Galvani, A., 2010. Decision making with regard to antiviral intervention during an influenza pandemic. Medical Decision Making 30 (4), E64–E81. Shlay, J., et al., 2004. Comparison of sexually transmitted disease prevalence by reported condom use: errors among consistent condom users seen at an urban sexually transmitted disease clinic. Sexually Transmitted Dis. 31 (9), 526–532. Stone, E., et al., 1989. Correlates of condom failure in a sexually active cohort of men who have sex with men. J. Acquired Immune Defic. Syndr. Human Retrovirology 20 (5), 495–501. van den Driessche, P., Watmough, J., 2002. Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 180, 29–48.