\section{Introduction}
\label{gcontlong.sec:intro}

Contagion, in its most general sense, is the spreading of an entity
or influence between individuals in a population, via direct or
indirect contact.  Contagion processes
therefore arise broadly in the social and
biological sciences, manifested as, for example, the spread of
infectious diseases~\citep{murray2002,daley1999,anderson1991a,brauer2001,diekmann2000,hethcote2000} and
computer viruses, the diffusion of innovations~\citep{coleman1966,valente1995,rogers2003}, 
political upheavals~\citep{lohmann1994}, and the dissemination of religious
doctrine~\citep{stark1996,montgomery1996}.  
Existing mathematical models of contagion, while motivated in a 
variety of ways depending on the application at hand,
fall into one of only two broad categories, where the critical distinction
between these categories can be explained in terms of the interdependencies between
successive contacts; that is, the extent to which the effect of an
exposure to a contagious agent is determined by the presence or absence
of previous exposures.

The standard assumption in all mathematical models of infectious
disease spreading (for example, the classic SIR model~\citep{kermack1927,murray2002}), 
and also in some models of social contagion~\citep{goffman1964,daley1965,bass1969}, 
is that there is no interdependency between contacts; rather, the infection probability is
assumed to be independent and identical across successive
contacts.  All such models therefore fall into a category that we call 
\textit{independent interaction models}.  
By contrast, what we call \textit{threshold models} assert that an individual can
only become infected when a certain critical number of
exposures has been exceeded, at which point infection becomes highly probable.
The presence of a threshold corresponds to interdependencies of an especially strong
nature: contacts that occur when an individual is near its threshold
are extremely consequential while others have little or no
effect. Threshold models are often used to describe social contagion
(e.g., the spreading of fads or rumors), where individuals either
deterministically~\citep{schelling1973,granovetter1978,watts2002} or
stochastically~\citep{bikhchandani1992,banerjee1992,morris2000,brock2001}
``decide'' whether or not to adopt a certain behavior based in part or in
whole on the previous decisions of others.

An alternative way to think about the interdependence of successive events
is in terms of memory: threshold models implicitly assume the presence
of memory while independent interaction models assume (again implicitly) that the
infection process is memoryless~\footnote{
  Note that independent interaction models typically incorporate a state of immunity,
  which reflects a type of memory, but not the kind with which we are chiefly concerned 
  here---namely, the memory of a contagious influence or entity \textit{prior} to an infection occurring.
}.
Neither class of model, however, is
able to capture the dynamics of contagious processes that possess an
intermediate level of interdependency, or equivalently a variable
emphasis on memory.  Furthermore, the relationship between
interdependent interaction models, threshold models, and any possible
intermediate models is unclear.  Motivated by these
observations, our model seeks to connect threshold and independent interaction models
both conceptually and analytically, and explore the classes of
contagion models that lie between them.  Such an
analysis is clearly relevant to problems of social contagion, in which
memory of past events obviously plays some role, but one that may be
less strong than is assumed by most threshold models.  However, it may
also be relevant to biological disease spreading models, which, to our knowledge,
have previously not questioned the assumption of independence between
successive exposures.  While the independence assumption is indeed
plausible, it has not been demonstrated empirically, and little enough
is understood of the dynamics of immune system response that the
alternative---persistent sub-critical doses of an infectious agent
combining to generate a critical dose---cannot be ruled out.
Furthermore, memory effects are known to be inherent to certain kinds of 
immune system responses, such as allergic response.  Hence if, as we indeed show, only a
slight departure from complete independence is required to alter the
corresponding collective dynamics, then our approach may also shed light
on the spread of infectious diseases.


\section{Descripton of Model}
\label{gcontlong.sec:model}

\begin{figure*}[tbp]
  \centering
  \includegraphics[width=\textwidth]{fig001_SIRtransitions6.pdf}
  \caption{
    (A) 
    Representation of how an individual's dose at time $t$ is determined.
    At each time step, each individual, \textit{regardless} of its state, contacts 
    one other randomly chosen individual.
    The probability that individual $i$ contacts an infective is then
    $\phi_t$, i.e., the current fraction of infectives in the
    population.  If the contact is infected
    then with probability $p$, $i$ is exposed to the contagious entity
    and receives a dose $d$ drawn from a fixed distribution $f$.
    Otherwise and also if the contacted individual is not infected
    in the first place (occurring with probability $1-\phi_t$),
    $i$ receives a dose of zero size.
    (B) 
    Individual $i$ then updates its dose count $D_{t,i}$ by `forgetting' the dose it received $T+1$ time
    steps ago and incorporating the current dose [\Req{gcontlong.eq:sumT}].
    (C)
    Transition probabilities for individual $i$ cycling through
    the three states S (susceptible), I (infected), and R (recovered).
    If $i$ is in the susceptible state, it becomes infected 
    with probability 1 once its dose count $D_{t,i}$
    exceeds its threshold $d_i^\ast$,
    otherwise it remains susceptible.  If $i$ becomes infected,
    then whenever $i$'s dose count
    drops below its threshold, it has a probability
    $r$ of recovering in each time step.  Once $i$ is in state R
    (where it is immune to infection), it becomes susceptible again 
    with probability $\rho$, again in each time step.
    Note that if $r=\rho=1$, infected individuals whose dose
    count falls below their threshold immediately return
    to the susceptible state.
  }
  \label{gcontlong.fig:SIRtransitions}
\end{figure*}

Our model, aspects of which we have reported in brief
elsewhere~\citep{dodds2004pa}, comprises a population of $N$
individuals, each of which is assumed to occupy
one of three states: \textit{susceptible} (S);
\textit{infected} (I); or \textit{removed} (R). 
Figure ~\ref{gcontlong.fig:SIRtransitions} 
provides a schematic representation of our model which
we describe as follows.
At each (discrete) time step $t$,
each individual $i$ comes into contact with another individual chosen
uniformly at random from the population.
If the contact is infected---an event that occurs with probability $\phi_t$, the fraction
of infectives in the population---then with probability $p$, $i$ receives
a `dose' $d$ drawn from a fixed dose-size distribution $f$; else
$i$ receives a dose of size $0$. We call a successful transmission of a
positive dose an \textit{exposure} and $p$ the \textit{exposure probability}. 
Individuals carry a
\textit{memory} of doses received from their last $T$ contacts and we
denote the sum of individual $i$'s last $T$ doses ($i$'s \textit{dose count}) 
at the $t$th time step by
\begin{equation}
  \label{gcontlong.eq:sumT}
  D_{t,i} = \sum_{t'=t-T+1}^{t} d_{t',i}.
\end{equation}
If $i$ is in the susceptible state, then it becomes infected once $D_{t,i}$ exceeds
$i$'s \textit{dose threshold} $d_i^\ast$, where $d_i^\ast$ is drawn from a given
distribution $g$ (dose thresholds do not change with time). Note
that we differentiate exposure from infection, the latter being the
possible result of one or more exposures and only occurring once a
susceptible individual's threshold has been equaled or
exceeded.  Having become infected, an individual remains in state I until its dose
count drops below its threshold, at which point it recovers with
probability $r$ at each time step.  Once recovered, an individual
returns to being susceptible with probability $\rho$, again at
each time step.  

The probability $P_{\rm inf}$ that a susceptible individual 
who comes into contact with $K$ infected individuals in $T$
time steps will become infected
is therefore given by
\begin{equation}
  \label{gcontlong.eq:Pinfbare}
  P_{\rm inf}(K)
  =
  \sum_{k=1}^{K}
    \binom{K}{k}
    p^k
    (1-p)^{K-k}
    P_k,
\end{equation}
where $K=1,\ldots,T$, and 
\begin{equation}
  \label{gcontlong.eq:phifixgendosedstar0}
  P_k
  = 
  \int_{0}^{\infty} \dee{d^\ast}
  g(d^\ast)
  \int_{d=d^\ast}^{\infty} \dee{d}\, f^{k\star}(d)
