The Azimuth Project
Blog - stochastic epidemic-type-1D-models

This is a blog article in progress, written by Jacob Biamonte. See also the network theory series. To see discussions of the article as it was being written, visit the Azimuth Forum. For the final polished version, go to the Azimuth Blog.

If you want to write your own article, please read the directions on How to blog.

For those of you following the network theory series, you already know we’ve been talking a lot about how process change with time. For example, how a virus might spread, or how a population of fish might depend on time. In most, though not all, of our posts we have assumed that “everything interacts with everything else” all of the time. What do I mean by this? Say you have a large population of rabbits, and a predator (a wolf), we typically assume that every rabbit in this population has an equal ‘opportunity’ to encounter the predator and be eaten. In other words, the predator is interacting with all of the rabbits all of the time. We know in real life this is not the case. The predator will only interact with a small number of rabbits at any given time: this depends heavily on spatial location.

Of course, throwing away spatial degrees of freedom is just one of many approximations that can be made. And as you’d expect, there are clearly many situations when these spatial degrees of freedom are not very relevant to a problem at hand. For example, the classic thought example is if you stir the contents of a test tube fast enough so that all of the contents could be considered as interacting all of the time. Such a situation is typically referred to as a ‘mean-field’ approximation.

Today we don’t want to talk about the mean-field approximation, since we’ve done a lot of that already. What we’ll do today is add physical location to our models. So far, we have not done this when considering stochastic Petri nets (or in the words of a chemist, reaction networks). We have however considered stochastic walks on graphs, for example in Part 16. A graph of course could be interpreted as depicting a physical situation, in which a system changes state on a 2D plane.

Today what we are going to do is study epidemics. We’ve done this before, through the lens of a commonly used model, called the SIR model (see Part 3). I hope you’ll agree that epidemics are no laughing matter, but that does not stop me from telling bad jokes.

  • SIR model is perhaps the only model important enough to have been knighted.

In contrast to what we’ve done before, today we assume that the epidemic will spread in 1D space. The way I like to think about this is all about a bad pet shop. I like to think of a line of rabbits in cages. On any given day, say that most of the rabbits start out healthy but susceptible to catching a cold from an infected rabbit. If none of the rabbits in any of the cages are sick, then things go on just fine. However, if one of the rabbits is sick, then there is a chance that this rabbit could infect a neighbouring rabbit. And the epidemic could spread further through the chain of cages. Today we are going to talk more about this situation, in more detail and talk about some recent work which has proposed some methods of exact solution for just these exact sort of scenarios.

Stated in another way, in this article, we’re going to explain how a member of a bunny chain like this

Susceptible rabbits SIR

can get sick like this

Single infected rabbit SIR

and then how stochastic mechanics takes over and infects the whole batch like so:

infected rabbits SIR

If you’re worried, please don’t, since they can recover:

Susceptible rabbits SIR

Susceptible rabbits SIR

SIR model
SIR model
SIR model

The recent work

Recently Blake Stacy pointed out the following recent preprint.

  • [WMM11] H. Thomas Williams, Irina Mazilu and Dan Mazilu, Stochastic epidemic-type model with enhanced connectivity: exact solution. (2011). arXiv:1108.5135v1

Abstract. 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 Schuetz et al.. 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 paper [WMM11] has clear relevance to the types of things that have been talked about in the network theory series. As mentioned, the SIR model was considered by John Baez in the petri net field theory framework in Network Theory Part 3.

The purpose of this page is to cast this recent work [WMM11] into our framework, and see what we can learn. We will label our sections precisely as was done in [WMM11] and conduct our analysis of the paper sequentially.

In the introduction, the authors motivate the problem and cite the following reference.

  • G.M. Schutz, M. Brandau and S. Trimper, Exact solution of a stochastic susceptible-infectious-recovered model, Phys Rev. E, 78: 061132 (2008). 10.1103/PhysRevE.78.061132.

