## Abstract

The success of an infectious disease to invade a population is strongly controlled by the population's specific connectivity structure. Here, a network model is presented as an aid in understanding the role of social behaviour and heterogeneous connectivity in determining the spatio-temporal patterns of disease dynamics. We explore the controversial origins of long-term recurrent oscillations believed to be characteristic of diseases that have a period of temporary immunity after infection. In particular, we focus on sexually transmitted diseases such as syphilis, where this controversy is currently under review. Although temporary immunity plays a key role, it is found that, in realistic small-world networks, the social and sexual behaviour of individuals also has a great influence in generating long-term cycles. The model generates circular waves of infection with unusual spatial dynamics that depend on focal areas that act as pacemakers in the population. Eradication of the disease can be efficiently achieved by eliminating the pacemakers with a targeted vaccination scheme. A simple difference equation model is derived, which captures the infection dynamics of the network model and gives insights into their origins and their eradication through vaccination. Illustrative videos may be found in the electronic supplementary material.

## 1. Introduction

Developing strategies for controlling the dynamics of epidemics as they spread through complex population networks is now an issue of great concern (e.g. Hethcote & Yorke 1984; Anderson & May 1991; Wasserheit & Aral 1996; Earn *et al*. 2000; Eames & Keeling 2002; Grenfell *et al*. 2002; Eubank *et al*. 2004; Keeling & Eames 2005; Riley 2007; Stone *et al*. 2007; Olinky *et al*. 2008). Future progress depends on gaining a better theoretical understanding of the spatial dynamics of disease spread, including the effects of a population's social contact structure and its network topology (see Lloyd & May 2001; Eames & Keeling 2002; Eubank *et al*. 2004; Keeling & Eames 2005; Riley 2007). Here, we show how these factors control epidemic spread and, in the process, formulate a novel aggregated targeted vaccination scheme.

We are particularly interested in diseases that confer temporary immunity to individuals after recovery from infection. This is typical for diseases such as pertussis, influenza and human respiratory syncytial virus, and some sexually transmitted diseases (STDs) such as syphilis. In terms of population dynamics, the temporary immunity is understood to give rise to recurrent epidemic oscillations (see Levin *et al*. 1989) that can have a period of several years for pertussis (Rohani *et al*. 1999) and certain strains of influenza (Ferguson *et al*. 2003), to decadal oscillations in the case of syphilis (see St Louis & Wasserheit 1998; Grassly *et al*. 2005). In simple terms, the epidemic cycles arise because of a delayed susceptible–infectious–recovered–susceptible (‘SIRS’) process in which susceptible individuals become infected, recover with temporary immunity, but then eventually return to the susceptible pool after a time delay when immunity wears off. The loss of immunity allows the susceptibles in the population to gradually build up until sufficient in number to fuel the next disease outbreak.

Here, we demonstrate our findings using two models and taking the disease of syphilis as a representative example. Syphilis is an STD caused by the bacterium *Treponema pallidum*. It has often been called ‘the great imitator’ because so many of the signs and symptoms are indistinguishable from those of other diseases. For yet unknown reasons, syphilis still poses an epidemic threat. For example, health officials reported over 36 000 cases of syphilis in 2006 in the USA. Lately, there is an even stronger desire than ever to eliminate the syphilis threat, since, even though it is a curable disease, it increases the chances of HIV infection. There is an estimated two- to fivefold increased risk of acquiring HIV if exposed to that infection when syphilis is present. Numerous myths, tales and literature are tied to syphilis; many of them are spurred by the relation of syphilis to some of the greatest geniuses and their careers, as described in Hayden (2003). Another curiosity relates to the many different names this disease has received over the years of its known existence. For example, it had a series of ‘national names’ (e.g. the ‘French disease’ in Italy and Germany, and the ‘Italian disease’ in France). These national names are due to the disease often being present among invading armies or sea crews, due to the high incidence of unprotected sexual contact with prostitutes. This may suggest that, indeed, a few ‘short-cut links’ in a small-world model represent short-route connections to far regions.

Grassly *et al*. (2005) suggested that the almost decadal oscillations seen in long-term syphilis datasets from the USA stem from the temporary immunity of this disease. Their argument is supported by the fact that gonorrhoea, which lacks temporary immunity, fails to show the same cycles in long-term datasets. This view, however, is controversial and the CDC (Douglas 2005) has countered that trends in the US syphilis epidemiology follow parallel changes in population-wide high-risk sexual behaviour (see also St Louis & Wasserheit 1998). Most probably, it is the combined presence of temporary immunity and social behaviour that is responsible for the recurrent waves of syphilis epidemics. The modelling approach described here makes it possible to investigate and assess the impact of these different but important factors. Indeed, our work demonstrates that both these factors have a crucial influence on the periodicity of syphilis outbursts.

Complex networks (or graphs) provide an important means for investigating the effects of social behaviour in population models of disease spread. Individuals are represented as nodes of a graph and edges are placed between any two individuals should there be an infection route between them (see Eames & Keeling 2002; Keeling & Eames 2005; Riley 2007). A *random* Erdos–Renyi network is formed if there is an equal probability *q* of a connection between any two individuals (Newman 2003). A *regular* and tightly clustered network structure is obtained if an individual is able to infect only his/her nearest neighbours. The random and clustered-regular graphs might be considered as two endpoints of a spectrum. Watts & Strogatz (1998) developed a scheme that allows the construction of networks that interpolate anywhere between these two endpoints. This is achieved by introducing a proportion of *p* random ‘short cuts’ between nodes in a regular graph. Only relatively few short cuts are required (*p*<0.1) to create ‘small-world’ networks that have the often realistic qualities of both a high degree of clustering and, at the same time, relatively high overall network connectivity introduced by long-range connections (i.e. via short cuts).