\end{equation}
Both $P_{\rm inf}(K)$ and $P_k$ are important quantities in our model.
The quantity $P_k$ is the expected fraction of a population that will be infected
by $k$ exposures, where the distribution $f^{k\star}(d)$, the $k$-fold
convolution of $f$, is the probability that the sum of $k$ doses
will be equal to $d$.  
The infection probability $P_{\rm inf}(K)$ gives, in effect, a
\textit{dose response curve}~\citep{haas2002} averaged over all members
of the population and also the distribution of dose sizes (where we
note that $K$ contacts with infected individuals will result in $k$
actual exposures with probability $\binom{K}{k} p^k (1-p)^{K-k}$).
Figure~\ref{gcontlong.fig:doseresponse30} shows examples of dose
response curves, calculated from \Req{gcontlong.eq:Pinfbare}
for four configurations of the model. The plots correspond to (A)
independent interaction, (B) deterministic threshold, and, in both (C) and (D),
stochastic threshold models.  For the independent interaction example, a single
exposure is needed to generate an infection and so exposures
effectively act independently.  The deterministic threshold example
incorporates uniform dose sizes and thresholds, and 
when the probability of an exposure is set to $p=1$, the
response becomes deterministic (individuals are always infected when
their dose count is met or exceeded and never otherwise); but now
the threshold can only be exceeded by multiple infections.
The two
stochastic cases generalize the deterministic case by allowing (C)
dose sizes to be heterogeneous, and (D) both dose sizes and thresholds
to be heterogeneous.

\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.49\textwidth]{fig002_figdoseresponse31_noname.pdf}
  \caption{
    Examples of dose response, i.e., 
    $P_{\rm inf}(K)$, the probability an individual becomes infected due to
    $K$ exposures in the last $T (=12)$ time steps 
    [see \Req{gcontlong.eq:Pinfbare}].
    The plots correspond to 
    (A) independent interaction (or disease-like) models with homogeneous
    dose sizes and thresholds 
    [$p=0.25$, dose distribution $f(d) = \delta(d-1)$,
    and threshold distribution $g(d^\ast)=\delta(d^\ast-1)$];
    (B) deterministic threshold models with homogeneous
    dose sizes and thresholds 
    [$p=1$, $f(d) = \delta(d-1)$, and $g(d^\ast)=\delta(d^\ast-5)$];
    (C) stochastic threshold models with
    heterogeneous doses and homogeneous thresholds
    [$p=0.9$, doses distributed lognormally with unit
    mean and standard deviation $0.5$, and $g(d^\ast)=\delta(d^\ast-5)$];
    and
    (D), same as (C) but with heterogeneity in thresholds introduced via
    a lognormal distribution with mean of 5 as per (C) and standard deviation of 10.
  }
  \label{gcontlong.fig:doseresponse30}
\end{figure}

We explore the behavior of our model with respect to
three qualitative types of dynamics: (1) permanent removal ($\rho=0$)
dynamics, analogous to so-called SIR models in mathematical
epidemiology in which individuals either die or acquire permanent
immunity; (2) temporary removal ($1>\rho>0$) dynamics, analogous to
SIRS models where recovered individuals become susceptible again after
a certain period of immunity; and (3) instantaneous replacement ($\rho=1$)
dynamics, analogous to SIS models, where infected individuals
immediately become susceptible again upon recovery. Chicken pox, for
example, would correspond to SIR-type contagion, while the common cold
resembles the SIRS case.  Because of its simplicity,
we obtain the majority of our analytic results for the somewhat special SIS case
($\rho=1$).  However, our main findings for the SIS case have analogs in
the more complicated SIRS and SIR cases which we investigate with numerical
simulations.  Furthermore, while the
assumption of instantaneous re-susceptibility is probably not
appropriate in the case of contagious biological agents (where
recovery is generally associated with some period of immunity), it may
well be approximately true for contagious social influences, such as
``social smoking,'' where an individual, having quit, can restart
immediately. A summary of the main parameters of the model and their
definitions is provided in Table~\ref{gcontlong.tab:defns}.

\begin{table}[htbp]
  \begin{center}
    \begin{tabular}{c|l}
      \hline
      $T$ & length of memory window \\
      $p$ & probability of exposure given contact with infective\\
      $r$ & probability of moving from infected to recovered state\\
      $\rho$ & probability of moving from immune to susceptible state \\
      $f(d)$ & distribution of dose sizes $d$ \\
      $g(d^\ast)$ & distribution of individual thresholds $d^\ast$\\
      $\bar{d^\ast}$ & uniform threshold of homogeneous population \\
      $\phi_t$ & fraction of population infected at time $t$ \\
      ${\phi^{\ast}}$ & steady-state fraction of population infected \\
      \hline
    \end{tabular}
    \caption{
      Summary of main model parameters and definitions.
      }
    \label{gcontlong.tab:defns}
  \end{center}
\end{table}

We structure the remainder of the paper as follows. In
section~\ref{gcontlong.sec:homog}, through analysis and simulations,
we examine in detail the SIS $(\rho=1)$ version of the model for a
homogeneous population, where by ``homogeneous,'' we mean that all doses
are equal and of unit size [i.e., $f(d)=\delta(d-1)$], and that all individuals
have the same threshold [i.e., $g(d^\ast) =
\delta(d^\ast-\bar{d^\ast})$]. For homogeneous populations, we find
that only two universal classes of dynamics are possible: (1)
\textit{epidemic threshold dynamics}~\footnote{
  In this paper, we use the word `threshold' in two terms:
  `threshold models' and `epidemic threshold models.'
  At the risk of some confusion, we have done so to maintain consistency
  with the two distinct literatures we are connecting, sociology and mathematical epidemiology.
  Threshold models of sociology refer to individual level thresholds whereas
  the term epidemic threshold refers to the critical reproduction number of a disease,
  above which an outbreak is assured.
},
according to which initial
outbreaks either die out or else infect a finite fraction of the
population, depending on whether or not the infectiousness $p$ exceeds
a specific critical value $p_c$; 
and (2) \textit{critical mass dynamics}
according to which a finite fraction of the population can only ever
be infected in equilibrium if the initial outbreak size itself
constitutes a finite ``critical mass.''  Although homogeneity is
a restrictive assumption for biological or social
contagion, it provides a useful special case in that it illuminates
the basic intuitions required to understand the more general,
heterogeneous case.  In section~\ref{gcontlong.sec:heterog}, we relax
the homogeneity assumption and move to the richer and more realistic
case of arbitrarily distributed dose sizes and individual thresholds.
Here we find that three universal classes of contagion models are
possible.  In addition to the epidemic threshold and critical
mass classes that carry over from the homogeneous case, we also find an
intermediate class of \textit{vanishing critical mass dynamics} in
which the size of the required critical mass diminishes to zero for
$p<1$. Furthermore, we determine where the transitions between these
classes occur and also the conditions required for more complicated
kinds of contagion models to arise. Subsequently, in
section~\ref{gcontlong.sec:SIR}, we explore the SIRS and SIR versions
of the model, finding behavior that in many ways resembles that of
the simpler SIS case.  In the SIR case, for example, it is necessarily
true that all epidemics eventually burn themselves out (because for
$\rho=0$, the removed condition is an absorbing state).  However, we
find that the presence of memory may cause an epidemic to persist for a
surprisingly long time.
In section~\ref{gcontlong.sec:concl}, we conclude our
analysis, discussing briefly the implications of our findings for
stimulating or retarding different contagious processes.  Finally, in
the appendices we provide detailed derivations of the analytical
results from sections~\ref{gcontlong.sec:homog}
and~\ref{gcontlong.sec:heterog}.

\section{Homogeneous SIS contagion models}
\label{gcontlong.sec:homog}

\subsection{Epidemic threshold models}
\label{gcontlong.subsec:epithreshmodels}

We begin our analysis with a simple non-trivial case
of our model, for which we assume that individuals are:
identical (i.e., the population is homogeneous);
have no memory of past interactions $(T=1)$; 
and, upon recovery from infection, immediately
return to the susceptible state (i.e., $\rho=1$). 
In this limit our model coincides with the SIS version of the
traditional Kermack-McKendrick model~\citep{kermack1927,murray2002}, as
individuals with no memory necessarily become infected upon exposure
to a single infected individual [$g(d^\ast)=\delta(d^\ast-1)$].  For a
specified recovery time $(r\leq 1)$ (which again is identical across
all individuals), the fraction of infected individuals at time $t$,
$\phi_{t}$, evolves according to
\begin{equation}
  \label{gcontlong.eq:r<1k1T1pre}
  \phi_{t+1} = p \phi_{t}
  + \phi_t(1 - p\phi_t)(1-r).
\end{equation}
The first term on the right is the fraction of individuals
newly infected between time $t$ and $t+1$, regardless of
whether or not they were infected beforehand (the model
allows for individuals to recover and be reinfected
within one time step).  The second term 
is the fraction of individuals who were infected in the 
preceding time step, were not infected between time $t$ and $t+1$, and did not recover.