Abstract. The susceptible-infectious-recovered (SIR) model describes the evolution of three species of individuals which are subject to an infection and recovery mechanism. A susceptible S can become infectious with an infection rate β by an infectious I type provided that both are in contact. The I type may recover with a rate γ and from then on stay immune. Due to the coupling between the different individuals, the model is nonlinear and out of equilibrium. We adopt a stochastic individual-based description where individuals are represented by nodes of a graph and contact is defined by the links of the graph. Mapping the underlying master equation onto a quantum formulation in terms of spin operators, the hierarchy of evolution equations can be solved exactly for arbitrary initial conditions on a linear chain. In the case of uncorrelated random initial conditions, the exact time evolution for all three individuals of the SIR model is given analytically. Depending on the initial conditions and reaction rates β and γ, the I population may increase initially before decaying to zero. Due to fluctuations, isolated regions of susceptible individuals evolve, and unlike in the standard mean-field SIR model, one observes a finite stationary distribution of the S type even for large population size. The exact results for the ensemble-averaged population size are compared with simulations for single realizations of the process and also with standard mean-field theory, which is expected to be valid on large fully connected graphs.

There is a related paper on arXiv by these authors, that seems to be the same paper, with a slightly different title (this one is free to download, for those without subscriptions to APS).

  • Gunter M. Schütz, Marian Brandau, Steffen Trimper, Exact solution of a stochastic SIR model. (2011). 0806.4440

This is the work we will review today. Before we do this, we will present a quick overview of related ideas we’ve considered to date.

The SIR model

The traditional SIR model consists of a fixed number of individuals TT split into three classes. (note, the authors use NN but we use this elsewhere in this series to represent the number operator).

  • SIR (classes: S - susceptible/ I - infected/ R - recovered), SIR epidemic model.

  • In Network Theory 3 John Baez explained the SIR model. Network Theory Part 3

The SIR model is said to have been invented by Kermack and McKendrick in 1927:

  • W. O. Kermack and A. G. McKendrick, A Contribution to the mathematical theory of epidemics, Proc. Roy. Soc. Lond. A 115 (1927), 700-721.

At any particular time, the average number of individuals in each class is given by

  • S, I, R with S+I+R=TS+I+R=T

The time evolution of the three classes is governed by the following set of coupled non-linear differential equations.

dSdt=βSI \frac{d S}{dt} = -\beta S I
dIdt=βSIαI \frac{d I}{dt} = \beta SI -\alpha I
dRdt=+αI \frac{d R}{dt} = + \alpha I
  • The first of these equations describes the reduction of SS via infection.

  • The second shows II increasing via infection of the susceptibles and decreasing by spontaneous recovery

  • The third shows the increase of RR due to the spontaneous recovery.

These are called the “rate-equations”. They account for how an average number of things changes with time. The other option is to consider the master equation, which considers the complete history of a stochastic process. Most people suspect that you arrive at the rate-equation through a large number approximation of the master equation. However, very recently this view was smashed with a counter intuitive result by John actually. This correspondence is indeed replaced with a more subtle result which is physically relevant in several regimes.

  • John Baez, Quantum techniques for reaction networks. [provide good link! since the one on the networks page is not working!]

TODO: explain the above better!

Ok, so back to these rate equations. The system can be solved numerically, yielding the time dependence of each class of individuals. Analytical solutions have several advantages, and are not known for this model in general. Special cases admitting exact solutions are interesting, and is what we will consider here.

The Petri net for the model is given as follows.

SIR model

From this Petri net, we can arrive at the stochastic-Hamiltonian, using our all-too-familiar prescription.

h=β(I +I +S I S +S I +I )+α(R +I I +I ) h = \beta(I^+I^+S^-I^- - S^+S^-I^+I^-) + \alpha (R^+I^- - I^+I^-)

Here we consider an initial state of a system. We will represent this as a formal power series in indeterminates i,s,ri, s, r as

ψ= n,m,k=0γ nmki ns mr k \psi = \sum_{n,m,k=0} \gamma_{nmk}i^n s^m r^k

We recall that the creation (I +,S +,R +)(I^+, S^+, R^+) and destruction operators (I ,S )(I^-, S^-) act as

