THE PROVINCE OF British Columbia, Canada, has experienced a sustained outbreak of heterosexually transmitted infectious syphilis since 1997, focused in a 10-block area of Vancouver's Downtown Eastside, where drugs, prostitution, HIV, AIDS, and hepatitis are common. Annual new syphilis cases increased from 10 to 17 during the period of 1992 to 1996 and to 126 by 1999 (Figure 1). The rate of cases per 100,000 population increased from less than 0.5 to more than 3.0. In order to control the outbreak, a mass treatment intervention was attempted in early 2000, which resulted in an immediate drop in the number of cases. Despite the initial success of the program in reducing the number of cases by 40%, the number of cases of syphilis reported in the year 2001 was unexpectedly higher (177) than in the period before the intervention. We queried whether the observed increase was the result of mass treatment or whether it was the natural continuation of the previous upward trend that had begun in 1997.

To answer this question, we excluded the data for the years 2000 and 2001 and then projected the number of cases in 2001 by fitting several curves to the pretreatment data, using different techniques (only a few curves are shown in Figure 1). The analysis shows that the projected values for 2001 with any technique are significantly different from the observed value (*P* < 0.001 for all scenarios). Several factors may have contributed to this difference, such as imported cases and heightened public and physician awareness, resulting in better case-finding. Awareness has been high throughout the outbreak, however, and we found no evidence that imported cases played a greater role in 2001. We also found no other reasons for the high number of cases, except for the mass treatment program itself. If syphilis mass treatment results in a rebound in cases beyond what would have occurred naturally, this phenomenon would be important to document. We decided to investigate this further.

To assess the potential impact of mass treatment on syphilis dynamics, we developed a compartmental model that incorporates demographic and epidemiologic aspects as well as the sexual activity level of the population under study. Models developed earlier to address this issue are simple either in form, i.e., they include few compartments, ^{1,2} or in demographic structure, i.e., they are not stratified by age and sexual activity level together. ^{3} It is worth mentioning that the complexity of a model comes at the cost of introducing more parameters in the system for which there are no reliable estimates. Uncertainty in parameter values is particularly high for behaviorally related issues, such as the mixing patterns of different subgroups within the population.

Despite these uncertainties, we were able to make general conclusions in terms of the impact of an intervention program on syphilis dynamics. In the next section we introduce the mathematical model to describe the syphilis transmission dynamics. In section three, we present the results of the simulations, using demographic data corresponding to the sexually active population in Vancouver. Finally, in the last section, we present a discussion of model conclusions.

#### Methods

The compartmental model we used is shown schematically in Figure 2. The figure shows the transfer diagram within a sexually active population. In the beginning everyone is assumed to be susceptible and therefore belongs to the susceptible class, *S*. After being exposed to an infectious person there is a probability β that a susceptible individual acquires infection during a sexual partnership. At this moment the person is not yet infectious, so that the individual moves to the exposed class, *E*, at the rate λ, the force of infection. After this incubation period the infected individual becomes infectious and moves toward the infectious class, *I**1*, at the rate ς, which is inversely proportional to the incubation period.

This class (*I**1*) represents the population with primary syphilis. Ideally, all members of this class should have syphilis diagnosed, but in reality diagnosed cases of primary syphilis account for only a fraction of this population. In the absence of treatment, the disease will advance to the next stage, secondary syphilis. The individual moves to the *I**2* class, at the rate of γ. The person in *I**2* is more likely to pass on infection than a person in *I**1*. The patient then is cleared of infection and after spending some time in the temporary immune class, *P*, moves back to the susceptible class or enters a latency period. Eventually, tertiary disease can emerge after an extended period of latency. The *I**3* class includes the fraction of the population in whom secondary syphilis but not the underlying infection has resolved.

Garnett et al ^{3} divide this population into smaller subgroups according to the replication rate of *Treponema pallidum*. For the purposes of our study, those subgroups can be merged into one compartment, *I**3*. Since our goal is to assess the impact of mass treatment on syphilis dynamics, we delineated the class, *T*, corresponding to those who receive treatment. Since mass treatment may affect the entire population, we sketched a link from all compartments to the *T* class to represent the flow of the treated people (dotted lines). The immunity acquired after treatment wanes with time, and people move from the *T* class to the susceptible class, *S*. The inclusion of the *P* and *T* classes enables us to assess the impact of temporary immunity after natural clearance of the infection and after treatment, respectively, a notion about which there is no consensus in the literature.