When considering the population dynamics of STDs, it is important to take into account that some individuals spread the disease to a much greater extent than others. In this way, social behaviour and sexual promiscuity govern the heterogeneity of the contact structure in the population. This contrasts with standard mean-field differential equation models, which are based on the assumption of a ‘randomly mixing’ population and lack a heterogeneous contact structure. However, there is no unanimous agreement on how the contact structure of the network should be fixed. Barbási & Albert (1999) and Liljeros *et al*. (2001) have argued that ‘scale-free’ networks, whose nodes have a power-law connectivity distribution, are the most appropriate for STDs. Lloyd & May (2001), on the other hand, suggest that such a formulation is unnecessarily exaggerated. Here, we follow Eames & Keeling (2002) and Keeling & Eames (2005), who suggest that small-world networks are reasonable approximations of sexual networks in restricted communities (e.g. an active neighbourhood of single people in a central city). This choice is discussed in more detail shortly.

The paper is structured as follows. First, we describe the SIRS network model and its spatio-temporal dynamics, when formulated for diseases with temporary immunity. For representative parameters, the model exhibits expanding circular waves of infection, some of which are generated by unusual ‘pacemaker centres’. The pacemakers are studied in detail, and their important role suggests a practical disease control strategy based on targeted vaccination. We show that by vaccinating or quarantining the regions around pacemaker centres, the disease can be eradicated. The vaccination scheme is tested on various more realistic modifications of the basic model. We then formulate a simple difference equation model that captures some of the main features of the heterogeneous network. The results obtained from both models support the claim that both the temporary immunity and the social structure/behaviour have equal responsibility for the existence of periodicity in the disease outbursts.

## 2. The network SIRS model

The network model is based on a two-dimensional lattice of individuals (or nodes). Each node on the lattice is occupied and connected to *k*-nearest neighbours oriented in each of four directions (north, south, east and west, with diagonal connections excluded). That is, each node is initially connected to *K*=4*k* nearest neighbours, with *k*=3 unless otherwise specified. The horizontal and the vertical edges of the lattice are ‘glued’ together creating a 2-torus. Then, with a probability *p*, each of the *K*=4*k*=12 nearest neighbours of each of the edges in the lattice is randomly rewired to an arbitrary node, not permitting duplications. These rewired connections, or short cuts (Watts & Strogatz 1998; Kuperman & Abramson 2001), may extend to far regions of the network.

The parameter *p* controls the population's connectivity structure. *p*=0 corresponds to nearest-neighbour contacts only, and where clustering is at its maximum. Small *p* in the range of 0<*p*<0.1 corresponds to a small-world network (each individual has a certain amount of nearest-neighbour contacts+a small proportion of distant contacts, short cuts). Large *p*>0.4 is qualitatively equivalent to a randomly mixing population with minimal clustering (Kuperman & Abramson 2001). As *p* is the probability of a short cut, it may be viewed as an index of population mobility. Alternatively, it may be interpreted as an index of social behaviour such as sexual promiscuity in the case of STDs, given the manner in which it controls the overall network connectivity and clustering (Watts & Strogatz 1998).

Disease dynamics follow the classical SIRS formulation (e.g. Anderson & May 1991; Kuperman & Abramson 2001; Eames & Keeling 2002; Murray 2002; Grassly *et al*. 2005), with susceptible individuals (*S*) having a probability *q* of becoming infected when linked to an infected individual (*I*). Infected individuals eventually recover from the disease after a fixed time period, *τ*
_{I}, and are conferred temporary immunity. After a time period of *τ*
_{R} time units, immunity wears off and recovered individuals (*R*) return once again to the susceptible pool (*S*) closing the SIRS loop.

This is implemented on the network using a two-dimensional cellular automata (CA), SIRS spatial model. At time *t*, an individual at the (*i*, *j*)th location of the lattice has the state *x*
_{i,j}(*t*), which is *S*, *I* or *R*. The model is based on the following transition rules:
2.1

Infections are transmitted to susceptible individuals with a probability *q*, if they are connected to an infective via a nearest neighbour or a short cut. Thus, the probability that a susceptible at place (*i*, *j*) and at a time *t* becomes infected is
, where
is the total number of infectious contacts of this individual at the time *t*; 1≤*i*≤*N* and 1≤*j*≤*M*. In (2.1), *t*
_{0} is assumed to be the time at which an individual at position (*i*, *j*) first becomes infectious, *τ*
_{I} is the number of time steps an individual remains infectious, *τ*
_{R} is the number of time steps an individual stays in the recovered class and *τ*
_{0}=*τ*
_{I}+*τ*
_{R}. The proportion of *S*(*t*), *I*(*t*) and *R*(*t*) individuals are calculated over the lattice and their dynamics are followed as a function of time.

To help fix ideas, we focus on two representative parameter settings, which are as follows.

Parameters used in a general theoretical model taken from Kuperman & Abramson (2001). The infectious period is fixed at

*τ*_{I}=4 time units and a recovery period of*τ*_{R}=9 time units (see simulations in figure 1 and in video 1 in the electronic supplementary material).Parameters associated with syphilis epidemics as based closely on the study of Grassly