I (i ns mr k):= i(i ns mr k)=(i n1s mr k)n I^- (i^n s^m r^k) := \partial_i (i^n s^m r^k) = (i^{n-1} s^m r^k)n

and as

I +(i ns mr k)=(i n+1s mr k) I^+ (i^n s^m r^k) = (i^{n+1} s^m r^k)

We’re reviewed all of these properties in much more detail throughout this series.

Looking at these equations, it’s useful to consider exactly what we mean by “spacial degrees of freedom”. Here, a term such as i ns mr ki^n s^m r^k has no reference to spacial location for an individual.

In real life, we know very well that epidemics have a strong spacial dependence in most cases. (Examples where this is not true would include a scenario where you place all of your rabbits inside of the same cage and they all interact approximately equally).

To motivate how we will consider position in our model let’s consider now just two rabbits. These two rabbits are initially separate. One of them is infected and the other is susceptible. We could write their state as

ψ 0=si \psi_0 = si

When the rabbits are allowed to interact, they are exposed to a situation omitting the following transition.

silongmapstoi 2 si \longmapsto i^2

Of course, the healthy rabbit could get lucky meaning that it would not become infected. Either of these two scenarios is possible. Assuming the healthy rabbit is infected after some time with probability pp the final state of the system would then be

ψ 1=pi 2+(1p)si \psi_1 = p i^2 + (1-p)si
  • But how do we determine pp?

In quantum mechanics, the interaction strength between two systems is typically determined by some sort of first principles law (e.g. dipole–dipole interaction etc.). Although we can use quantum techniques to model stochastic systems, such as epidemic modelling, it’s certainly not as clear how to determine what the interaction strength should be. For this one must turn to empirical data, leaving first principles account to be inaccurate. This is indeed much more difficult and must often by handled on a case-by-case basis when modelling epidemics.

Chemical reaction network analysis

In this subsection, we use

The chemical reaction network is given from the transition equations.

S+I2I S+I \rightarrow 2I
IR I \rightarrow R

The set of species is

{S,I,R} \{S, I, R\}

We will now calculate the deficiency of the resulting network.

  • The complexes. We let v 1:=(1,1,0) v_1 := (1,1,0)^\top, v 2:=(0,2,0) v_2 := (0,2,0)^\top, v 3:=(0,1,0) v_3 := (0,1,0)^\top and v 4:=(0,0,1) v_4 := (0,0,1)^\top denote the four complexes.

  • The stoichiometric subspace of the whole network is given by span{v 2v 1,v 4v 3}\text{span}\{v_2 - v_1, v_4 - v_3\} which has dimension 2.

  • The total number of linkage classes is two. One for each connected component. These are (i) v 1v_1 and v 2v_2 connecting, (ii) v 3v_3 and v 4v_4 connecting.

  • The deficiency of the network is the number of complexes (C=4C=4) subtract the dimension of the stoichiometric subspace (dim =2= 2), subtract the number of linkage classes (l=2l=2), which is zero.

  • (Remark) The SIR model is not weakly reversible.

  • (Result) The SIR model has deficiency zero, but is not weakly reversible, so cannot admit a steady state. This follows from the deficiency zero theorem.

We apply a result from the following paper.

  • Feinberg, M. (1987). Chemical reaction network structure and the stability of complex isothermal reactors: I. The deficiency zero and deficiency one theorems. Chemical Engineering Science 42 (10): 2229–2268. doi:10.1016/0009-2509(87)80099-4.

From Feinberg page 2252.

  • “r networks of deficiency zero that are not weakly reversible we know from the Deficiency Zero Theorem, part (i), that there can be no positive steady states at all; and so, for such networks, it is pointless to ask whether there can be multiple steady states within a positive stoichiometric compatibility class. Networks of nonzero deficiency present a very different situation. For them, weak reversibility is not a precondition for the existence of a positive steady state. Thus, for networks of nonzero deficiency, the issue of multiple positive steady states becomes a serious one even in the absence of weak reversibility.’‘ [Feinberg]

Dual neighbour model and its quantum mechanical representation