We also divided the population into groups according to age and sexual activity levels, for each of which we assigned a flow diagram like the one indicated in Figure 2. For each gender, we have 5 age groups (15–19,20–29,30–39,40–49, and 50–59 years), and 3 sexual activity levels (active, very active, and core group). The mixing pattern between these subgroups is important.

To calculate the mixing matrix between a subgroup of one gender and all subgroups of the opposite sex, we took the approach suggested by Anderson et al. ^{4} Whereas Anderson and colleagues divided the population into different age groups only, we have incorporated both the age groups and sexual activity levels. We now describe how we calculated the mixing matrix for the heterosexually active population and how we calculated the force of infection for disease transmission.

We normalized the total population so that it is equal to unity. The Greater Vancouver Regional District (GVRD) population is divided fairly evenly between opposite sexes (679,734 males versus 675,490 females, i.e., 1:1.006 ratio, in the year 2001) so we represent each gender as 0.5. We represent in Appendix 1 the elements of the preference matrix for males (females) as *p**isjν* (*p**′isjν*), which defines the probability that a person of one gender in age class *i* and sexual activity level *s* chooses a sexual partner of the opposite sex in age class *j* and sexual activity level *v*. After the pattern of Anderson, ^{4} we define two special mixing patterns: within-group mixing and proportional mixing. We next build a third and more general mixing pattern, preferred mixing, that can create a spectrum of different patterns by changing two parameters introduced in this scheme.

In the within-group mixing, males of age *i* and sexual activity level *s* tend to have sex with partners of the same age and the same sexual activity level. In proportional mixing, a member of each gender chooses partners of the opposite gender in a specific subgroup in proportion to the representation of that subgroup in the total population of the opposite sex. Finally, in the preferred mixing, we consider the preference of older men for younger women, against a background of proportional mixing and within group mixing.

We introduce the parameters Ψ and Φ (0 ≤ Ψ, Φ ≤ 1, such that 0 ≤ 1 − Ψ − Φ ≤ 1) to create the desired combination of those three patterns. Ψ is the weight factor of proportional mixing in the overall pattern, Φ the weight factor for within-group mixing, and 1 − Ψ − Φ the weight factor for the preference of males for younger females within the same sexual activity level (see Appendix 1 for details). In this scenario the youngest male age group has proportional mixing with female subgroups. We then use the estimated forces of infection to solve the dynamic equations for the syphilis transmission, which are shown in Appendix 2.

##### Parameter Values

For the simulations presented here, we chose five distinct age groups: 15–19, 20–29, 30–39, 40–49, and 50–59 years old. The last age group does not have significant impact on the results, but it was included in the model because some cases were observed at this age group among males. We also chose three different sexual activity levels: active, very active, and core group. This stratification of sexual activity has been widely used in the literature. ^{1,5}

These groups are different in terms of the number of partners that they change per unit of time (e.g., per year). The active group corresponds to those who have only one or a few new partners per year; the very active group corresponds to those who have more than a few new partners per year; and the core group includes the most sexually active individuals among these three groups. The core group includes sex workers and other highly sexually active persons. For most of the parameters, we chose the baseline values from Garnett et al, ^{3} whose study report includes an elaborate literature review and parameter estimation.

Table 1 shows the baseline parameters used in this study. We set the population of each age group in such a way that it mimics the Vancouver demographic distribution. The number of partners that we chose as baseline parameters is shown in Table 2.

Table 1 Image Tools |
Table 2 Image Tools |

We placed 30% of the total population in the sexually active group, while 5% and 1% were placed in the very active and core groups, respectively. The rest of the population remains monogamous and does not contribute to the transmission of the disease. The fraction of different sexual activity levels for each age category is also shown in Table 2. In the following section, we do a sensitivity analysis by varying different parameters from their baseline values. We show that despite uncertainty in these values, general conclusions are obtained.

We also set the duration of treatment so that it mimics the mass treatment program implemented in the Vancouver Downtown Eastside area. ^{6} In that program, two doses of azithromycin were delivered to at-risk residents in the defined geographic area over two 6-day periods separated by 4 weeks. The failure of mass treatment and the substantial increase in syphilis cases in the year 2001 might be explained by the fact that two doses of antibiotics were not enough to effectively control the transmission. In order to test this hypothesis, we added a third 6-day delivery of mass treatment in our simulation, 4 weeks after the second dose. The results presented in the next section correspond to this three-dose program, although the conclusions remained valid in our simulations for a two-dose program as well.

#### Results