Setting $\phi_t={\phi^{\ast}}$ in \Req{gcontlong.eq:r<1k1T1pre},
we find the stable fixed points of the model as a function of $p$ are given by
\begin{equation}
  \label{gcontlong.eq:r<1k1T1}
{\phi^{\ast}} = 
\left\{
  \begin{array}{ll}
    0 & 
    \mbox{for} \quad 0 \le p \le r, \\
    \frac{1 - r/p}{1-r} & 
    \mbox{for} \quad r < p \le 1.
  \end{array}
\right.
\end{equation}
Furthermore, a single unstable set of fixed points is found along
${\phi^{\ast}}=0$ for $r < p \le 1$.  Thus, the standard memoryless SIS
model exhibits an \textit{epidemic threshold}~\citep{murray2002} at
$p=p_c$, as displayed graphically in
Figure~\ref{gcontlong.fig:gc_rvar_k1_Tvar}: for $p \le p_c = r$ no
infection survives (${\phi^{\ast}} \equiv 0$) while for $p > p_c = r$, a
stable, finite fraction of the population will become
infected (${\phi^{\ast}} > 0$). In the language of dynamical systems theory,
the epidemic threshold is a transcritical bifurcation~\citep{strogatz1994}
which lies at the intersection of two fixed point curves,
where the stable curve changes to unstable at
the intersection and vice versa (in statistical mechanics,
such behavior is referred to as a second-order or continuous
\textit{phase transition}, where $(p_c,0)$ is called the
\textit{critical point}~\citep{goldenfeld92}).
For the choice of parameters above, one
branch of the transcritical bifurcation is the $p$-axis which is stable
to the left of $p_c$ and unstable to the right (imagining the
unphysical extension of ${\phi^{\ast}}$ to ${\phi^{\ast}}<0$, the other branch may
be seen to extend from below the $p$-axis where it is unstable to
above the $p$-axis where it is stable).

\begin{figure}[tbp!]
  \centering
  \includegraphics[width=0.49\textwidth]{fig003_figgc_r0p50_k1_T2_paper2_002_noname.pdf}                                        
  \caption{
    Fixed point curves for an example epidemic threshold model.
    The stable fixed point curves shown are from 
    simulation (circles) and \Req{gcontlong.eq:r<1k1T>1}
    (blue line) for $\bar{d^\ast}=1$, $T=2$, and $\tau=1/r-1=1$.
    Unstable fixed points are indicated by the dashed red line.
    The trajectories of three initial conditions are shown
    to illustrate how the level of an epidemic evolves for different values of $p$.
    The epidemic threshold $p_c=1/3$ is as predicted by 
    $p_c = (T + \tau)^{-1}$, \Req{gcontlong.eq:r<1k1T>1pc}.
    For the simulation, the population size $N=10^5$, the
    number of time steps $N_t=10^3$, each data point
    represents the average value of ${\phi^{\ast}}(p)$ over the last $N_s=100$
    time steps, and the initial condition is that all are infected.
  }
  \label{gcontlong.fig:gc_rvar_k1_Tvar}
\end{figure}

We observe this kind of bifurcation structure (i.e., where the rising branch has
positive slope and comprises stable equilibria) to be 
a robust feature with respect to a range of parameter choices, and
we classify all such models as \textit{epidemic threshold models}.
As we show below, the same qualitative equilibrium behavior is present in 
all homogeneous models with dose thresholds set at $d^\ast=1$,
even with arbitrary memory $T \ge 1$ and recovery rate $r \le 1$. 
While some details of the dynamics do change as $T$ and $r$ are varied,
the existence of a single transcritical bifurcation depends only on the
assumption that all individuals exhibit ``trivial
thresholds:'' $g(d^\ast)=\delta(d^\ast-1)$.  

First, allowing memory to be arbitrary ($T>1$) but keeping $r=1$,
we observe that the individuals who are infected at some time $t$ are
necessarily those who have experienced at least one infectious event
in the preceding $T$ time steps; thus we obtain the following implicit
equation for ${\phi^{\ast}}$,
\begin{equation}
  \label{gcontlong.eq:r1k1T>1}
  {\phi^{\ast}} = 1 - (1 - p{\phi^{\ast}})^T.
\end{equation}
As with the $T=1$ case above, the equilibrium behavior exhibits a
continuous phase transition and therefore also falls into the epidemic
threshold class. While we can no longer find a general, closed-form
solution for ${\phi^{\ast}}$, \Req{gcontlong.eq:r1k1T>1} can be rearranged
to give $p$ as an explicit function of ${\phi^{\ast}}$:
\begin{equation}
  \label{gcontlong.eq:r1k1T>1_inv}
  p = {\phi^{\ast}}^{-1} [ 1 - (1-{\phi^{\ast}})^{1/T} ].
\end{equation}
Taking the limit of ${\phi^{\ast}} \rightarrow 0$ in either
\Req{gcontlong.eq:r1k1T>1} or \Req{gcontlong.eq:r1k1T>1_inv}, we find
$p_c=1/T$. 
Thus, receiving at
least one exposure from the last $T$ contacts is analogous to the
zero memory ($T=1$), variable $r$ case above where recovery occurs on a time scale $1/r
\simeq T$.

Next, we also allow the recovery rate to be arbitrary $(r<1)$. 
To find the fixed point curves, we modify
\Req{gcontlong.eq:r1k1T>1} to account for the fraction of individuals
that have not experienced a single exposure for at least the last $T$
time steps but have not recovered from a previous infection.  We first
write down the probability that an individual last experienced an
infectious event $m$ time steps ago and has not yet recovered.
Denoting the sequence of a positive unit dose followed by $m$ 0's as
$H_{m+1}$, we have
\begin{equation}
  \label{gcontlong.eq:k=1_b}
  P(\mbox{infected}|H_{m+1}) = p{\phi^{\ast}} (1 -p{\phi^{\ast}})^{m}
  (1-r)^{m-T+1},
\end{equation}
The first term on the right hand side of the \Req{gcontlong.eq:k=1_b}
is the probability of a successful exposure; the second term is the
probability of experiencing no successful exposures in the subsequent
$m$ time steps; and the final term is the probability that, once the
memory of the initial single exposure has been lost after $T$ time
steps, the individual remains infected.  Since we are only concerned
with individuals who have `forgotten' the source of the infection, we
have $m \ge T$.  Summing over $m$, we obtain the total probability
that an individual was infected at least $T$ time steps ago and has
not yet recovered:
\begin{equation}
  \label{gcontlong.eq:k=1inf}
  \sum_{m=T}^{\infty} P(\mbox{infected}|H_{m+1}) = \frac{p{\phi^{\ast}}
  (1-p{\phi^{\ast}})^T (1-r)}{1 - (1-p{\phi^{\ast}})(1-r)}
\end{equation}
Adding this fraction to the right hand side of
\Req{gcontlong.eq:r1k1T>1} then gives
\begin{eqnarray}
  \label{gcontlong.eq:r<1k1T>1}
  {\phi^{\ast}} & = & 1 - (1-p{\phi^{\ast}})^T \left[
    1 -
    \frac{p{\phi^{\ast}}(1-r)}
    {1 - (1-p{\phi^{\ast}})(1-r)}
  \right], \nonumber \\
  & = & 1 - \frac{r (1-p{\phi^{\ast}})^T }
    {1 - (1-p{\phi^{\ast}})(1-r)}.
\end{eqnarray}
Figure~\ref{gcontlong.fig:gc_rvar_k1_Tvar} shows a comparison
between the above equation and simulation results
for $T=2$, $r=1/2$, and $\bar{d^\ast}=1$.

Taking the limit of small ${\phi^{\ast}}$ in \Req{gcontlong.eq:r<1k1T>1}
we find the epidemic threshold to be
\begin{equation}
  \label{gcontlong.eq:r<1k1T>1pc}
  p_c = \frac{r}
    {1 + r(T-1)} = \frac{1}{T + 1/r - 1}.
\end{equation}
Checking the special cases of the preceding calculations,
we find $p_c=r$ when $T=1$ and $p_c=1/T$ when $r=1$.
Denoting the mean time to recovery
of an infected, isolated individual by $\tau$,
we observe $\tau=1/r-1$ and \Req{gcontlong.eq:r<1k1T>1pc} becomes
\begin{equation}
  \label{gcontlong.eq:r<1k1T>1pc2}
  p_c = \frac{1}{T + \tau}.
\end{equation}
The time scales $T$ and $\tau$ can thus both be thought of as 
corresponding to two different kinds of memory, 
the sum of which---the total characteristic time scale of memory in the
model---determines the position $p_c$ of the epidemic threshold.
Qualitatively, therefore, all homogeneous models that possess trivial individual
thresholds exhibit the same kind of equilibrium dynamics.  Varying the
thresholds, however, produces equilibrium behavior of a quite distinct
nature, as we show in the next section.

\subsection{Critical mass models}
\label{gcontlong.subsec:critmassmodels}

\begin{figure}[tbp!]
  \centering
  \includegraphics[width=0.49\textwidth]{fig004_figgc_r1_k3_T12_paper2_002_noname.pdf}
  \caption{
    Comparison between simulation (circles) and theory (lines, \Req{gcontlong.eq:r=1T>1k>1}) for
    $T=12$, $\bar{d^\ast}=3$, and $r=1$. 
    For $\bar{d^\ast}>1$, the homogeneous SIS contagion model
    exhibits a saddle-node bifurcation.
    Shown are the non-zero stable and unstable points.
    For the theoretical curves, solid lines represent stable points
    and dashed ones unstable points.
    All points on the line ${\phi^{\ast}}=0$ are stable points.
    The arrows show trajectories of the system for
    three example initial conditions.
    Simulations details are as per Figure~\ref{gcontlong.fig:gc_rvar_k1_Tvar}.
    The location of the unstable fixed point curve is determined by binary search.
  }
  \label{gcontlong.fig:gc_r1_k3_T12}
\end{figure}

When individuals of a homogeneous population require more than one
exposure to become infected---that is, 
when $g(d^\ast) = \delta(d^\ast-\bar{d^\ast})$ and $\bar{d^\ast}>1$---the closed-form
expression for the fixed points ${\phi^{\ast}}$ in the case of $r=1$ is
\begin{equation}
  \label{gcontlong.eq:r=1T>1k>1}
  {\phi^{\ast}} = 
  \sum_{i=\bar{d^\ast}}^{T}
  \binom{T}{i}
  (p{\phi^{\ast}})^{i} (1 - p{\phi^{\ast}})^{T-i}.
\end{equation}
We now observe a fundamentally different behavior of the model: the
transcritical bifurcation that is characteristic of epidemic threshold
models is absent and is replaced by a saddle-node
bifurcation~\citep{strogatz1994}, an example of which is illustrated in
Figure~\ref{gcontlong.fig:gc_r1_k3_T12}.  We call models whose
equilibrium states are determined in this manner \textit{critical mass
models} because, as indicated by the arrows in
Figure~\ref{gcontlong.fig:gc_r1_k3_T12}, an epidemic will not spread
from an infinitesimal initial outbreak, requiring instead a finite
fraction of the population, or ``critical mass''~\citep{schelling1978},
to be infected initially.  Saddle-node bifurcations (or backwards bifurcations) 
have also been observed in a number of unrelated epidemiological 
multi-group models, arising from differences between groups and inter-group 
contact rates~\cite{hadeler1995,hadeler1997,dushoff1998,kribs-zaleta2000,greenhalgh2000,kribs-zaleta2002}.

Although it is not generally possible to obtain an expression for ${\phi^{\ast}}(p)$
from \Req{gcontlong.eq:r=1T>1k>1}, we are able to write down
a closed-form expression involving the position of the 
saddle-node bifurcation $(p_b,{\phi_b^{\ast}})$ 
(see Appendix~\ref{gcontlong.app:saddlenodeextra} for details):
\begin{equation}
  \label{gcontlong.eq:r=1T>1k>1diff2}
  0 = 
  \sum_{i=\bar{d^\ast}}^{T} 
  \binom{T}{i}
  z^{i-\bar{d^\ast}} (1 - z)^{T-i-1}
  [i -1 - z(T-1)],
\end{equation}
where $z = p_b {\phi_b^{\ast}}$.
Using~\Req{gcontlong.eq:r=1T>1k>1diff2}, we solve for
$z$ by standard numerical means and then use~\Req{gcontlong.eq:r=1T>1k>1}
to obtain ${\phi_b^{\ast}}$ and hence $p_b$.
Figure~\ref{gcontlong.fig:gc_bipts_r1} shows positions of saddle-node
bifurcation points computed for a range of values of $T$ and
$\bar{d^\ast}$. In all cases, as $\bar{d^\ast}$ increases, the
bifurcation point moves upward and to the right of the fixed point
diagram.  However, for small values of $\bar{d^\ast}$ and $T$,
we are able to determine $(p_b,{\phi_b^{\ast}})$ exactly.
For example, for $\bar{d^\ast}=2$ and $T=3$, we find that $p_b=8/9$ and
${\phi_b^{\ast}}=27/32$ and that the bifurcation is parabolic.  

\begin{figure}[tbp!]
  \centering
    \includegraphics[width=0.49\textwidth]{fig005_figgc_bipts_r1_002_noname.pdf}
    \caption{
      Bifurcation points for $\bar{d^\ast}>1$, ${\phi^{\ast}}>1$, and $r=1$,
      determined numerically using equations~\Req{gcontlong.eq:r=1T>1k>1}
      and~\Req{gcontlong.eq:r=1T>1k>1diff2}.  The data shown are for
      $T=3$ ($\bigcirc$), 
      $T=6$ ($\Box$),
      $T=12$ ($\triangleleft$),
      $T=24$ ($\triangleright$),
      and
      $T=96$ ($\vartriangle$).
      The solid lines guide the eye for the range $\bar{d^\ast}=2$ to $\bar{d^\ast}=T-1$
      with $\bar{d^\ast}$ increasing with $p_b$.
      The dashed lines connect to the transcritical bifurcation points 
      ($\diamond$) observed for $\bar{d^\ast}=1$.  
      Note that no bifurcations occur when $\bar{d^\ast}=T>1$ (the sole fixed point is ${\phi^{\ast}}=1$ when $p=1$)
      and when $\bar{d^\ast}=T=1$ (the fixed points lie along the line $p=1$ and $0 \le \phi \le 1$).
      }
    \label{gcontlong.fig:gc_bipts_r1}
\end{figure}

As with epidemic threshold models, the analysis for critical mass
dynamics can be generalized to include a non-trivial recovery rate
$r<1$.  Now, however, we are unable to find a closed-form
expression for the fixed points for general $T$ and
$\bar{d^\ast}$.
The difficulty in making such a computation lies in finding
an expression for the number of individuals whose dose counts
are below threshold but are still infected since they have not yet recovered.
For the $\bar{d^\ast}=1$ case, this was straightforward since 
the only way to stay below threshold
was to experience a sequence of null exposures. 
Nevertheless, for $r<1$, we can formally modify the expression for ${\phi^{\ast}}$ 
given in \Req{gcontlong.eq:r=1T>1k>1}:
\begin{equation}
  \label{gcontlong.eq:Gamma2}
  {\phi^{\ast}} = \Gamma(p,{\phi^{\ast}};r,T) + \sum_{i=\bar{d^\ast}}^{T}
  \binom{T}{i} (p{\phi^{\ast}})^{i} (1 - p{\phi^{\ast}})^{T-i},
\end{equation}
where the additional term $\Gamma(p,{\phi^{\ast}};r)$ accounts for the
proportion of below threshold individuals who have not yet recovered.
For small values of $T$ and $\bar{d^\ast}$, exact expressions for
$\Gamma(p,{\phi^{\ast}};r)$ can be derived and then ${\phi^{\ast}}$ can be solved
for numerically. In Appendix~\ref{gcontlong.app:exactdstar=2}, we
consider two special cases, $\bar{d^\ast}=2$ for $T=2$ and
$T=3$, that illustrate the process of constructing expressions for
$\Gamma(p,{\phi^{\ast}};r)$. 
These results allow us to explore the movement of the fixed points with decreasing $r$.
Figures~\ref{gcontlong.fig:gc_T3_k2_bif_theorycomp2} A and B show
comparisons over a range of $r$ between the solutions of
\Req{gcontlong.eq:Gamma2} and simulations for $T=2$ and
$T=3$ respectively, confirming that the agreement is excellent.

\begin{figure}[tbp!]
  \centering
  \includegraphics[width=0.49\textwidth]{fig006_figgc_T3T2_k2_bif002_noname.pdf}
  \caption{
    A: Theory versus simulation for $T=3$, $\bar{d^\ast}=2$, and varying $r$.  
    Solid lines represent theoretically derived stable fixed points,
    dashed lines represent unstable ones, and the squares indicate bifurcation points
    (using Eqs.~\req{gcontlong.eq:Gamma2} and~\req{gcontlong.eq:Gammafull}).
    From left to right, we have $r=0.01, 0.05, 0.10, 0.15, \ldots, 1.00$.
    The circles correspond to simulation data for the case $r=0.2$
    and matches the theoretical curves.  
    B: Theoretical fixed point curves for $T=2$, $\bar{d^\ast}=2$, and varying $r$.  
    Bifurcations appear for $r < 0.3820 \pm 0.0001$.  Otherwise, $(p,{\phi^{\ast}})=(1,1)$ is
    the only fixed point of the system with ${\phi^{\ast}}>0$.  The values of $r$ are as per the $T=3$ case in A.
    Circles correspond to simulation results for $r=0.1$.
    For all values of $r$, both of the above systems possess a line of stable 
    fixed points described by $0 \le p < 1$ and ${\phi^{\ast}}=0$, i.e., the $p$-axis.
  }
  \label{gcontlong.fig:gc_T3_k2_bif_theorycomp2}
\end{figure}

\section{Heterogeneous SIS contagion models}
\label{gcontlong.sec:heterog}

\subsection{Distributions of doses and thresholds}
\label{gcontlong.subsec:dosethreshdist}

In real populations, both in the context of
infectious diseases and also social influences, individuals
evidently exhibit varying levels of susceptibility.  
Furthermore, contacts between infected and susceptible
individuals can result in effective exposures of variable size,
depending on the individuals in question, as well as the nature of 
their relationship and the circumstances of the contact (duration,
proximity, etc.).  Our homogeneity assumption of the preceding section
is therefore unlikely to be justified in any real application.  
Furthermore, as we show below, 
the more general case of heterogeneous populations, while preserving much of the
structure of the homogeneous case, also yields a new class of dynamical
behavior; thus the inclusion of heterogeneity not only makes the model
more realistic, but also provides additional qualitative insight.

While in principle all parameters in the model could be assumed to
vary across the population, we focus here on two
parameters---individual threshold $d^\ast$ and dose size $d$---which
when generalized to stochastic variables, embody the variations in
individual susceptibility and contacts that we wish to
capture. We implement these two sources of heterogeneity
as follows.  In the case of thresholds, each individual
is assigned a threshold drawn randomly from a specified probability
distribution $g(d^\ast)$ at $t=0$ which remains fixed for all $t$.
In effect, this assumption implies that individual characteristics remain roughly invariant on the time
scale of the dynamics, rather than varying from moment
to moment.  By contrast, in order to capture the unpredictability of circumstance,
we assign dose sizes stochastically,
according to the distribution $f(d)$, independent both of time
and the particular individuals between whom the contact occurs.

Again considering first the more tractable $r=1$ version of the SIS
model, we first examine the effect of allowing dose size $d$ to vary while
holding thresholds $d^\ast$ fixed across the population.  In the
homogeneous case, $k$ exposures of a susceptible to infected
individuals resulted in a dose count of $k$, but the result can
be more complicated when $d$ is allowed to vary continuously. Note that $d^\ast$ also no
longer need be an integer. First, we calculate the probability that a
threshold will be exceeded by $k$ doses. As the distribution of dose
size is now some arbitrary function $f$, we have that the probability
distribution of the sum of $k$ doses is given by the $k$-fold
convolution
\begin{equation}
  \label{gcontlong.eq:kfoldconv}
  \mathbf{P}\left(
    \textstyle{\sum_{j=1}^{k} d_j = d}
    \right)
    =
    f \ast f \ast \cdots \ast f (d) = f^{k\star}(d).
\end{equation}
The probability of exceeding $d^\ast$ is then
\begin{equation}
  \label{gcontlong.eq:kfoldconvexceed}
  \mathbf{P}\left(
    \textstyle{\sum_{j=1}^{k} d_j \ge d^\ast}
    \right)
    =
    \int_{d=d^\ast}^{\infty} \dee{d} \, f^{k\star}(d).
\end{equation}
Since, for $r=1$, individuals are only infected when their dose
count exceeds their threshold, we have that the steady-state 
fraction infected is given by
\begin{equation}
  \label{gcontlong.eq:phifixgendose}
  {\phi^{\ast}} =
  \sum_{k=1}^{T}
  \binom{T}{k}
  (p{\phi^{\ast}})^{k}
  (1-p{\phi^{\ast}})^{T-k}
    \int_{d=d^\ast}^{\infty} \dee{d} \, f^{k\star},
\end{equation}
where we have averaged over all possible ways
an individual may experience $1 \le k \le T$ exposures in $T$
interactions.

Next, in order to account for any variation $d^\ast$,
we must incorporate another layer of 
averaging into~\Req{gcontlong.eq:phifixgendose} as follows:
\begin{align}
  \label{gcontlong.eq:phifixgendosedstar}
  & {\phi^{\ast}} =  \nonumber  \\
  & \int_{0}^{\infty} \! \! \! \dee{d^\ast}
  g(d^\ast) 
  \sum_{k=1}^{T}
  \binom{T}{k}
  (p{\phi^{\ast}})^{k}
  (1-p{\phi^{\ast}})^{T-k} 
 \int_{d=d^\ast}^{\infty} \! \! \! \dee{d}\, f^{k\star}(d), \nonumber  \\
  & = 
  \sum_{k=1}^{T}
  \binom{T}{k}
  (p{\phi^{\ast}})^{k}
  (1-p{\phi^{\ast}})^{T-k} 
  P_k,
\end{align}
where $P_k$ is defined by \Req{gcontlong.eq:phifixgendosedstar0}.
An important insight that can be derived immediately from 
\Req{gcontlong.eq:phifixgendosedstar} is that all information
concerning the distributions of dose sizes and thresholds is
expressed via the $\{P_k\}$; hence the details of
the functions $f$ and $g$ are largely unimportant.  
In other words, many pairs of $f$ and $g$
can be constructed to give rise to the same $\{P_k\}$ and hence the same
fixed points.  For example, any desired $\{P_k\}$ can be 
realized by a uniform distribution of unit doses $f(d) = \delta(d-1)$,
along with a discrete distribution of thresholds
\begin{equation}
  \label{gcontlong.eq:univPk}
  g(d^\ast) = P_1 \delta(d^\ast-1) + \sum_{k=1}^{T-1}(P_{k+1}-P_{k}) \delta(d^\ast-k-1).
\end{equation}

We observe that for the homogeneous case where doses and thresholds
are fixed at $1$ and $\bar{d^\ast}$ respectively [i.e., $f(d) =
\delta(d-1)$ and $g(d^\ast) = \delta(d^\ast-\bar{d^\ast})$], we have
$f^{k\star}(d) = \delta(d-k)$ and the expression for $P_k$
[\Req{gcontlong.eq:phifixgendosedstar0}] simplifies to $P_k = 0$ if $k
< \bar{d^\ast}$ and $P_k = 1$ if $k \ge \bar{d^\ast}$. Substituting
these conditions into \Req{gcontlong.eq:phifixgendosedstar}, we
recover our previous expression \Req{gcontlong.eq:r=1T>1k>1} for the
$r=1$ homogeneous case.

\subsection{Universal classes of contagion}
\label{gcontlong.subsec:universal}
 
In the homogeneous version of the model, we determined the existence
of two classes of dynamics---epidemic threshold and critical
mass---with the former arising whenever $\bar{d^\ast}=1$, and the latter
when $\bar{d^\ast} \ge 2$.  In other words, in homogeneous populations,
the condition for differentiating between one class of behavior and
another is a discrete one.  Once heterogeneity is introduced, however, we
observe a smooth transition between epidemic threshold and critical
mass models governed by a continuous adjustment of the distributions
$f$ and $d$ (or equivalently the $\{P_k\}$).  One consequence of this now-continuous transition is the
appearance of an intermediate class of models, which we call
\textit{vanishing critical mass} (see
Fig.~\ref{gcontlong.fig:threeclasses}). As with pure critical mass
models, this new class is characterized by a saddle-node bifurcation,
but now the lower unstable branch of fixed points crosses the
$p$-axis at $p_c$; in other words, the required critical mass
``vanishes.''  The collision of the unstable branch of the saddle node
bifurcation with the horizontal axis also effectively reintroduces a
transcritical bifurcation, in the manner of epidemic threshold models.
However, the transcritical bifurcation is different from the one
observed in epidemic threshold models because the rising branch of the
bifurcation has negative, rather than positive, slope and comprises
unstable, rather than stable, fixed points.  Vanishing critical mass
dynamics are therefore qualitatively distinct from both previously
identified classes of behavior, exhibiting important properties of
each: for $p<p_c$, they behave like critical mass models; and for
$p>p_c$, they behave like epidemic threshold models, in the sense that
an infinitesimal initial seed can spread.

Our generalized model therefore exhibits behavior that falls into
one of only three universal classes: class I (epidemic threshold), class II
(vanishing critical mass), and class III (critical mass).  As we show
below, more complicated fixed point curves exist (i.e., curves possessing two or more saddle-node
bifurcation points) but nevertheless
belong to one of these three universal classes, since together they include all
possible behaviors of the fixed point curves near $p=p_c$.

\begin{figure*}[tbp!]
  \centering
  \includegraphics[width=0.98\textwidth]{fig007_figgc_lp_classes3_noname.pdf}
  \caption{
    Fixed point curves for the three main classes of contagion models produced
    by the model, along with transitions between classes.
    All curves are determined by numerical solution
    of~\Req{gcontlong.eq:phifixgendosedstar} with an error tolerance
    of $10^{-10}$.
    Here, $T=12$, individual doses are lognormally distributed with mean 1
    and a variance of $0.6$ for the underlying normal distribution,
    and thresholds are homogeneously distributed with
    $g(d^\ast) = \delta(d^\ast-\bar{d^\ast})$.
    The main plots correspond to 
    (A) $\bar{d^\ast}=0.4$ (class I),
    (B) $\bar{d^\ast}=1.5$ (class II),
    and
    (C) $\bar{d^\ast}=3$ (class III).
    The insets plots (D) and (E) respectively show fixed point curves for models
    at the class I-class II transition (i.e., when $P_1=P_2/2$ which occurs
    when $\bar{d^\ast}=0.8600\ldots$) and the
    class II-class III transition (i.e., when $p_c=1$ which occurs when 
    $\bar{d^\ast}=1.9100\ldots$).
    In all plots except (E), the intersection between the fixed point curve
    and the $p$-axis is a transcritical bifurcation.  In plots
    (B), (C), and (E), the second bifurcation is a saddle-node bifurcation.
    }
  \label{gcontlong.fig:threeclasses}
\end{figure*}

In addition to identifying three universal classes of behavior, it is
also possible to specify the conditions that govern into which class
any particular choice of model parameters will fall. For $r=1$, we can
calculate when the transitions between universal classes occur in
terms of the parameters of interest; that is, the $\{P_k\}$ (viz $f$ and $g$) and
$T$. This exercise amounts to locating the transcritical bifurcation
and determining when it collides with the saddle-node bifurcation.
To locate the transcritical bifurcation, we examine the fixed point equation as ${\phi^{\ast}}
\rightarrow 0$.  Since, from \Req{gcontlong.eq:phifixgendosedstar},
\begin{equation}
  \label{gcontlong.eq:phifix->0}
  {\phi^{\ast}} = 
  T p {\phi^{\ast}} P_1
  + O({\phi^{\ast}}^2),
\end{equation}
we have
\begin{equation}
  \label{gcontlong.eq:transc}
  p_c = \frac{1}{T P_1}.
\end{equation}
The position of the transcritical bifurcation is therefore determined
by the memory length $T$ and the fraction of individuals who are
typically infected by one exposure, $P_1$.  We see immediately that
the transition between class II and class III contagion models occurs
when $p_c=1$; that is, when
\begin{equation}
  \label{gcontlong.eq:transII-III}
  P_1 = 1/T.
\end{equation}

Recalling the homogeneous case, we see that when $\bar{d^\ast} = d =1$, $P_1=1$,
giving $p_c=1/T$ as before. It necessarily follows that when
$\bar{d^\ast}>d=1$, $P_1=0$ and therefore $p_c=\infty$, thus confirming our
earlier finding that the homogeneous model with $\bar{d^\ast}>1$ is always
in the pure critical mass class (i.e., the lower, unstable fixed point
curve of the saddle node bifurcation must have ${\phi^{\ast}}>0$ for all
$p_b \le p \le 1$ since it only reaches ${\phi^{\ast}}=0$ in the
limit $p \rightarrow \infty$).  We also see that technically all models possess
a transcritical bifurcation somewhere along the $p$-axis, even though
it may be located beyond $p_c>1$ (when $P_1=0$, it lies at $p_c=\infty$).

In order to locate the transition between classes I and II, we determine
when the transcritical and saddle-node bifurcations are coincident,
i.e., when $\dee{{\phi^{\ast}}}/\dee{p}=\infty$ at $(p,{\phi^{\ast}})=(p_c,0)$.  
In other words, we calculate
when the fixed point curve emanating from $p=p_c$ is at the transition
between having a large positive slope (class I) and a large negative
slope (class II). We find the condition for this first transition is
\begin{equation}
  \label{gcontlong.eq:transIcond}
  P_1 = P_2/2,
\end{equation}
where details of this calculation are provided in Appendix~\ref{gcontlong.app:transI-IIsec}.

The condition of \Req{gcontlong.eq:transIcond} is a
statement of linearity in the $\{P_k\}$, though importantly only for $k=1$ and $k=2$.
Providing $p_c<1$, if $P_1 > P_2/2$ (i.e., sublinearity holds) a
contagion model is class I whereas if $P_1 > P_2/2$ (i.e.,
superlinearity holds) it is class II.  This condition means that class II
contagion models arise when, on average, two doses are more than twice
as likely as one dose to cause infection.

When this condition is satisfied exactly,
the system's phase transition is a
continuous one as per those of class I, but of a different
universality class: class I systems exhibit a linear scaling near the
critical point whereas the scaling when $P_1=P_2/2$ is ${\phi^{\ast}}
\propto (p-p_c)^{1/2}$ (in Appendix~\ref{gcontlong.app:transI-IIsec},
we show that a sequence of increasingly specific
exceptions to this scaling exist depending on the extent of linearity present in the $\{P_k\}$).  

Thus, for $r=1$, the condition that determines whether a given
system is described by an epidemic threshold model, a critical mass model, or a member of 
the intermediate class of vanishing critical mass models, 
depends only on $T$, $P_1$, and $P_2$---a surprising result
given that \Req{gcontlong.eq:phifixgendosedstar} clearly
depends on all the $\{P_k\}$.
A summary of the three basic system types, the transitions
between them, and the accompanying conditions is given in
Table~\ref{gcontlong.tab:systemsummary}.

\begin{table}[tbp]
  \centering
  \begin{tabular}{l|l}
    \hline
    Model class: & Conditions: \\
    \hline
    I: Epidemic threshold &
    $P_1 > P_2/2$ and $P_1 > 1/T$ \\
    I-II transition & 
    $P_1 = P_2/2$ and $P_1 > 1/T$ \\
    II: Partial critical mass & 
    $P_1 < P_2/2$ and $P_1 > 1/T$ \\
    II-III transition & 
    $P_1 < P_2/2$ and $P_1 = 1/T$ \\
    Pure critical mass &
    $P_1 < P_2/2$ and $P_1 < 1/T$ \\
    \hline
  \end{tabular}
  \caption{
    Summary of basic states of the $r=1$, heterogeneous version of model along
    with the corresponding parameter ranges.
    The $P_k$ quantities depend on the distributions of dose size
    and individual thresholds, and are given in~\Req{gcontlong.eq:phifixgendosedstar0}.
  }
  \label{gcontlong.tab:systemsummary}
\end{table}

For $r<1$, the position of the transcritical bifurcation,
\Req{gcontlong.eq:transcr}, generalizes in the same manner as
\Req{gcontlong.eq:r<1k1T>1pc2}. We find
\begin{equation}
  \label{gcontlong.eq:transcr}
  p_c = \frac{1}{P_1(T + \tau)},
\end{equation}
where we recall that $\tau = 1/r - 1$.  As we will see 
in Section~\ref{gcontlong.sec:SIR}, the above statement is
also true for $\rho<1$ and thus stands as a completely general
result for the model.  We therefore have a condition for the
transition between class II and class III contagion models for a given
$T$ and $r$, analogous to that for the $\rho=1$ case
in~\Req{gcontlong.eq:transII-III}: $P_1 = 1/(T+\tau)$.  For $r<1$, the
condition for the transition between class I and class II contagion
models is more complicated both in derivation and form.  We observe
that as $r$ is decreased, class III models must at some point
transition to class II models and class II models
will eventually become class I models, where we can 
determined the former tranisition by setting
$p_c=1$ in \Req{gcontlong.eq:transcr} and solving for $r$:
\begin{equation}
  \label{eq:gcontlong.eq:transII-IIIr}
  r = \frac{P_1}{1+P_1 - P_1T}.
\end{equation}
In \Req{eq:gcontlong.eq:transII-IIIr} we have assumed $P_1>0$ and $P_1T<1$, since the $r=1$ limit is, by assumption, 
a class III contagion model.
Fig.~\ref{gcontlong.fig:gc_T3_k2k1_bif_theorycomp} presents an
example of a model transitioning from Class III through Class II to Classs I.

\begin{figure}[tbp!]
  \centering
  \includegraphics[width=0.49\textwidth]{fig008_figgc_T3_k2k1_bif_theorycomp2b_002_noname.pdf} 
  \caption{
    Example of transitions between all three classes of
    contagion models occurring as the probability of recovery $r$ changes.
    Shown are theoretically derived fixed point curves for a range of $r$
    for $T=3$ and a heterogeneous population with threshold distribution 
    $g(d^\ast) = 0.3 \delta(d^\ast-1) + 0.7 \delta(d^\ast-2)$.
    From left to right, the curves correspond to $r=0.01, 0.05, 0.10, 0.15, \ldots, 1.00$.
    Solid lines indicate stable fixed points,
    dashed lines unstable ones, and the squares mark saddle-node bifurcation points.
    The $p$-axis is a line of fixed points,
    and for all fixed point curves, dashed lines indicate unstable fixed points
    and solid lines stable ones. 
    For $r>0.75$, the system belongs to class III (see \Req{eq:gcontlong.eq:transII-IIIr}).
    Fixed points lying on the $p$-axis are not indicated in the plot.
    For each instance of the model, points on the $p$-axis
    to the left of $p_c$ (the intercept of the non-zero fixed point curve with
    the $p$-axis) are stable and points to the right are unstable.
    The expression used here to determine the fixed points is formed by
    appropriately weighting a combination of
    Eqs.~\req{gcontlong.eq:r<1k1T>1}, \req{gcontlong.eq:Gamma2}, and \req{gcontlong.eq:Gammafull}.
  }
  \label{gcontlong.fig:gc_T3_k2k1_bif_theorycomp}
\end{figure}


\subsection{Composite classes of dynamics}
\label{gcontlong.subsec:composite}

\begin{figure}[tbp!]
  \centering
  \includegraphics[width=0.49\textwidth]{fig009_figgc_manybif_thy_comb04_noname.pdf} 
  \caption{
    Fixed point curves involving multiple bifurcations.  In all four examples,
    $T=20$ , $r=1$, and dose size is uniformly held at unity.
    The different curves are obtained by adjusting the threshold distribution
    $g$ which in turn leads to changes in the $\{P_k\}$ [see~\Req{gcontlong.eq:phifixgendosedstar0}].
    The curves correspond to 
    (A) $g(d) = 0.2\delta(d-1) + 0.8\delta(d-6)$;
    (B) $g(d) = 0.15\delta(d-1) + 0.4\delta(d-5) + 0.45\delta(d-12)$;
    (C) $g(d) = 0.075\delta(d-1) + 0.4\delta(d-2) + 0.525\delta(d-12)$;
    and
    (D) $g(d) = 0.3\delta(d-3) + 0.7\delta(d-12)$.
    Note that example (D) consists of two separate fixed point curves.
    The curves were found numerically by solving for $p({\phi^{\ast}})$ 
    using~\Req{gcontlong.eq:phifixgendosedstar}.  
    Solid and dashed lines indicate stable and unstable fixed points respectively.
  }
  \label{gcontlong.fig:manybif}
\end{figure}

Although our main results (three universal classes of behavior and
the conditions that govern the transitions between the three) involve
examples of contagion models with at most two bifurcations, other more
complicated kinds of equilibrium behavior are possible.
Fig.~\ref{gcontlong.fig:manybif} shows four examples of what can
happen for particular distributions of $d^\ast$ across the
population. In each example, $T=20$ and $r=1$, all doses are of unit
size, and the population is divided into either two or three
subpopulations with distinct values of $d^\ast$. The main features of
each system are captured by the number and locations of the
saddle-node bifurcations. As was the case for the homogeneous model [see
\Req{gcontlong.eq:r=1T>1k>1diff2}], we are able to find an expression
for $z=p_b{\phi_b^{\ast}}$:
\begin{equation}
  \label{gcontlong.eq:gdiff_eq}
  0 
   =  \sum_{k=1}^{T}
  \binom{T}{k}
  P_k
  z^{k-2}
  (1-z)^{T-k-1}
  [k-1 -z(T-1)],
\end{equation}
where the details of this calculation are provided in
Appendix~\ref{gcontlong.app:saddlenodeextra}. In principle,
\Req{gcontlong.eq:gdiff_eq} could be analyzed to deduce which
$\{P_k\}$ (and hence which $f$ and $g$) lead to
what combination of bifurcations.  
While substantially more complicated, and thus beyond the
scope of this paper, such a classification scheme would be
a natural extension of our
present delineation of the model into three universal classes of
contagion models, based on the behavior near $p=p_c$.  One simple
observation, however, is that \Req{gcontlong.eq:gdiff_eq} is a polynomial of
order $T-2$ and hence a maximum of $T-2$ saddle-node bifurcations may
exist.  From our investigations this outcome seems unlikely as nearby
bifurcations tend to combine with or overwhelm one another;
in other words, subpopulations with sufficiently distinct $d^\ast$ are required to
produce systems with multiple saddle-node bifurcations.  Furthermore,
the distribution of dose sizes $f$ seems unlikely to be
multimodal for real contagious influences or entities,
and the way it enters into the calculation of the $\{P_k\}$
[\Req{gcontlong.eq:phifixgendosedstar0}] reduces its effect in
producing complicated systems.  Thus, the number of distinct bifurcations
is limited and strong multimodality in the threshold distribution $g$
appears to be the main mechanism for producing systems with more
than one saddle node bifurcation.
As a first step in this extended analysis of the model, we derive 
in Appendix~\ref{gcontlong.app:saddlenodeextra} the condition
for the appearance of two saddle-node bifurcations (i.e., one forward
and one backward).

\section{SIRS and SIR Contagion Models}
\label{gcontlong.sec:SIR}

As mentioned in section~\ref{gcontlong.sec:model}, the SIS
class of behavior that we have analyzed exclusively up to now is a
somewhat special case of the general contagion process as it assumes
that recovered individuals instantly become re-susceptible.  This
assumption renders the SIS case particularly tractable, and we have
taken advantage of this fact in the preceding sections to make some
headway in understanding the full range of equilibrium behavior of the
model.  However, it remains the case that very few, if any, infectious
diseases could be considered to obey true SIS-type dynamics, as almost all
recovery from infection tends to be associated with some finite period
of immunity.  Any purportedly ``general'' model of contagion ought
therefore to be analyzed in a wider domain of the associated removal
period, and any corresponding classes of behavior labeled ``universal''
ought to withstand the introduction of at least some period of
immunity to re-infection. Thus motivated, we now extend our previous
analysis to systems where individuals experience temporary (SIRS, $0 <
\rho < 1$) or permanent removal (SIR, $\rho=0$). We present
some preliminary results for each of these cases in turn, relying now exclusively on 
numerical simulation.

\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.49\textwidth]{fig010_figgcSIR_logn_070b_002_noname.pdf}
  \caption{
    The effect of reducing $\rho$ (the probability of an immune individual becoming
    susceptible) for a contagion model that is class II when $\rho=1$.
    Here, $N=10^5$, $T=6$, $r=1$, dose sizes are fixed at unity,
    and $g(d^\ast) = 0.25 \delta(d^\ast-1) + 0.75 \delta(d^\ast-2)$.
    From left to right, the plots correspond to 
    $\rho=1$, $0.15$, and $0.09$.
    In all cases, $p_c$, the position of the $p$-axis intercept, is independent of $\rho$.
    Apart from the intercept at $p=p_c$ and the point $(p,{\phi^{\ast}})=(1,1)$, all other points of the curve move
    to the right as $\rho$ decreases.
    Consequently, as $\rho$ is reduced,
    all class II contagion models will at some point become members of class I.
    Contagion models in class I and class III remain unchanged in their nature.
    In finding the upper stable branch of ${\phi^{\ast}}(p)$, all individuals are initially infected.
    A binary search is then used to detect the position of the lower unstable branch.
  }
  \label{gcontlong.fig:rhoclasstrans}
\end{figure}

For SIRS contagion, we observe that the position of the transcritical
bifurcation $(p_c,0)$ does not change as $\rho$ is reduced from 1;
however, all non-zero fixed points move in the positive $p$ direction.  
Because they remain in the removed state for a longer time, individuals in systems with lower $\rho$
spend relatively less time infected than those in systems with higher $\rho$
(this is in contrast to the effect of reducing $r$, 
which prolongs the time individuals are infected,
thereby causing fixed points
to move in the negative $p$ direction).
Thus, contagion models
belonging to class I and class III in the $\rho=1$ special case
remain in their respective  classes as $\rho$ decreases.  
However, as shown in Fig.~\ref{gcontlong.fig:rhoclasstrans},
class II models will transition to class I for some $\rho < 1$.


\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.49\textwidth]{fig011_figgcSIR_logn_050comb3_002_noname.pdf}
  \caption{
    Time taken for an initially universal infection to die out in the SIR case ($\rho=0$).
    Moving from bottom to top in both plots, $T=3$, 4, 5, and 6.
    The model parameters used here are $d^\ast=1$, unit dose size, $r=0.2$,
    and population size $N=10^4$.   For all systems, $p=\phi=1$ is a
    fixed point and hence $t_{\rm max}=\infty$ at $p=1$.
    The main plot shows the rapid increase in $t_{\rm max}$ for $p \rightarrow 1$.
    The inset shows the behavior of $t_{\rm max}$ as a function of  $(1-p)$
    as plotted on a double logarithmic scale.
    The lines have slopes of $-(T-1)$ indicating that infections are
    strongly persistent with $t_{\rm max} \propto (1-p)^{-(T-1)}$.
  }
  \label{gcontlong.fig:SIR_logn_050}
\end{figure}

For the SIR version of the model, recovered individuals cannot return
to the susceptible state, and no fraction of the population remains
infected in the infinite time limit; hence we can no longer speak of
non-zero fixed point curves, and the fraction of individuals infected and
recovered, $\phi(t)$ and $R(t)$, become
the relevant objects of study.  Some defining quantities are then
the maximum fraction of the population infected at any one time, $\max_t \phi(t)$, the fraction
eventually infected, $1-\phi(\infty)-R(\infty)$, and the \textit{relaxation time}
required for the epidemic
to die out.  We focus on the latter here which we denote by $t_{\rm max}$.

We observe that when individuals possess a memory of doses (i.e., $T>1$), $t_{\rm max}$
diverges as $p \rightarrow 1$.  
Figure~\ref{gcontlong.fig:SIR_logn_050} shows $t_{\rm max}(p)$
for four sample systems with $T=3$, 4, 5, and 6.  We find the
divergence of $t_{\rm max}$ near $p=1$ to be well approximated by
\begin{equation}
  \label{gcontlong.eq:SIR_tmax}
  t_{\rm max} \propto (1-p)^{-(T-1)},
\end{equation}
where fits are shown in the inset of Fig.~\ref{gcontlong.fig:SIR_logn_050}.
Equation \req{gcontlong.eq:SIR_tmax} also 
shows that for a fixed $p$, the relaxation time increases exponentially with 
length of memory $T$ for all $p$, i.e., 
\begin{equation}
  \label{gcontlong.eq:SIR_tmax2}
  t_{\rm max} \propto e^{\alpha T},
\end{equation}
where $\alpha=-\ln{(1-p)}>0$.

Thus, when $\rho=0$, epidemics,
while not ever achieving a non-zero steady state as in the $\rho>0$ case,
can persist for (arbitrarily) long periods of time.
The introduction of memory, which allows
infected individuals to maintain their dose count 
above their threshold by repeatedly infecting 
each other, creates an SIR model with strikingly different
behavior to the standard memoryless SIR model.

\section{Concluding remarks}
\label{gcontlong.sec:concl}

Our aim here has been to develop and analyze in detail a model of
contagion that incorporates and generalizes elements of contagion models from
both the social sciences and epidemiology.
A key feature of the current model is
that interdependencies between successive exposures are introduced in
a natural way by varying the length of memory that an individual
maintains of past exposures.  Contagion models incorporating memory correspond to
standard notions of contagion in the social sciences (although these
models rarely discuss the role of memory explicitly), while
memory-less contagion models correspond to traditional models of
disease spreading.  Our model suggests, however, that in reality these
two kinds of contagion models may not be entirely distinct.  We
venture the possibility that some infectious diseases may spread in a
fashion similar to social contagion processes.  For example,
two exposures to an agent sufficiently close in time may infect an
individual with higher probability than would be 
expected if the exposures acted as independent events.
The response of some individuals to
allergens where separate doses accumulate in the body may be an
example of a disease with memory.  Although allergies are not
contagious, they demonstrate that such a phenomenon is biologically
plausible---a possibility that has largely remained unexamined in the
microorganismal dose-response literature~\citep{haas2002}.

The main result of our analysis is the identification of three
universal classes of contagion dynamics, along with precise conditions
for the transitions between these classes.  Given the
complexity of the model, these conditions are surprisingly simple---at
least in the SIS case ($r=\rho=1$, see
Table~\ref{gcontlong.tab:systemsummary})---and present us with
quantities such as $T$, $P_1$, and $P_2$ that may in principle be
measurable for real epidemics.  Furthermore, the dependence of the
transition conditions only on $P_1$, $P_2$, and $T$, rather than on the full
details of the underlying distributions of thresholds $(g(d^\ast))$
and doses $(f(d))$, suggests a new and possibly useful level of
abstraction for thinking about contagious processes; that is,
measuring individuals in terms of their dose-response and
characterizing a population in terms of its $\{P_k\}$.

For the more complicated and general cases of $r<1$ (finite recovery
period) and $\rho<1$ (finite immunity period), we have confirmed that
the same basic three-class structure persists, and determined the
position of the transcritical bifurcation $p_c$,
\Req{gcontlong.eq:transcr}, that is one of the two quantities needed
to specify the conditions for transitioning between classes. The other
condition, derived from calculating the slope of the fixed point curve
as it passes through the transcritical bifurcation, merits further
attention.  All of these conditions, however, are ultimately dependent
on the distributions $f$ and $g$.  Our analysis of the model suggests,
for example, that composite fixed point diagrams involving more than
one saddle point node can only result from multimodality in $g$, the
distribution of thresholds. Exactly how the details of these two
distributions affect $P_1$ and $P_2$ would be worth further
investigation.

Our model suggests that some epidemics may be prevented or enabled
with slight changes in system parameters (if feasible). For example,
knowing that a potentially contagious influence belongs to class II and that $p$
is just below $p_c$ would indicate that by increasing inherent
infectiousness (i.e., $p$) or by creating a sufficiently large enough
base of infected individuals, the contagion could be kicked off with
potentially dramatic results. Alternatively, by increasing $r$ or
reducing $T$ or $\rho$, the possibility of undesirable epidemics may be reduced,
as for all these adjustments fixed point curves are generally moved in the
direction of higher values of $p$.

Finally, we have focused exclusively in this paper on a
 mean field analysis of the model, which
is to say, we have made the standard
assumption that individuals in a population mix uniformly and at
random.  A natural generalization would
be to consider the model's behavior for a networked population of
individuals.  Other simulation possibilities would be to consider
distributions of $T$, $p$, $r$, and $\rho$. Finally, additional technical
investigation of this model would also include analysis of the 
closed-form expression for saddle-node bifurcations given in
\Req{gcontlong.eq:gdiff_eq}, and a derivation of analytic expressions
for the model when $\rho<1$ for small $T$ and $d^\ast$.  We hope that
our preliminary investigations into this interesting and reasonably
general class of contagion models will stimulate other researchers to
pursue some of these extensions.

\acknowledgments

We are grateful for discussions with Duncan Callaway, Daryl Daley, Charles Haas,
and Matthew Salganik.
This research was supported in part by the National
Science Foundation (SES 0094162), the Office of Naval Research, Legg
Mason Funds, and the James S.\ McDonnell Foundation.