In contract to the fully general SIR model, in the paper [WMM11] the authors consider a stochastic one-dimensional adaptation of the SIR model. They show that this model admits exact solutions.

  • The model they consider is a 1D spin-1 model with periodic boundaries, with a three-body (pseudo) Hamiltonian. The Hamiltonian only acts on three sites in a row, e.g. operators O i1,i,i+1O_{i-1, i, i+1}.

In other words, each node is in contact with only two other nodes, and the rate of infection of susceptible individuals increases when neighbours are infected. Not all transitions are allowed. The authors define the interaction of the considered model to be

β:ISIIII \beta: ISI \rightarrow III
λ:ISSIIS \lambda: ISS \rightarrow IIS
λ:SSISII \lambda: SSI \rightarrow SII
α:IR \alpha: I \rightarrow R
  • Translation. Terms such as ISI is a state defined on three sites next to each other in the line. III means three infected members, right next to each other, etc.

Mean field deterministic model

The paper briefly describes mean field equations. We won’t consider this here at this stage.

Master equation

They let CC represent a state (S,I,R)(S,I,R) of each of NN individuals. Here NN is adjacent to 11 from periodic boundaries. The authors then state.

  • The master equation expresses the rate of change of the probability P(C,t)P(C,t) of finding the system in configuration CC at time tt as the rate of transfer of probability into CC from other configurations less the rate at which CC passes probability into the others.
dP(C,t)dt= CC{r[CC]P(C,t)r[CC]P(C,t)} \frac{d P(C,t)}{dt} = \sum_{C'\neq C} \{r[C'\rightarrow C]P(C',t) - r[C\rightarrow C']P(C,t)\}
  • The transition rate r[CC]r[C\rightarrow C'] is the probability per unit time that configuration CC changes into a different configuration CC'.

In quantum mechanics notation, a general state is then expressed as

|P(t)= CP(C,t)|C |P(t)\rangle = \sum_C P(C,t) |C\rangle

where P(C,t)P(C,t) is the probability that the system would be found in configuration CC at time tt.

They then cast the master equation into a quantum form

ddt|P(t)=H|P(t) \frac{d}{dt} |P(t)\rangle = - H |P(t)\rangle

and say that the pseudo-Hamiltonian HH has matrix elements.

C,H(C)=r(CC),CC \langle C', H(C)\rangle = -r(C\rightarrow C'), C'\neq C
C,H(C)= CCr(CC) \langle C, H(C)\rangle = \sum_{C'\neq C} r(C\rightarrow C')

with the known formal solution

|P(t)=e Ht|P(0) |P(t)\rangle = e^{-H t}|P(0)\rangle

In this formalism the expectation value of M(C)M(C) is given as