In this section, we present the results of the model simulation as we varied the parameter values. Figure 3 shows the time variation of the total number of primary and secondary cases. The same pattern can be seen for each stage of the disease and each gender, separately. After a rapid initial increase as a result of the introduction of a few infected cases into a totally susceptible population, the curve reaches an endemic equilibrium. After the implementation of mass treatment, this equilibrium is disturbed, and the curve moves down after the delivery of each dose (Figure 3B).

The important feature of this curve is that as soon as the mass treatment is terminated, a rebound takes place. Despite expectations for a sharp decline after mass treatment, the number of new syphilis cases surpasses the previous endemic equilibrium level and settles to an endemic equilibrium level, as if no treatment program had taken place. The rebound predicted by simulation is consistent with the pattern observed in cases in Vancouver after the mass treatment initiative. The extent of the rebound that takes place depends on the parameter values. For some parameters such as the duration of acquired temporary immunity after infection or after treatment, we have no clear estimate. Nevertheless, for any set of parameter values the simulated curve results in such a rebound.

Therefore, despite uncertainty in the estimation of some parameter values, the model reflects the observed data over a range of parameter estimates. Similar rebound behavior has also been reported by Korenromp et al, ^{7} who used a stochastic model to study the transmission of STDs in a rural African population. However, they did not include migration, and their assumptions on partner preferences are different from ours.

Figure 4 shows the simulated time variation of the fraction of the susceptible population. While the population of the active and very active subgroups (any age) remains almost entirely susceptible, a significant proportion of the core group is distributed among other compartments. This group is also the main cause for the rebound in cases after mass treatment. Because of mass treatment, many of those belonging to the core group who were otherwise exposed or infected cases (symptomatic or asymptomatic) return to the susceptible pool. This phenomenon in turn leads to an outbreak-like increase in the number of new infected cases, not only among the core group but also among the very active group in contact with core group members.

Simulations also suggest that if the entry (exit) rate to (from) the population was zero (closed system), under certain conditions (e.g., a repeated very-high-coverage mass treatment) and with assumed postintervention immunity for a period of time (θ ≠ 0), this approach could effectively control the transmission of syphilis (Figure 5A). For an open system (entry/exit rate ≠ 0), rebound is inevitable. The higher the entry/exit rate, the earlier the rebound would take place (Figure 5B). It comes as no surprise, therefore, that we observed a rebound in syphilis cases after 1 year. According to the municipal statistics, the migration rate within Downtown Eastside for the period of 1991 to 1996 was an astounding 68.9% (11.48% per year), while the change in population size for the same period was 2.4% (0.4% per year). This represents a high degree of mobility for this population. Since this population is involved in drugs, sex, and crime, the real moving rate may be even higher.

Figure 6A shows that independent of the mass treatment coverage rate, rebound will occur after mass treatment, as long as infected cases remain in the population and the entry rate to the population is not zero. In fact, according to our model, higher coverage would initially decrease the incidence of new cases to a lower level but cause a larger outbreak, later. This is in agreement with the previous hypothesis that mass treatment, while decreasing the number of infected cases, will increase the size of the susceptible pool. Since the disease is not eliminated, there is a greater likelihood that susceptibles will contract the disease for a second time, thereby temporarily increasing the incidence in a rebound fashion.

Figure 6B demonstrates that the probability of transmission during each sexual partnership will not alter the simulation results significantly, even if we replace the values for β_{1} and β_{2} from those indicated in Table 1 to those proposed by Garnett et al. ^{3} Numerical results also show that while the number of partnerships per year for the core group has an impact on the length of effectiveness of mass treatment, it will not change the overall pattern of the posttreatment rebound (Figure 7). The change in the number of sex partners for noncore groups has even less impact on disease transmission.

In these simulations sexual preferences in the sexually active population play an important role. We assumed that the mixing pattern is composed of three different preferences, namely, proportional mixing, within-group mixing, and preferred mixing—that is, men prefer sex with women at a younger age group. In our model, these three preferences can be mixed at different proportions, resulting in a spectrum of patterns. The real mixing pattern is closer to one or a few combinations within that spectrum, and the exact proportion in unknown. Despite this uncertainty, one can still draw general conclusions about the effectiveness of a mass treatment program for the entire spectrum.

To investigate mixing patterns, we introduced two new parameters, Ψ and Φ, as the fraction of the male population who prefer proportional mixing and within-group mixing, respectively. Therefore, 1 − Ψ − Φ would be the fraction of the male population that prefers to have sex with younger women. The female mixing pattern can then be calculated with Equations 4 and 5 (Appendix 1). We varied these two parameters in order to cover the entire spectrum. We chose Ψ = 0.2 and Φ = 0.3 as the baseline value for these parameters and then compared all other combinations with these baseline values. The choice of baseline values was completely arbitrary, and we by no means argue that the baseline pattern mimics the real mixing pattern in the Vancouver Downtown Eastside.