*et al*. (2005; see simulations in figure 2 and in video 2 in the electronic supplementary material). In the latter case,*τ*_{I}=1 time unit, which is taken to correspond to half a year. The recovery time varies randomly and uniformly in the range*τ*_{R}∈{8, 9, 10, …, 16} time units corresponding to a period of 4–8 years of immunity.

To add realism, several other features were also incorporated. In some simulations, a birth–death process was introduced at rate *μ*. That is, a proportion of *μ* new nodes were introduced every year as well as deleted every year to mimic the birth–death process. This requires updating the network structure appropriately to allow for the new and removed nodes. In figure 2, *μ*=0.01 per half-year time step, which is equivalent to a rate of *μ*=0.02 per year. Provision was made for the possibility that a proportion *φ* of individuals fail to gain immunity after infection (similar to Grassly *et al*. 2005). Thus, *φ*=0 corresponds to all nodes passing through an SIRS loop, while *φ*=1 corresponds to all nodes exhibiting SIS dynamics (no individual can acquire immunity). Figure 2 shows a proportion of *φ*=0.2 of the population (20%) gaining no immunity per year.

The reproductive number *R*
_{0}, i.e. the average number of secondary infections resulting from a single infected individual, is of great importance and controls the dynamical behaviour of syphilis. In our simulations, *R*
_{0} is kept in the range of 1.5–3, as considered to be realistic (CDC 1997–to date; Grassly *et al*. 2005). For the network model with *K*-nearest neighbours (2.1), *R*
_{0}≈*qKτ*
_{I} (e.g. Keeling & Grenfell 2000). We chose a range of values for the probability of infection, *q*, but considered *q*≃0.3–0.4 (e.g. as in figure 2) to be a reasonable estimate, while *τ*
_{I} is approximately half a year (similarly to Grassly *et al*. 2005). Thus, in order to maintain *R*
_{0} in a realistic regime, the average number of contacts per node, *K*, must be high and lie somewhere between *K*=7 and 20. In the work presented here, we chose either *K*=8 or 12, as further supported by the analysis performed on the difference equation model (see §3). Note that, although this number of contacts might appear to be large, it is close to half of the number of the contacts assumed in mean-field models (e.g. as in Grassly *et al*. 2005). The latter models also must keep overall promiscuity rates high to achieve a realistic *R*
_{0}.

Eames & Keeling (2002) and Keeling & Eames (2005) advocate a small-world network SIR/SIS model for an STD rather than a scale-free network. These authors compare the results of two population network models: first, with a fixed number of connections per node on a regular network (*k*-regular) and, second, with a variable number of connections per node on a small-world network. It is important to emphasize that these two network models have basic differences: the *k*-regular model is a homogeneous, ordered, regular model with *only* nearest-neighbour connections (in addition to having *exactly k* links per node); the second model has a varying number of links per node and its main feature is that it is a *small-world network* (as most of the assigned links are nearest neighbours and only a few are far links, i.e. short cuts). Our network model presented here is constructed slightly differently, but is essentially similar to the Eames and Keeling network model being a small-world network and having a variable number of contacts per node.

### 2.1 Recurrent circular waves and pacemaker centres

We have found that, in the small-world regime, models of type (2.1) exhibit concentric waves, which correspond to periodic oscillations similar to those seen in the time series of syphilis data. Some of these waves are recurrent in time and spatial position, and hence give rise to unusual pacemaker centres. Indeed, numerical simulations show that, for 0.001<*p*<0.12, the model (2.1) exhibits spatial oscillations with expanding circular waves of infection travelling through the lattice (figure 1
*a*). Some of these waves are recurrent, both spatially and temporally. The latter are generated by *pacemakers* (figures 1
*b* and 2
*a*) that form at connectivity centres—localized areas denser in short cuts. The waves grow in size about the pacemakers as the infection spreads radially. When infected individuals recover, the interior of the growing wave boundary becomes a fresh pool of susceptible individuals. At the end of the cycle, a distant infectee short cuts through the network to reinfect the wave's focal pacemaker, enabling it to perpetuate. The cycle allows recurrent spatial waves to propagate with a fixed period, *T*. For the syphilis parameters (see figure legend 2), *T*≈10 years, as observed in US syphilis datasets (see CDC 1997–to date; Girvan *et al*. 2002). Similar spatial wave patterns had been observed in a number of biological contexts including epidemiology (Grenfell *et al*. 2001), ecology (Blasius *et al*. 1999), neural networks (Lewis & Rinzel 2000) and theoretical studies of excitable systems (Greenberg & Hastings 1978).

We refer to a ‘pacemaker’ as a recurring concentric wave in the two-dimensional population space. New infectees are introduced via short cuts to the centre of this wave every period of the disease dynamics, thereby allowing the pacemaker centre to persist. The physical coordinates of the pacemaker centre remain at the same location of the lattice for very long time periods, while the disease keeps returning to this centre.

These pacemakers follow a pattern-formation mechanism. We observe that there is a minimal necessary amount of short cuts needed for the creation of a pacemaker centre. In addition, within the small-world range, a pacemaker wave region has more short cuts than a non-pacemaker wave region. The initiation of a pacemaker also requires that the infection is able to spread in both the horizontal and the vertical directions. That is, the development of the infection from an initial state should progress in the two directions spanning a plane (i.e. in an X or an L shape). This condition follows the rule that heterogeneities are needed for the creation of spirals in excitable media (see Greenberg & Hastings 1978). For these reasons, having an aggregation of short cuts in a small region enhances the likelihood of creating a pacemaker. On the other hand, having too large an aggregation of short cuts in a localized area of the lattice results in the opposite effect—disease extinction will occur in this localized area due to a synchronization effect (see Earn *et al*. (1998) and below).

