Stochastic epidemic-type model with enhanced connectivity: exact solution
We present an exact analytical solution to a one-dimensional model of the Susceptible-Infected-Recovered (SIR) epidemic type, with infection rates dependent on nearest-neighbor occupations. We use a quantum mechanical approach, transforming the master equation via a quantum spin operator formulation. We calculate exactly the time-dependent density of infected, recovered and susceptible populations for random initial conditions, and compare our results with a low connectivity SIR model reported by Schütz et al. schutz (). Our results compare well to those of previous work, validating the model as a useful tool for additional and extended studies in this important area. Our model also provides exact solutions for the n-point correlation functions, and can be extended to more complex epidemic type models.
The study of cooperative evolution of multi-agent systems is a part of many sciences, from biology and social science to physics, chemistry and engineering. The methods of statistical physics are often employed, and simple models are building blocks in this quest to understand complexity. From the point of view of non-equilibrium statistical physics, biology is an exciting area of investigation. After all, every living organism is an example of a far-from-equilibrium system, and stochastic processes are ubiquitous in biological systems.
Epidemic-type models abound in the literaturebook reviews (), original SIR (), from very simple ones that capture the basic rules of the infection mechanism, to very complex models that account for spatial spread, age structure and the possibility of immunization complex models (). It is interesting to see how some of these epidemic-type models have been applied successfully in other fields as well, such as social sciences (voter models, rumor spreading models) liggett (),privman () or computer science (the spread of a virus in a computer network) computer virus (). Some models are deterministic, following a set of evolution equations with given initial conditions solved using the mean field theory approach, while others are stochastic and studied using methods such as the Langevin equation, the Fokker-Planck equation and computer simulations. Some recent numerical studies of epidemic-type models can be found in numerical studies ().
Despite numerous studies and approximation schemes, exact solutions for epidemic models are rare. A study that sparked our interest was published in 2008 by Schütz et al. schutz (). This presents an exact solution for a stochastic one-dimensional SIR (susceptible/infected/recovered) epidemic model. The method used is a quantum mechanical formulation of the master equation in terms of second quantized operators. The authors define cluster functions that describe the behavior of susceptibles adjacent to infected individuals at the cluster boundaries. They derive and solve exactily a set of coupled evolution equations for these functions. Their exact solution shows the significant difference between low connectivity and high connectivity SIR models, and the role of fluctuations. Fluctuations are built into the exact solutions, but are missing in the mean field approach.
Working towards an exact solution for an SIR model with higher connectivity, we study a variation of a one-dimensional SIR model in which we define different rates of infection depending on the number of infected neighbors. A susceptible individual with two infected neighbors will have a different probability of being infected than a susceptible neighboring only one infected person. As in the traditional SIR model, we assume the possibility of recovery. (To differentiate it from other one-dimensional SIR models, we refer to our model as the dual neighbor model.) This model also can be cast as a non-conservative voter model, with the three classes of individuals defined as S - undecided, I - biased, R - decided.
We here present an exact solution of the dual-neighbor SIR model. We employ a quantum mechanical approach to the problem, using the cluster function method used by Schütz schutz (). The steady-state solution depends on initial populations of susceptible and infected individuals. It has fluctuations built-in, and has a stationary state different from the low connectivity SIR model. Although the overall trend of the solution is similar to that of the low connectivity model, there are significant differences as well.
In Section 2, we define our model and its quantum mechanical representation. Next, we present the cluster function method and derive the evolution equations for the cluster functions and particle densities that fully resolve the model (Section 3). We conclude with an analysis of our solutions, summarize of our work, and suggest some interesting open questions (Section 4).
Ii Dual Neighbor Model and its quantum mechanical representation
The traditional SIR model consists of a fixed number of individuals N split into three classes: susceptible, infected, and recovered. In a high connectivity (mean field) model, each ”node” of the network is represented by an individual in one of the three classes, in contact with every other node. A susceptible becomes infected with rate when in contact with an infected; an infected individual recovers spontaneously with rate ; recovered individuals cannot change. At a particular time the average number of individuals in each class is represented by , , and , with .
The time evolution of these classes is governed by a set of coupled differential equations:
The first equation describes reduction of via infection; the second shows increasing via infection of susceptibles, and decreasing by spontaneous recovery; and the third shows the increase of due to spontaneous recovery. This system of coupled nonlinear differential equations can be solved numerically, yielding the time dependence of each class of individuals, but without any information regarding correlations between individual nodes. This model exhibits smooth time dependence of the class populations, without any statistical fluctuations.
In contrast, we propose a stochastic one-dimensional model in which each node is in contact with only two other nodes (as if arrayed along a line), and the rate of infection depends on the number of infected neighbors. Representing linear sequences of neighboring individuals via strings of symbols (e.g. representing a susceptible to the left of an infected to the left of another susceptible), the dynamics of this model is defined as:
The first process describes the mechanism of infection of a susceptible neighboring two infecteds; the next two processes describe infection with a different rate when the susceptible has only one infected neighbor; and the last process describes spontaneous recovery with yet another rate. We can pick the values of , and based on the application of the model. To model actual epidemics, for example, we would assign a higher rate of infection when a susceptible is in contact with two infecteds compared to just one, thus . Note that our model does not allow a succeptible individual with a recovered on one side and an infected on the other ( or ) to become infected. As a model of disease transmission this represents an immunity provided by a recovered neighbor; in a voter model, it suggests that an undecided voter remains so if flanked by one biased voter and one who is decided. Practically, it says that the rate of infection in and configurations are negligible compared to other infection rates and . Making such an assumption produces a Hamiltonian amenable to exact solution of the rate equations. As we shall exhibit, features of the results of this model validate it as a reasonable first approximation to real-world situations.
The related mean field (deterministic) model is governed by the following differential equations for the average populations:
Compared with the traditional SIR model, the dynamics of this system is governed by three-point interactions between the S and I type individuals. This system also can be solved numerically for the time dependence of , and .
Proceeding beyond the mean field approximation, the time evolution of our model is best described by the master equation, expressing conservation of probability within a continuous-time dynamics. We let represent a configuration, giving the state (, , ) of each of the individuals. Proximity is defined by labeling each individual with an index and assuming individual is adjacent to individual (periodic boundary conditions.) The master equation expresses the rate of change of the probability of finding the system in configuration at time as the rate of transfer of probability into from other configurations less the rate at which passes probability into others glauber ():
The transition rate is the probability per unit time that configuration changes into a different configuration .
Utilizing Dirac notation, we represent a configuration as , and use standard methods schutz1 () to build a vector representation for a general state as a probabilistic superposition of all configurations of a system:
where is the probability that the system will be found in configuration at time . This allows the master equation to be re-written as
where the pseudo-Hamiltonian has matrix elements:
A formal solution to Eq. 3 can be written as . In this formalism, the expectation value (at time t) of a physical quantity that has value for configuration is
where is the operator corresponding to the observable and the state is the sum (with weight 1) of all configurations. Time dependence of obeys
where is the commutator of the operator with the Hamiltonian.
The operators for our model can be written in terms of creation and annihilation operators, following the rules for constructing a quantum Hamiltonian as presented in schutz (). We define the , to be the creation and annihilation operators for a type particle at site , and , to be the creation and annihilation operators for a type particle at site . A site with neither an infected nor a susceptible is assumed to contain a recovered. Operators at different sites commute, and operators of the same species at the same site anticommute. The number operator () can only take values 0 and 1, and its expectation value equals the density of type () particles. In this notation,
We can rescale the time variable in such a way that the rate (assumed non-zero) becomes one. With that assumption, and defining and , the Hamiltonian obeys
Iii Exact solution: derivation and analysis
The Hamiltonian function (Eq. 8) is invariant relative to translation along the string of sites, so if we assume a translationally invariant initial state, all future configurations also have this property. Translation invariance makes these clusters independent of their starting point . Because a susceptible cannot change into an infected unless they are neighbors, clusters change only at their edges.
The equations for the cluster functions are derived by calculating the respective commutators of the cluster operator products with the pseudo-Hamiltonian .
For , this gives
The equations of motion for the particle density of susceptibles and infecteds are found in a similar fashion:
In order to solve for the particle density of the susceptibles and infecteds, we thus need first to find the cluster functions.
iii.1 Solution for cluster functions
Introducing an operator defined by the property we can rewrite Eq.14 (for ) as
which has a formal solution
Representing the initial density of susceptibles and infecteds as and , we describe an initial state such that . From this,
Eq. 12 for G(1)can be rewritten as
and can be integrated by standard means to yield
Using a similar method, we find the solution for to be
where in order to express the solution more compactly, we have defined a relaxation time .
iii.2 Particle densities
The pieces are in place now to find the time dependence of the average number of susceptibles and infecteds. We assume that at there are no recovered individuals, therefore . Combining Eq.’s 15, 18 and 19 we find
where we have defined a second relaxation time . Integration yields
where the integration constant turns out to be the stationary value of the density of susceptibles for large times. It is evaluated by enforcing :
This is non-zero because a susceptible can be ”trapped” with a recovered on either side, and be no longer susceptible to infection. For the density of the infecteds, we combine Equations 16, 18 and 19 to get:
where is an integration constant to be evaluated using :
Note that the steady-state value for is zero, due to the process of spontaneous transformation of an infected into a recovered, regardless of its neighbors.
The time dependence of the average number of recovereds can be easily found from
iii.3 Analysis of solutions
In an actual epidemic, major concerns are the time evolution of the number of infected individuals, and their maximum number. Basic questions include: at what time is the peak of the infection reached; what are the factors that control the maximum number of infecteds; and how long does it take for a population to fully recover? With these questions in mind, we examine the exact solution for the infected population, Eq. 24. The solution depends on the initial fraction of infecteds, , the parameter representing the relative rate of infection with two infected neighbors vs. that for one infected neighbor, and representing the relative rate of spontaneous recovery vs. the rate of infection with a single infected neighbor.
We point out two values of that deserve special attention, and . Physically, means that the the probability of infection is the same regardless of the number of infected neighbors. This most closely matches the previous SIR model schutz (). The case corresponds to an infection rate proportional to the number of infected neighbors, perhaps most closely aligned with the case of a more realistic three-dimensional medical infection model.
Two typical time behaviors of the density of infecteds are shown in Fig. 1. Figure 1(a) shows the density rising from an initially small value to a maximum, and thereafter relaxing exponentially via spontaneous recovery until all infected individuals have disappeared. If the initial infection rate is high and/or the spontaneous recovery rate is high compared to infection rates, however, the model can exhibit time behavior that shows no peak, as in Fig. 1(b).
Figures 2(a) and 2(b) show how infection peak height and peak time vary with the initial state infection density. The peak height shows an increase as increases, leveling off as it must as one approaches an initial state with the majority of the population already infected. The smaller the initially infected population, the longer delay until the maximum number of infecteds is experienced. For values close to one, the peak will be early, if it appears at all. Because we have rescaled time in such a way that the infection rate of a susceptible with a single infected neighbor is , the proportional delay of the peak as a function of this single neighbor infection rate is hidden in these graphs.
Figure 3 shows how the spontaneous recovery rate controls whether a peak in appears or not. As expected, for a high spontaneous recovery rate the infected population begins to decay immediately once recovery is ”turned on,” as might be the case at the beginning of an immunization effort in the population. As the recovery rate drops, the figure shows the appearance of a peak that grows in size and in delay time.
Figure 4 shows the comparison of our solution (continuous curve) with results of the low connectivity model obtained by Schütz et al. (dashed curve) and the mean field approximation (dotted curve) for two representative sets of parameters. There are noticeable differences among the three solutions when becomes large in comparison to (Fig. 4(a)). These differences disappear for values on the order of or smaller than , as shown in Fig. 4(b). We notice that the peak of the infection in the mean field approximation is higher and happens at a later time than the one predicted by our model.
The various behaviors exhibited in the figures well represent effects seen in actual situations of disease spread mollison (). This serves to commend the model as a reasonable approach to modeling such behavior, and the exact analytical solutions enable precise predictions about features of the solutions.
In this paper we analyzed a modified one dimensional stochastic SIR model, the dual neighbor model, with infection rates dependent on the nearest-neighbor occupation. Using a quantum mechanical approach and the cluster function method introduced in Schütz et al. schutz (), we found an exact solution for the mean densities of susceptibles, infecteds and recovereds as a function of time and initial conditions. The quantum mechanical approach is a powerful analytical tool, because it leads to exact solutions not only for the mean number of infecteds, susceptibles and recovereds, but also for n-point correlation functions.
We analyzed our solution for the density of infecteds in various parametric regimes, and also compared our solution with the low connectivity model reported in schutz (). The qualitative behaviors exhibited by this model are straightforward and not unexpected. On the other hand, the details that this solution presents could lead to proposals for moderating the effects of an epidemic. This suggests the utility of further investigation of similar models, perhaps in two dimensions.
Although this model was introduced as an epidemic model, it can also be extended for other areas of study as well, such as voter problems, computer virus dynamics, or surface deposition. It can also be generalized to include correlated initial conditions.
We hope to further our study to find exact solutions for more complicated (and realistic) SIR models that include time delay of infection, the possibility of immunization, and also the reappearance of the disease (R can become S). Unfortunately the cluster function method fails for this latter (SIRS) model, and other mathematical avenues must be pursued.
I. M. wants to express special thanks to the Kavli Institute for Theoretical Physics for hospitality and financial support. This research was supported in part by the National Science Foundation under Grant No. PHY05-51164.
- (1) Mathematical epidemiology, ed. F. Brauer (Springer, 2008); O. Diekmann, J. A. P. Heesterbeek: Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation (Wiley Series in Mathematical and Computational Biology, 2000); J. D. Murray, Mathematical Biology I: An Introduction, (Springer-Verlag, Berlin 2002.)
- (2) W. O. Kermack and A. G. McKendrick, Proc. R. Soc. Edinburgh, 115: 700 (1927.)
- (3) P. Yongzhen et al. Applied Mathematical Modelling. 35: 3866 (2011); T. K. Kar and A. Batabyal, Biosystems, 104: 127 (2011); S. Zou et al, Phys. Rev. E , 83: 056121 (2011.); M. E. J. Newman, Phys. Rev. E , 66: 016128 (2002.)
- (4) T. M. Liggett Stochastic Interacting Systems: Contact, Voter, and Exclusion Processes (New York: Springer 1999.)
- (5) Nonequilibrium statistical mechanics in one dimension, ed. V. Privman (Cambridge University Press, Cambridge, England, 1997.)
- (6) J. C. Piqueira, and V. O. Araujoa, Applied Mathematics and Computation 213: 355 (2009.)
- (7) David R de Souza et al. J. Stat. Mech. Theor. Expr., P03006 (2011); T. Tome and R. M. Ziff Phys. Rev. E , 82: 051921 (2010); T. Tome and M. J. de Oliveira, J. Phys. A: Math. Theor., 44: 095005 (2011.)
- (8) G.M. Schütz, M. Brandau and S. Trimper, Phys Rev. E, 78: 061132 (2008.)
- (9) R. J. Glauber, J. Math. Phys. 4: 294 (1963); M. Doi, J. Phys. A , 9: 1465 (1976.)
- (10) I. Peschel, V. Rittenberg, and U. Schulze, Nucl. Phys. B, 430: 633 (1994.)
- (11) G. M. Schütz, in Phase Transitions and Critical Phenomena, Vol. 19, ed C. Domb and L. Lebowitz (Academic Press, London, 2001.)
- (12) D. Mollison Mathematical Biosciences , 107: 255 (1991.)