Equation A4 Image Tools |
Equation A5 Image Tools |

Figure 8A shows that if Ψ is varied, while Φ = 0.3 is kept constant, the rebound takes place faster, as the preference of men for younger females becomes increasingly dominant. Should the preference for younger females be absent, the rebound would take place at a much slower pace than we observed in the 2001 data. The rebound for this particular case occurs outside the range of the time axis in Figure 8A and was not shown in the figure. The reason for this variation is that, by adding the preference toward younger females, we are linking the two groups of most active men (20–29 and 30–39 years) to the two groups of most active women (15–19 and 20–29 years). Comparison with real data suggests that the preference toward younger females is an important behavioral pattern in shaping the sexual network in the target population.

We now repeat the same task but keep Ψ = 0.2 and vary Φ (Figure 8B). The results show that the greater the preference of men for younger females, the earlier the rebound takes place. If there was no such preference, the transmission would be curbed within each subgroup, which in turn would result in a lower incidence level, even during the before-treatment equilibrium. The weak coupling between those groups due to proportional mixing prevents each subgroup from becoming completely independent of others. If we keep the preference for younger females constant, then varying Ψ (and therefore Φ) would result in milder variation than we observed in the last two cases. As can be seen in Figure 8A, the dominance of a proportional mixing pattern results in a lower incidence. The importance of the behavioral pattern in the sexually active population is also clearly shown. To build more reliable models for public health policy, behavioral patterns need to be identified and considered. Models based on network theory may prove useful in this regard. ^{8}

We tested the response of our syphilis transfer diagram to the variation of other parameter values besides the mixing pattern. Simulations suggest that despite uncertainty about the duration of immunity after treatment, the low incidence period after mass treatment would eventually be followed by a rebound. The same argument holds if we vary the duration of possible temporary immunity after clearance of the secondary syphilis. It is important to underline that the same pattern can be seen both for ξ + α = 0 and θ = 0, which correspond to the absence of immunity after acquisition of secondary syphilis and postintervention immunity, respectively. Therefore, under the assumption that naturally acquired immunity may not exist, the results remain robust.

We varied all other parameters introduced in the transfer diagram depicted in Figure 2 and observed the same pattern presented in previous figures. The model suggests, therefore, that mass treatment does not have long-term benefit in the control of the syphilis transmission inside a sexually active population with new recruitments and where mixing is heterogeneous.

#### Discussion

In this study we used a compartmental model to explore the effectiveness of a mass treatment intervention in the control of syphilis transmission within a sexually active population during a heterosexual outbreak. The goal of our study was to assess the long-term impact of the mass treatment program that took place in Vancouver's Downtown Eastside in early 2000. Despite the initial partial success of this strategy, the syphilis case data for the year 2001 showed no lasting effect on syphilis transmission. In fact, a significant rebound occurred in 2001. Our model demonstrates that this unexpected increase may be the direct result of the mass treatment intervention in a sexually active population where new recruitments and different behavioral patterns coexist. This conclusion has important implications for the syphilis elimination effort currently under way in the United States. ^{9}

We tested many scenarios by varying all the parameters in our model. We showed that despite uncertainty in the precise value of a number of parameters within the host population, general conclusions could be drawn about the effectiveness of certain strategies. The results underline the importance of different behavioral patterns in the study of the target population. Ignoring various patterns may result in unrealistic conclusions about the elimination of syphilis within a host population. This suggests the need for more work in identifying those patterns within the target population.

Mathematical models are theoretical means of producing caricatures of real phenomena. Our purpose in this study was not to quantify the impact of mass treatment but to study such impact in a qualitative manner. Uncertainty in the precise value of parameters is an obstacle to a more quantitative approach. We had to make some assumptions in order to simplify our model. For instance, we assumed the preference for males to choose females only from the next-younger age group. This can be extended further to all younger groups. We also assumed that all new recruitments enter the susceptible class and that the entry rate for all demographic subgroups is identical.

While the demographic data used in this model represent the real age and sex profile of the Greater Vancouver area, the proportion of different sexual activity levels was chosen arbitrarily. Despite these simplifications, the model was robust in its ability to capture the key feature of the notification data and relate it to the mass treatment strategy. In conclusion, the model suggests not only that mass treatment has no benefit in long-term control of the syphilis transmission but also that it increases the number of cases soon after the intervention by returning the already infected population to the susceptible pool.