The patterns of waves observed from the spatio-temporal simulations of the network model have a twofold importance. First, they lead to a better understanding of the disease spread. The periodicity of the waves is an outcome of both the small-world structure and the temporary immunity: the reinfection at the centre of the wave after a certain amount of time is due to the intrinsic time delay introduced by the temporary immunity. (Without the time delay *τ*
_{R} (i.e. *τ*
_{R}=0), or if the contribution of the SIS loop is too large (*ϕ*>0.2), oscillations cannot occur.) On the other hand, that reinfection occurs at the centre of the wave is a direct outcome of the small-world structure. Second, the wave patterns have importance in that they suggest a novel spatial vaccination scheme, as described in §2.2.

### 2.2 A targeted vaccination scheme

The centrality of the pacemaker centres suggests a practical control strategy. We have found that, by vaccinating or quarantining the regions surrounding the pacemakers, the disease can be usually brought to a complete extinction within one period of the disease (in some cases, depending on the refinement of the vaccination algorithm and specific parameter values, two vaccination pulses are required). Thus, rather than the conventional scheme of immunizing approximately more than 85 per cent of the population to achieve herd immunity (Anderson & May 1991), it is only necessary to vaccinate groups enclosing the pacemakers. This requires vaccination of some 10–30% of the population (depending on the specific application and algorithm refinement; approx. 20% in most cases). Figures 1 and 2 show spatial snapshots of an infected population upon application of the vaccination scheme. In figures 1
*a*,*b* and 2
*a*, the characteristic circular waves of infection (red) are shown. Vaccination around the pacemakers leads to complete eradication of the disease. Pacemakers have such a large impact on the spatial dynamics that they are relatively easy to detect using a simple threshold algorithm that identifies recurring aggregations of infected individuals. The algorithm is based on the observation that pacemakers always appear (whenever infectives reside in their centre) approximately at the same spatial location and that they always exist near a minimum of the infectives' time series. Hence, the vaccination algorithm is based on identifying the pacemakers near a minimum of the infectives' time series and vaccinating a predefined region surrounding them (by eliminating this area from the simulation).

Once a pacemaker is identified, a small region enclosing the pacemaker is marked out and vaccinated by effectively removing these nodes from the simulation (the blue rectangles in figures 1
*c*,*d* and 2
*b*,*c*).