M= C,sP(C,t)M(C)=s|M^|P(t) \langle M \rangle = \sum_{C', s} P(C,t)M(C) = \langle s|\hat M |P(t)\rangle
|s= C|C |s\rangle = \sum_C |C\rangle
M^= CM(C)|CC| \hat M = \sum_C M(C)|C\rangle \langle C|

We have the standard equation of motion for the time dependence of the expected value of MM as

Mt=s|[M^,H]|P(t) \frac{\partial \langle M \rangle }{\partial t} = \langle s|[\hat M, H] |P(t)\rangle

Creation and annihilation formalism

  • Type S particles: a i a_i^\dagger (creates); a ia_i (annihilates)

  • Type I particles: b i b_i^\dagger (creates); b ib_i (annihilates)

  • A site with neither an S nor I is assumed R (recovered)

  • Operators at different sites commute.

  • Operators of the same species at the same site anti-commute.

  • The number operators A i:=a i a iA_i := a_i^\dagger a_i, B i:=b i b iB_i := b_i^\dagger b_i take only values zero or one.

In this notation, the Hamiltonian is written as (Equation 8)

H=β i[b i a iA i(1B i)][γB i1B i+1+A i1B i+1+B i1B i+1]+δ i[b iB i] - H = \beta \sum_i [b_i^\dagger a_i - A_i(1-B_i)][\gamma B_{i-1}B_{i+1}+A_{i-1}B_{i+1}+B_{i-1}B_{i+1}] + \delta \sum_i[b_i - B_i]


γ:=βλ \gamma := \frac{\beta}{\lambda}
δ:=αλ \delta := \frac{\alpha}{\lambda}

Exact solution: derivation and analysis

Cluster functions

They use a method of nn-point cluster functions.

K r(n):=A rA r+1A r+n1B r+n K_r(n) := \langle A_r A_{r+1} \cdots A_{r+n-1} B_{r+n} \rangle
G r(n):=B r1A rA r+n1B r+n G_r(n) := \langle B_{r-1} A_r \cdots A_{r+n-1} B_{r+n}\rangle
  • (Symmetry) The Hamiltonian is translationally invariant and we consider periodic boundaries so we are dealing with the symmetry group C NC_N. If the initial state is translationally invariant, so are all future configurations.

  • Translational invariance make the quantities K r(n)K_r(n) and G r(n)G_r(n) independent of their initial starting point rr.

  • The time dependence of the cluster functions are found from the equation of motion of the operators, which involves calculating the commutation relations between the Hamiltonian and cluster operator products.


dK(1)dt=s|[A rBr+1,H]|P(t)=δK(1)γG(1) \frac{d K(1)}{d t} = \langle s| [A_rB_{r+1}, H]|P(t)\rangle = -\delta K(1) - \gamma G(1)
dG(1)dt=s|[B r1A rBr+1,H]|P(t)=(γ+2δG(1)+2G(2) \frac{d G(1)}{d t} = \langle s| [B_{r-1}A_rB_{r+1}, H]|P(t)\rangle = -(\gamma+2\delta G(1) + 2G(2)

nn greater than 1

dK(1)dt=s|[A rA r+n1B r+n,H]|P(t)=K(n+1)(1+δ)K(n)G(n) \frac{d K(1)}{d t} = \langle s| [A_r\cdots A_{r+n-1}B_{r+n}, H]|P(t)\rangle = K(n+1)-(1+\delta)K(n)-G(n)
dG(1)dt=s|[B r1A rA r+n1B r+n,H]|P(t)=2G(n+1)2(1+δ)G(n) \frac{d G(1)}{d t} = \langle s| [B_{r-1}A_r \cdots A_{r+n-1}B_{r+n}, H]|P(t)\rangle = 2G(n+1)-2(1+\delta)G(n)

Particle densities of I and S

The equations of motion for the particle density of susceptibles and infecteds are found in a similar fashion.

dA rdt=s|[A r,H]|P(t)=2K(2)γG(1) \frac{d A_r}{d t} = \langle s| [A_r, H]|P(t)\rangle = -2K(2) - \gamma G(1)
dB rdt=s|[B r,H]|P(t)=2K(2)+γG(1)δB r \frac{d B_r}{d t} = \langle s| [B_r, H]|P(t)\rangle = 2K(2)+ \gamma G(1)-\delta \langle B_r \rangle

Solution for cluster functions

Particle densities

Analysis of solutions



  • [WMM11] H. Thomas Williams, Irina Mazilu and Dan Mazilu, Stochastic epidemic-type model with enhanced connectivity: exact solution. (2011). arXiv:1108.5135v1
  • [NTP] Network Theory Part 1-7, 2, 3, 4, 5, 6, 7

Homepages of researchers

The paper was published as

  • H. T. Williams, I. Mazilu, D. A. Mazilu, “Stochastic epidemic-type model with enhanced connectivity: exact solution”, Journal of Statistical Mechanics: Theory and Experiment, 1742-5468, P01017 (2012).


H T Williams et al J. Stat. Mech. (2012) P01017 doi:10.1088/1742-5468/2012/01/P01017

Scrap for later

It is likely best at this stage to define a basis, for the spin-1 particles.

|S,|I,|R |S\rangle, |I\rangle, |R\rangle

We then have

S,S=I,I=R,R=1 \langle S,S\rangle =\langle I,I\rangle =\langle R,R\rangle = 1


S,R=I,R=I,S=0 \langle S,R\rangle =\langle I,R\rangle =\langle I,S\rangle = 0

category: experiments