For the theoretical values of Kuperman & Abramson (2001; settings (i)), used in the simulations of figure 1, the ratio *τ*
_{I}/*τ*
_{R}=4/9 is relatively large, hence the pacemakers are large in size and few in number (usually 1–3). In some cases, the scheme is able to remove all pacemakers after one application of vaccinating approximately 10 per cent of the population. However, as the ‘pacemakers’ compete with one another, there are cases where only the main pacemaker(s) is removed and the secondary pacemaker(s) may appear in the next period. Eradication then requires a second application of the vaccination in the area of the remaining pacemaker(s). As shown in figure 1
*c*,*d*, in such cases, it typically requires vaccinating a total of approximately 20 per cent of all individuals over both applications to bring the disease to total extinction.

For the same model with syphilis parameters (settings (ii)), the ratio *τ*
_{I}/*τ*
_{R} varies in a simulation within the range {1/8, …, 1/16}. As this ratio is small and variable, the pacemakers generated are more numerous (usually 2–5) and smaller in size. With additional model realism (20% of the population not gaining any immunity, birth–death process), the pacemakers ‘tend to become less circular’ in shape (figure 2
*a*). As a consequence, it is necessary to vaccinate more areas, although each area is smaller in size. Nevertheless, the vaccination scheme works well and generally eradicates the disease in a single application, requiring vaccination of a relatively low percentage of individuals (figure 2). It is worth noting that the vaccination percentage required for the targeted scheme proposed here is always lower than would be required with random vaccination.

For the theoretical parameter values (i) based on Kuperman & Abramson (2001) and used in figure 1, a random vaccination scheme is successful in eradicating the disease only after vaccinating at least 80 per cent of the population, and hence is not that much of interest. However, for the syphilis parameter values, the random vaccination threshold is 43 per cent of the population (this is somewhat similar to the results presented in Zanette & Kuperman (2002) for an SIR model). A further gain can be achieved by combining random vaccination with the above targeted scheme. This may advantageously reduce the vaccination threshold to as low as 10–15% of the population for the syphilis parameters. Instead of vaccinating the entire area surrounding the pacemaker, it suffices to randomly vaccinate 60–95% of the area normally targeted, where the percentage vaccinated grows in proximity to the target's centre. Moreover, as in practice, it is difficult to obtain full coverage when vaccinating an entire population, or even a specific targeted group; combining our targeted vaccination scheme with some randomness adds realism and effectively lowers vaccination rates. For example, in figure 2, we vaccinated only an average of 86 per cent of each of the identified areas that enclose the pacemakers. In more detail, the algorithm vaccinates up to 98 per cent in the core of the pacemaker, where the high clustering of repeatedly infected individuals resides, and as low as 60 per cent (randomly chosen) in the outskirts of the pacemaker. These changes did not reduce the effectiveness of the vaccination scheme.

To pin down the key pacemakers, we use the following simple threshold algorithm. We monitored the total number of infectives as they change in time over a simulation run. After the first minimum was reached and the number of infectives just begins to rise again, we located localized clusters of infectives by searching for areas in the lattice that contain a few infectives in a predefined small square (usually areas of 20×20 lattice points). Usually, the threshold for the number of infectives clustered in one such square was set at three. The simulation was then monitored over another two epidemic cycles and the suprathreshold squares just located were checked after each minimum. If, as before, such a square contains more than a threshold number of infectives (usually, set to at least two), for all three minima, the surrounding area (of approx. 30 nodes per each side of each infective) is identified as a key pacemaker, and is vaccinated.

As one of our primary interests in the analysis of the simulation model was to understand the role of short cuts, to some extent this came at the cost of allowing for variability in node degree and studying its implications. Nevertheless, we found that variability in node degree is indeed possible and the waves and spatial pacemakers are still observed. Thus, in one scheme, targeted vaccination of pacemakers was successful even with a variation of 2–15 contacts per node (mean 11.4). However, for more excessive levels of variability, model simulations generated many small competing pacemakers rather than the two or three key pacemakers found for the nearest-neighbour lattice. A large number of pacemakers can reduce the effectiveness of the targeted vaccination scheme. In such cases, vaccinating all pacemakers sometimes resulted in vaccinating 90 per cent of the population.

### 2.3 Synchronization and disease extinction

Worthy of comment is the model's behaviour for larger values of *p*, typically *p*>0.1, outside the small-world regime, and corresponding to high levels of sexual promiscuity in the case of STDs. Counter-intuitively, the epidemic consistently dies out abruptly due to the appearance of large-scale synchronized epidemics—a well-known cause of disease extinction (e.g. Earn *et al*. 1998; Blasius *et al*. 1999; Kuperman & Abramson 2001). The synchronization manifests with the formation of large spatial aggregations of infected individuals. Upon recovery, these infectives gain temporary immunity for a lengthy time period. Thus, the areas that once contained aggregations of infectives become exhausted of susceptibles and there is no possibility for an epidemic to sustain—it soon dies out.

Although it might at first seem unusual, this has the implication that, for society at large, grand sexual promiscuity (large *p*) has the potential to eliminate STDs, such as syphilis, altogether after approximately two decades of consistent behaviour (see figures 3 and 4
*c*, as well as video 3 in the electronic supplementary material, for further illustration). The same would be true for populations in which individuals who consistently have few proximate sexual partners for the correspondingly low *R*
_{0} would ensure that the infection-free equilibrium is stable. Intriguingly, the most conducive conditions for the persistence of such STDs appear to be the small-world structure similar to the varying manifestation of sexual promiscuity seen in western society over the last centuries.

### 2.4 Bifurcations in p: the impact of short cuts

The effect of the parameter *p*, the proportion of short cuts, on the disease dynamics may be assessed from the bifurcation diagram in figure 3. The figure plots the range in the number of infectives (maximum and minimum values) for any given *p*, when the model is run using the standard syphilis parameters. The following bifurcation scenario takes place: for *p*=0, the disease goes to extinction; for small *p*<*p*
_{c}=0.001, an endemic equilibrium is reached in which there is a relatively small proportion of infectives 0<*I*
^{*}≪1. Extensive simulations show that this critical value *p*
_{c}>0 occurs at the point where the clustering coefficient begins to descend as *p* increases (see Watts & Strogatz 1998). For 0.001<*p*<0.13 (approx.), there is a limit cycle of radius depending on *p* and hence noticeable oscillations in the number of infectives; for *p*≳0.13, the disease goes to extinction owing to the synchronization effect discussed above. The ‘limit-cycle region’ is the region where pacemaker centres develop and within it lies the region where the targeted vaccination scheme works very effectively. This region is an open strip of *p* values in the small-world regime. Note that, for both *p*=0 and large *p*, the disease rapidly goes to extinction. Thus, despite the period of temporary immunity built into this model, oscillations in *I* vanish for either very small or relatively large values of *p* (i.e. outside the small-world regime). This further demonstrates the importance of social structure in the formation of the recurrent peaks observed in the real time series of US syphilis cases.

## 3. Difference equation model

We formulate the following difference equation model to help gain insights into the network model's dynamics. Let *S*
_{t}, *I*
_{t} and *R*
_{t} be the proportion of susceptible, infective and recovered individuals in a large population at time *t*, respectively. Again, let *τ*
_{I} be the time period an individual remains infectious and *τ*
_{R} the period an individual remains immune. If we assume that *τ*
_{I}=1, as in the case for syphilis, the proportion of recovered individuals can be described by the sum
, and thus
3.1
where *τ*
_{0}=*τ*
_{I}+*τ*
_{R}. The model formulation (2.1) is well known (e.g. Girvan *et al*. 2002) and is used for our new difference equation model below (which includes some spatial information) as follows. We suppose that each individual has on average *K* connections including those to its nearest neighbours as well as short cuts. For the average individual, denote by
the number of *nearest neighbours* that are infected at time *t* and denote by
the number of *short-cut links* that point to infected individuals. Then,
3.2

3.3

Now, consider two classes of nodes. The first class is distinguished by the property that its nodes have no short-cut connections. Let *q*
_{K} be the probability that such a node is infected by a nearest neighbour in a given time step. The probability can be estimated from long-term simulations. The second class of nodes is distinguished by the property that all nodes have at least one short cut. Let *q*
_{p} be the probability that a typical node in this class is infected in a given time step. In practice, *q*
_{p}>*q*
_{K}, owing to the important role short cuts play in spreading the epidemic through the network. For example, simulations of the lattice model (equation (2.1)) under syphilis parameters show that *q*
_{p}≈4*q*
_{K}. For the model parametrized with the theoretical values taken from Kuperman & Abramson (2001), *q*
_{p}≈10*q*
_{K}. Incorporating this important observation in the difference equation leads to the following model of the *SIRS* dynamics on a complex contact network:
3.4
where the effective *q*, or *q*
_{eff}, is defined through the relation
3.5
and where *q*
_{K}<*q*
_{p}≪1. Note that equation (3.4) captures both the time delay dynamics resulting from the temporary immunity of the disease and the social effects of the short cuts. Moreover, one sees explicitly how the effective probability of infection *q*
_{eff} is composed of the distinct contributions—the nearest neighbours *q*
_{K} and the short cuts *q*
_{p}.

The dynamics of the model (3.4) depend on *p* in a manner that is very similar to the more complex contact network model (2.1). Figure 4 shows the results based on setting *q*
_{K}=0.11 and *q*
_{p}=0.65. For *p*=0, a small endemic equilibrium is reached (figure 4
*a*); for *p* in the small-world range, sustained oscillations arise (figure 4
*b*); and for large *p*, the disease is eradicated (figure 4
*c*) due to a synchronization effect. Comparing figures 2
*d* and 4
*b*, one can see exactly the same type of dynamics with a similar period of approximately 10 years as observed in the syphilis datasets (see CDC 1997–to date; Grassly *et al*. 2005).

### 3.1 Stability analysis

The equilibrium solutions of model (3.4) are found by solving the following equation for *I*
^{*}:
3.6
It is easily seen that the infection-free equilibrium *I*
^{*}=0 is always a solution. A stability analysis (based on the linearization of equation (3.4)) reveals that the infection-free equilibrium is stable when the following inequality holds:
3.7
Note that
3.8

An epidemic outbreak is only possible if the model parameters are such that the infection-free equilibrium is unstable (*R*
_{0}>1). It follows from equation (3.7) that, for the syphilis parameters in the small-world regime, this can occur only if the average number of contacts, *K*, is at least *K*=8. Thus, by choosing the parameter *K*≥8, at least some of the *p* values correspond to an unstable infection-free equilibrium, in the range covering the small-world regime. The instability region of the infection-free equilibrium grows with *p* and *K*, where for *K*=12, for example, the infection-free equilibrium is never stable for the syphilis parameters. This is visualized in figure 5
*a*, in which we plot the bifurcation curve of the infection-free equilibrium as a function of the parameters *K* and *p*. The infection-free equilibrium is stable for all values of parameters below the (lower solid) curve where *R*
_{0}<1 and unstable otherwise since *R*
_{0}>1. For example, it can be seen that, for *K*≤4, the infection-free equilibrium is stable in a wide strip of the parameter plane, containing the small-world regime.

A similar result is obtained by estimating the reproductive number *R*
_{0}, approximating it as the average number of secondary infectives produced by a typical infective individual in a sea of susceptibles. In the case presented here, where a single node may infect only those nodes it is linked to, the number of secondary infections may be approximated as
3.9
provided *q*
_{k} and *q*
_{p} are much less than 1. The above estimate gives a good approximation of the exact condition (3.7), as shown in figure 5
*a* (dashed curve).

For most parameter values, the difference equation (3.4) has a second endemic equilibrium in which *I*
^{*}>0. The bifurcation curves describing its stability in the (*p*, *q*
_{p}) parameter plane are shown in figure 5
*b*, for *K*=8 (dashed curve) and *K*=12 (solid curve). The curve indicates a Hopf bifurcation in the (*p*, *q*
_{p}) plane, where all other parameters are kept fixed. Appendix A provides technical details on the calculation of the curves. The endemic equilibrium exists for all parameter values for which the Hopf bifurcation curve exists. Below the bifurcation curve, the endemic equilibrium is stable and is unstable above it, where a stable limit cycle exists. The stable limit cycle thus coexists with the two unstable equilibria: infection free and endemic.

Hence, the bifurcation scenario, for the syphilis parameter values, for *K*=8, is as follows. For 0≤*p*≤0.0091, the infection-free equilibrium is stable and the disease goes extinct; for 0.0091<*p*≲0.079, the infection-free equilibrium loses its stability and an endemic equilibrium is born; for *p*≳0.079, the endemic equilibrium loses its stability through a Hopf bifurcation and a stable limit cycle is born (see the upper dashed Hopf bifurcation curve in figure 5
*b*). In the case of *K*=12, the infection-free equilibrium is unstable for all values of *p*. However, the endemic equilibrium is stable for 0≤*p*≲0.011. When *p*≳0.011, both the equilibria are unstable and, instead, a stable limit cycle is observed (see the lower solid Hopf bifurcation curve in figure 5
*b*).

As in the network model, the disease goes completely extinct for relatively large *p* values (for *K*=12, *p*
_{c}=0.43; for *K*=8, *p*
_{c}=0.87), despite the fact that the infection-free equilibrium is unstable (figure 4
*c*). The extinction should be attributed to the synchronization effect that takes place for these high *p* values. That is, the dynamics are such that a large proportion of the population becomes infected together and proceeds on to move to the recovered class together. The synchronization requires the initiation of a strong epidemic, implying that *R*
_{0} must be greater than unity, which explains why the infection-free equilibrium is unstable in this regime. Nevertheless, the disease becomes extinct due to the synchronization effect.

### 3.2 Vaccination

A *random vaccination* scheme may be incorporated into the difference equation model by replacing *S*
_{t} in equation (3.4) with (1−*v*)*S*
_{t}. The parameter *v* is the proportion of susceptibles vaccinated per time unit. Denote by *v*
_{e} the threshold proportion of vaccinated needed for disease extinction. A simple algebraic expression can be obtained for the extinction threshold by linearizing equation (3.4) about the infection-free equilibrium. Then, by using equation (3.7), disease extinction is reached if
3.10
For the parameter values of figures 4 and 5, the vaccination threshold is *v*
_{e}=0.530. Numerical simulations corroborate the existence of this threshold.

As the difference equation does not give any information regarding spatial patterns, it is impossible to apply the spatially oriented targeted vaccination scheme described for the CA model above. However, targeted vaccination schemes may nevertheless be explored by differentiating between vaccinating nearest neighbours and short cuts. This can be achieved by replacing the term *I*
_{t} in equations (3.1) and (3.4) with the term (1−*v*
_{K})*I*
_{t} for nearest neighbours and (1−*v*
_{p})*I*
_{t} for short cuts. The results reveal that the extinction threshold is very sensitive to, and is lowered dramatically by, the term *v*
_{p} for vaccinating short cuts, while the term for vaccinating nearest neighbours, *v*
_{K}, has little influence. Before applying vaccination, the model's infection-free equilibrium is never stable for the syphilis parameters (with *K*>8). However, targeted vaccination effectively reduces the number of the actual contacts of the key individuals, thereby reducing *R*
_{0}<1 in the vaccinated population.

## 4. Conclusion and discussion

Two models for studying the dynamics of diseases with temporary immunity in complex population networks have been proposed: a lattice model, which incorporates spatial information, and a difference equation model, which allows an analytic approach. The study focuses on the example of syphilis epidemics, which is a representative STD targeted to be eliminated in the USA, although with little success so far. The network model reveals that diseases with temporary immunity on a small-world contact network exhibit periodicity and waves of epidemics, some of which become pacemaker centres. It is shown that, by eliminating pacemakers through vaccination, the disease goes to extinction within one to two periods, where only approximately 20 per cent of the population requires vaccination. This is in contrast to standard vaccination programmes that set out to achieve herd immunity by vaccinating over 80 per cent of the population. The difference equation model allows further investigation of the Hopf bifurcation lying at the heart of the pacemaker phenomena. The two models complement each other, allowing a more profound view of the dynamics.

In reality, this control scheme may be implemented by vaccination, quarantine or a targeted education plan, depending on the disease and on the means available for its control. The main advantage of the vaccination methods proposed here is that they avoid the usual practice of vaccinating a large proportion of the population. In addition, vaccination is confined solely to relatively small and specified areas (figures 1
*c*,*d* and 2
*b*,*c*). In practice, it is always preferable to vaccinate as small a group as possible, as vaccination always carries a risk. Hence, it is advantageous to target only the relevant groups, already at risk. Note, however, that while other works refer to tracing infected individuals or the most connected individuals for applying a targeted vaccination scheme (e.g. Eames & Keeling 2002; Pastor-Satorras & Vespignani 2002; Zanette & Kuperman 2002), in the method proposed here no contact tracing is required, the pacemaker waves stem from the SIRS dynamics and the small-world structure in a natural and intrinsic way. Hence, these areas are easily identified as small areas where the infection appears repeatedly (figures 1 and 2).

In the case of syphilis, a vaccine is under development (e.g. Cullen & Cameron 2006), but there are already concerns regarding its safety. Thus, the model suggests that, if at all, vaccination of the carefully targeted centres of the population already at risk (e.g. in proximity to core groups of active individuals; see Thomas & Tucker 1996; Eames & Keeling 2002), rather than the full population, should still have the potential to result in disease extinction. Moreover, by treating only those individuals in high-risk pacemaker areas, it minimizes the application of the vaccination with its possible risks to the larger population. In the absence of a vaccine, education, early identification and the practice of safe sex are the available control measures (CDC 1997–to date; St Louis & Wasserheit 1998).

This work addresses the controversy as to whether syphilis epidemics recur approximately every 10 years because of the temporary immunity it endows to infected individuals (Grassly *et al*. 2005) and its associated time delay, or because of changing patterns in social behaviour (CDC 1997–to date; Douglas 2005). As shown here, both factors are crucial for recurrent syphilis epidemics. Thus, for example, oscillations cannot occur outside the small-world regime even in the presence of strong temporary immunity. For zero or very small *p* values (corresponding to none, or a very few short-cut links), epidemics cannot develop. Moreover, the analysis performed on the difference equation model reveals that, if *all* individuals in a small-world-type population network have only a few contacts, the infection-free equilibrium is stable. In addition, at the other end of the spectrum, it is pointed out that, for large enough *p* outside the small-world region (corresponding to many short-cut links), a synchronization effect occurs whereby epidemics are eradicated owing to exhaustion of susceptibles. By contrast, a society whose social behaviour approximates a small-world network with moderate heterogeneous levels of promiscuity would sustain the periodic recurrences of the syphilis epidemics approximately every 10 years.

The recent literature indicates another possible explanation for the syphilis oscillations. It has been proposed that similar stochastic epidemic oscillations may be an outcome of resonance-like effects (see McKane & Newman 2005; Risau-Gusman & Abramson 2007). This intriguing possibility requires further investigation. Finally, we conjecture that complete disease extinction is nevertheless achievable by a targeted vaccination scheme similar to that presented here. The targeting and vaccination of key individuals effectively reduces *R*
_{0} to less than unity in the vaccinated population, thereby leading to disease extinction.

Another issue to point out is that it seems most plausible that the nodes at a centre of the key pacemakers should be considered superspreaders, as discussed in Lloyd-Smith *et al*. (2005). By definition, a superspreader is ‘any case causing more infections than would occur in 99 per cent of infectious histories in a homogeneous population’. Owing to the structure of the models presented here, it is not possible to count how many infections are caused by any specific node per time step. This is because the exact source node of infection for any individual can never be identified; the formulation is based on probabilities. However, as discussed above, we know empirically that nodes with short cuts are infected at least four times more frequently than nodes without short cuts (for the theoretical model, nodes with short cuts are infected some 10 times more often). Hence, it is plausible that nodes with short cuts infect at least four times more than those without short cuts. In the small-world parameter regime, usually only approximately 10 per cent or less of the population has short cuts, so those are the only nodes that have at least four times larger probability of infecting someone than those who do not have short cuts (approx. 90% of the population or more). We have also shown numerically that it is sufficient to have a node with two short cuts to create a key pacemaker. So those individuals have at least eight times more probability of infecting a linked node, and they are less than 5 per cent of the population. Hence, it is most likely that the nodes at the centres of the key pacemakers are superspreaders, although a more accurate investigation of this plausibility is left for future work. It is interesting to note that, in this small-world framework, the whole structure is linked by short cuts within the whole population network, and thus short cuts might serve as a superspreading mechanism within the network. This is also left for further investigation.

## Acknowledgments

We thank Dana Torok for helping to develop the targeted vaccination scheme, and three reviewers for sharing their insights. We are grateful for the support of the James S. McDonnell foundation and the Israel Science Foundation.

## Appendix A. Technical details for the difference equation model

Here, we present the technical details of the stability analysis performed for the difference equation model (3.4). Consider the SIRS dynamics. Assume *τ*
_{I} is the amount of time units an individual spends in the infectious class, and that *τ*
_{R} represents the time units an individual later spends in the recovered class. As for syphilis *τ*
_{I}=1 (where a time unit represents six months), the proportion of recovered individuals can be described by the sum
. Thus, denoting by *S*
_{t} the proportion of susceptibles in the population at time *t*, we obtain
A1
where *τ*
_{0}=*τ*
_{I}+*τ*
_{R}. Set *I*
_{t} to be the proportion of infectives at time *t*, denote by *K* the number of connections an individual has on average and let *p* be the proportion of short-cut links an individual has among its *K* connections. Then, denote by
the number of infectious *nearest neighbours* an individual node has at time *t* and by
the number of infectious contacts via *short-cut links* an individual node has among its *K* connections at time *t*. Then,
A2
A3
Set *q*
_{p} as the probability of being infected via a short-cut link and set *q*
_{K} as the probability of being infected by a nearest neighbour, where *q*
_{K}<*q*
_{p}<1. Hence, the *SIRS* dynamics on a network with nearest-neighbour and short-cut connections can be described by
A4
The dynamics of (A 4) are at equilibrium for the solutions, *I*
^{*}, of the equation (A 5)
A5
It is easily seen that the infection-free equilibrium *I*
^{*}=0 is always a solution and that an endemic equilibrium *I*
^{*}>0 is a solution for most parameters relevant for syphilis.

Substituting *J*
_{t}=*I*
_{t}−*I*
^{*} and linearizing equation (A 4) about *I*
^{*}, we obtain
where

Now substitute *J*
_{t}=*J*
_{0}
*λ*
^{t} to obtain the characteristic polynomial
A6
where
A7
The stability of the infection-free equilibrium is derived from substituting *I*
^{*}=0 and requiring that |*λ*|=*R*
_{0}<1. The stability analysis of the infection-free equilibrium is presented in the main text (figure 5
*a*). Here, we provide details regarding the endemic equilibrium stability and the calculation of the Hopf bifurcation curve(s), presented in figure 5. This calculation is inspired by a calculation of the Hopf bifurcation of a mean-field difference equation model in Girvan *et al*. (2002).

Assume such a bifurcation exists and substitute *λ*=e^{iϕ} (as stability changes at |*λ*|=1) into equation (A 6). Separating the real and imaginary parts, we obtain two equations
A8
Substituting equations (A 7) into equations (A 8) and adding equation (A 5) results in a system of three equations in the seven variables: *q*
_{K}; *q*
_{p}; *p*; *K*; *I*
^{*}; *τ*
_{0}; and *ϕ*. Some of these parameters can be fixed to values relevant for syphilis: *τ*
_{0}=11; *K*=8 or 12; and *q*
_{K}=0.11. *ϕ* can be viewed as the frequency of the periodic solution emerging at the bifurcation, where the period of the limit cycle (when it exists) is *T*≈2*π*/*ϕ*. As for the syphilis parameter values, the period of oscillation is *T*≈20 time units, where 1 time unit=0.5 year=six months, we fix *ϕ*=0.3. Now, with these fixed parameter values, we use the Newton method to solve equations (A 8) and (A 5), using the syphilis parameter values as the initial guess, once for *K*=8 (dashed curve in figure 5
*b*) and once for *K*=12 (solid curve in figure 5
*b*). The results are plotted in the (*p*, *q*) plane (figure 5
*b*). In figure 5
*b*, the two bifurcation curves (one for *K*=8 (dashed) and one for *K*=12 (solid)) are presented, where below each curve an endemic equilibrium *I*
^{*}>0 is stable. At the bifurcation curve, the equilibrium loses its stability and a stable limit cycle is born so that, above the curve, a stable limit cycle coexists with two repelling (unstable) equilibria: the infection free and an endemic.

## Footnotes

- Received August 3, 2008.
- Accepted October 7, 2008.

- © 2008 The Royal Society