Exact equations for SIR epidemics on tree graphs

Exact equations for SIR epidemics on tree graphs

K.J. Sharkey, I.Z. Kiss, R.R. Wilkinson, P.L. Simon

Abstract

We consider Markovian susceptible-infectious-removed (SIR) dynamics on time-invariant weighted contact networks where the infection and removal processes are Poisson and where network links may be directed or undirected. We prove that a particular pair-based moment closure representation generates the expected infectious time series for networks with no cycles in the underlying graph. Moreover, this “deterministic” representation of the expected behaviour of a complex heterogeneous and finite Markovian system is straightforward to evaluate numerically.

1 Introduction

1.1 Background

The majority of epidemic models fall either into the category of stochastic models (Bailey 1975; Bartlett 1956) or into the category of deterministic differential equation-based models (Anderson and May 1991; Kermack and McKendrick 1927). These two strands developed largely independently for much of the twentieth century. Thus, an interesting question arises as to the precise mathematical connection between stochastic and deterministic models. Frequently, deterministic descriptions apply to large populations where the stochastic effects can be treated as negligible. For small populations we shall assume that it is the average or expected behaviour of the epidemic that we are hoping to replicate with “deterministic” descriptions. This average behaviour is a system characteristic that is fully specified by the system and its initial conditions.

The first epidemic models were based on the assumption that populations are evenly mixed, with each individual equally likely to interact with any other individual at any time (Heathcote 2000). A classic example of this type of model is the Susceptible-Infectious-Removed (SIR) compartmental model whereby individuals are classified according to being in one of these three states. It has been shown that for this type of mean-field model, the average of many stochastic simulations (the expected outcome of the stochastic model) converges to the solution of the “equivalent” mean-field deterministic model in the limit of an infinite population size and subject to strict conditions regarding the initialisation of the epidemic (Kurtz 1970, 1971; Simon and Kiss 2011).

More recently, a higher degree of realism has been introduced by considering stochastic models on contact networks where individuals are only able to contact a limited subset of the population. This enables significant heterogeneity to be incorporated, treating individuals as distinct entities with fixed connectivity to pre-allocated neighbours. While stochastic models are readily extended to incorporate such systems, deterministic descriptions have been more problematic. Several methodologies have been developed including pair-approximation models (Keeling 1999; Keeling and Eames 2005; Rand 1999), degree-based models (Pastor-Satorras and Vespignani 2001), and models based on the probability generating function (PGF) formalism which are applicable to configuration networks (Volz 2008) as well as the related edge-based compartmental modelling (Miller et al. 2012; Miller and Volz 2012). It has been observed (House and Keeling 2011) that these models are, at some level, equivalent and are all derived from similar principles of independence. Although comparison with simulation of stochastic models can sometimes be good, the basic link remains obscure.

Typically there are two idealised scenarios in which exact correspondence between stochastic models and solvable deterministic descriptions has been shown. Firstly, correspondence has been shown to sometimes occur in the limit of infinite populations for particular idealised graphs (Ball and Neal 2008; Decreusefond et al. 2012) which cannot be exactly realised in practice. It can also occur with some very simplified systems whose symmetry properties can be exploited to achieve reductions in the stochastic description (Keeling and Ross 2008; Simon et al. 2011).

Here we consider a recently introduced class of model, related to the pair-approximation models, which give an exact correspondence between a deterministic description and the stochastic model for SIR epidemics on finite, time-invariant networks. Pair-approximation models were introduced into network-based epidemic and ecological theory in the 1990s to describe large populations of interacting individuals (Matsuda et al. 1992; Sato et al. 1994; Harada and Iwasa 1994; Rand 1999; Keeling 1999). They are an example of a hierarchy of equations which are truncated at the second order by an approximation (truncation at the first order corresponds to mean-field). This type of hierarchy was first considered in statistical physics and is sometimes known as the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy (Kirkwood 1946, 1947; Born and Green 1946). Recently, related models have been considered at the level of individuals, variously called subsystem equations, moment dynamics equations, pair-based equations (Sharkey 2008, 2011; Baker and Simpson 2010; Markham et al. 2013). This method generates a solvable class of models which can encompass a significant amount of heterogeneity and enables a fundamental link with finite stochastic models (Sharkey 2008, 2011).

We consider a pair-based representation of Markovian SIR dynamics. We show that by considering subsystems at the level of pairs, a closure can be found that determines the expected infectious time series exactly for arbitrary network structures where the underlying graph is a tree and, in some special circumstances, for particular networks with cycles. We note that the recent, related message passing formulation of epidemics on contact networks developed by Karrer and Newman (2010) also enables an exact description of epidemic dynamics on finite tree graphs.

1.2 Statement of the main result

We consider an SIR compartmental model composed of individuals whose states are described at any given point in time by vectors and with respective components and , such that if individual is infectious ( otherwise) and if individual is susceptible ( otherwise). Transmission and recovery occur by Poisson processes with rate parameters and , respectively where is a “transmission” matrix with (time-independent) elements denoting the rate parameter for an infectious node infecting a susceptible node ( for all ) and where denotes the rate parameter for an infectious individual to recover, enabling individual-specific removal rates.

As shown by Sharkey (2011), for any transmission matrix and any nodes , the following differential equations are provably exact (consistent with the stochastic model):

(1)

where and denote the time-dependent probabilities (or equivalently the expected values of the indicator functions) for individual to be susceptible and infectious, respectively, and expressions of the form denote the time-dependent probability that individual is in state and individual is in state with a similar interpretation of terms of the form . Here and throughout, we adopt the dot notation to denote time derivatives. It follows that the expected population-level susceptible and infectious time series are given by and respectively.

Note that it is a short step (see Sharkey 2008) from (1) to the familiar population-level pair equations (Keeling 1999; Keeling and Eames 2005; Rand 1999), also proved independently by Taylor et al. (2012) for the susceptible-infectious-susceptible variant.

This system can be completed by formulating differential equations for the triples, quadruples, and so forth until we reach the full system size. This yields a self-contained system of differential equations that exactly determines the probabilities of each quantity given initial conditions. However, cascading these equations up to the full system size will usually result in a system that is impractical to solve due to its sheer size. This is why this system is typically closed at some level by introducing a functional relation approximating higher-order probabilities in terms of lower-order ones. One of the most frequently used closure relations can be written as

(2)

for the current context. Applying this closure relation to our system at the level of pairs we arrive at the following system:

where we use for susceptible and for infectious to emphasise that these are approximating differential equations based on the closure. When in the denominator is zero, we assume that the approximation takes the value zero.

In general, we consider networks (graphs) with directed and undirected edges. In what follows, we use the terminology “tree graph” to include graphs with directed edges where the underlying (equivalent undirected) graph is a tree. Our main aim is to show that when matrix represents a tree and the system is initiated in a pure system state (that is, one of the possible configurations has probability 1 at time ), then the system can be closed at the level of pairs such that the closure holds exactly. Specifically, solving the closed system above, we obtain the same values for all marginal and pairwise-joint probabilities present in the unclosed system: , , with similar equalities holding for the pairs.

In fact, we will prove the following theorem.

Theorem 1.1.

Let us assume the following:

  • The graph (transmission network) is a tree (the underlying graph has no cycles).

  • The initial condition is a pure state, i.e. the system is initially in one of its possible configurations with probability 1.

Then the following relations hold:

for all and for all with links towards and all with links towards : ;

for all and for all and with links towards : .

Remark.

This theorem also holds for mixed (probabilistic) initial system states provided that the initial probabilities of the states of individuals in the system are uncorrelated. However, in general, mixed initial states cannot be represented exactly.

The theorem will be formulated in a more general context stating that even higher-order closure relations are also exact.

Figure 1

Figure 1: a) An undirected tree indicating two nodes which we infect to initiate epidemics, with all other nodes initially susceptible. b) The mean (dots) of 100,000 stochastic simulations on the network with transmission rate across each link and removal rate for each node, with error bars denoting the 5th and 95th percentiles plotted together with the solution of (LABEL:0.2) (solid line) using the Matlab code published with Sharkey (2011)

shows the numerical solution of (LABEL:0.2) for a small network of 9 nodes where it is clear that it is accurate to within the precision visible on the graph. Matlab code for solving the system of equations (LABEL:0.2) is provided in Sharkey (2011). This code also works on networks which are not trees but is no longer exact in these cases. Cycles in the underlying graph of order three utilise the alternative closure which is believed to gain increased accuracy in most circumstances, but these do not occur in the tree graphs considered in the present work.

The structure of the paper is as follows. Section 2 introduces some notation which is needed to prove the result. This also contains an important theorem (Theorem 2.1) which specifies equations describing the probabilities of the states of arbitrary subsystems (the proof of this result is given in Appendix A). The relevant state space for our domain of a tree graph is then developed. Section 3 proves the main result, initially focusing on some special cases to help motivate and facilitate understanding of the main ideas and steps of the general proof in Section 3.5. The main ingredient for the general proof is Lemma 3.1 which is proved via Theorem 2.1. Theorem 3.1 then follows easily by induction from Lemma 3.1. The theorem as stated above is a simple corollary of Theorem 3.1. In Section 4 we discuss an application of the pair-based model to some special cases of graphs with cycles where it is also exact.

2 Formulating the full system

In this section we introduce a new notation which will assist in formulating the set of differential equations for the full system. In (1) we formulated the differential equations up to the level of pairs and we said that this could be continued up to the full system level. This will be done formally here. In order to make the method clearer, using our existing notation let us first evaluate the full set of equations for the undirected line graph with three nodes which we refer to as the open triple, depicted in Figure 2.

Figure 2: Open triple graph

Here we shall assume that the transmission rate parameter is across both links and that the removal rate is for all three nodes. Firstly we write all of the single node equations for this network. From (1):

(4)

and

(5)

We also need to specify the following equations for pairs:

(6)

Finally, at the triple level we have from the master equation (since the system has only three nodes):

(7)

In order to formulate the full system for an arbitrary graph, we introduce notation for the subsystem states.

2.1 Notation for system and subsystem states

In general, our stochastic system (which we denote by ) comprises of individuals, each of which may be in any of the , or states at any given time. In total, this corresponds to possible states. Denoting these system states by , , the probabilities for each state are given by the master equation (or Kolmogorov equations):

(8)

where denotes a constant matrix of Poisson rate parameters. The master equation completely describes our stochastic system using a set of ordinary differential equations. Our overall objective is to show that (LABEL:0.2) is implied by the master equation when represents a tree graph.

It is useful for us to define a general subsystem comprising of nodes in indexed by vector of length : , where we can assume . We assume that the network connections of the nodes of are a subset of the connections of .

Let denote the state of subsystem where and is a sequence of and symbols of length such that the state of node is . In terms of the notation of the previous section for subsystems of single nodes and pairs of nodes, we have , , , etc. We shall use these two notations interchangeably. We shall also sometimes treat indexing vectors such as as sets such that means that the node is in the subsystem .

In general, although we can specify the states of each node with this type of notation, an important ambiguity remains because information about the network structure is not included. To remove this ambiguity, the notation should normally be used in the context of a sketch of the relevant network structure or where the network structure is clear from the context of its use (as in (1)).

Let us now show how the differential equations of the different subsystem states can be formulated in general.

2.2 Differential equations for subsystems

Here we obtain differential equations describing the rate of change of the state of any subsystem. First we make some definitions.

Definition 2.1.

A neighbour of node is a node with a network link directed towards .

Definition 2.2.

denotes the set of neighbours of node . That is: .

Definition 2.3.

For the subsystem state , if node is infectious then:

Otherwise, .

Remark.

This operator changes the state of node in subsystem to if it is infectious. If node is susceptible or removed then it leaves the state unchanged.

Definition 2.4.

For the subsystem state of nodes, a subsystem of nodes can be generated as follows: Take and take a neighbour of outside of the subsystem with a network link towards , i.e. let , . If , then the generated subsystem state of nodes is given by the generating rule:

i.e. the subsystem is extended by an infected at node which is connected towards . If , then the generated subsystem state is given by:

i.e. the subsystem is extended by an infected at node which is connected towards and the state of node is changed from to .

To complete the definition, if then the operator leaves the subsystem unchanged. We also assume that for any state where there is no link from node to node in the transmission matrix , then the subsystem is also left unchanged.

Remark.

The generated order subsystem is obtained by replacing a susceptible or infectious node in the original subsystem by an arc such that the node of the arc is put in the place of the node and where the node of the arc is external to the subsystem.

Definition 2.5.

For the subsystem , if node is removed then:

Otherwise, .

Definition 2.6.

For any subsystem of nodes in state we define where and to have value 1 if and to have value zero otherwise:

Theorem 2.1.

The rate of change of the probability of a subsystem state is:

(9)

The proof of this theorem is a rather long diversion and can be found in Appendix A.

As an example of applying the theorem, we can use it to obtain the set of subsystem equations (1) by considering each equation in turn:

  • If the subsystem is a single susceptible individual , then so can only take the value where and , reducing (9) to:

    The first term on the second line of (9) is zero because , and the other terms are zero because and .

  • For an infectious individual we obtain:

    where the first term on the second line of (9) is zero because and the last term is zero because .

  • If the subsystem is the pair then the sum over is over and and , , , so:

    where the first line corresponds to and the second to .

  • If the subsystem is the pair then the sum is over and where , , and so:

    where both terms come from the first line of (9).

We have therefore obtained (1) in a slightly different notation (recall that ).

2.3 The state space for a tree graph

Here we build up a state space which is sufficient to describe a tree graph. We first make some definitions.

Definition 2.7.

An -motif is a subsystem of comprising of nodes and of network links such that it forms a weakly connected network.

Definition 2.8.

An -state is the state of an -motif.

The state space that we need to consider is built up inductively from the states of single nodes by considering the infection process. Starting with the infected states of the single nodes , , (9) shows that they depend on the 2-states , as described by the generating rule (Definition 2.4).

The differential equations for in turn contain the 3-states , and , . The differential equations for the 3-states contain 4-states and typically, the differential equations for -states contain -states for . This state generation process can continue until we reach -states which can only depend on other -states.

Note that this process always forms subsystems which are motifs and that the motif states can never include removed nodes.

Definition 2.9.

An out-neighbour of node is a node with a network link from towards it.

Proposition 2.1.

For a tree graph, if the out-neighbours of the nodes are all in an -motif, then this is true for all -motifs generated from this -motif.

Proof.

This follows easily from the definition of the generating rule (Definition 2.4). ∎

Definition 2.10.

Consider a tree graph and take the 1-motifs with nodes: . The “basic state space” is formed by these 1-states together with the set of motif states that can be iteratively generated from them using the generating rule (Definition 2.4).

Remark.

Due to the method of its construction, the state space gives a self-contained system of differential equations, i.e. the time derivatives of the probabilities of each motif state can be expressed in terms of the probabilities of other motif states in the state space. An example in the case of the open triple is given by the motifs in (4), (6) and (7).

Definition 2.11.

Consider a tree graph and the 1-states: and the 2-states with , i.e. . The “extended state space” comprises of these motif states together with the set of motif states that can be generated from them by repeated iteration of the generating rule.

Remark.

The extended state space is required to form the relevant closure relations. Due to the method of its construction, it is also self-contained.

Lemma 2.1.

Let . Then the out-neighbour of an node is an node in .

Proof.

Follows from Proposition 2.1. ∎

Lemma 2.2.

For a tree graph, the equation for the time derivative of the probability of an -state is given by:

(10)
Proof.

For these states we have . Additionally, when , the first term on the second line of (9) never arises because implies that an is connected to an node in r-state which contradicts Lemma 2.1. Therefore (9) reduces to (10). ∎

Let us now formulate the exact closure relations and prove our main result.

3 Closure relation and proof of the main result

The exactness of (LABEL:0.2) is straightforward to see provided that outbreaks of epidemics are always initiated with a single infected individual. We prove this first before considering the general case.

3.1 Proof for single initial infected

When infection is initiated on a tree graph at a single individual, infection must always proceed in linear chains. Consequently there is no possibility of the state illustrated in Figure 3

Figure 3: Shown is a state which cannot arise on a tree graph where there is only one initially infectious node

arising because an infection initiated at either or must pass through to get to the other node. Furthermore,

but since and consequently , we have:

reducing (1) to the following closed system:

Similar arguments show that this can be written in the form of (LABEL:0.2).

More generally, this argument also applies to any tree graph where there is at most one network path by which any susceptible individual in the network can become infectious from the initial configuration of infected individuals.

Before discussing the general proof for any tree graph with multiple initially infected individuals, we consider two very simple example networks which will serve to motivate and illustrate the method of proof.

3.2 Proof for an open triple

Here we consider the case for the open triple depicted in Figure 2. The equations for the probabilities of the basic state space are given in (4), (6) and (7). To form the relevant closure relations, we require the equations for the extended state space formed by the equations for together with (5),

(11)

Our objective is to close the system at the level of pairs using the closure relation (2), eliminating the need for differential equations describing triples (7), and show that the system remains exact. We note that the exactness of (LABEL:0.2) can be proved in this case along the lines of the previous argument by considering each possible initial condition separately; however, the approach discussed here will be more useful for understanding the general case.

We need to consider closures for the triples , and . Let us consider the closure:

This is exact if where

and . Taking the derivative of with respect to time gives

Substituting the relevant derivatives in from (5)-(7) and cancelling terms reduces this to

so:

Now it is easily verified that provided the system is initiated in a specific system state then . Consequently for all and the closure is exact.

By symmetry, it will suffice to consider one of the remaining two triples in (6). We wish to show that where

Here it is necessary to also use (11) for pairs of type SS in the extended state space. This closure is not established immediately, but there is a two-step process to establishing that which the reader can verify by analogy with the example of the star graph in the next section.

3.3 Proof for a star graph

We now consider the case of the undirected star graph with shown in Figure 4,

Figure 4: Star graph with nodes

where again we assume that the strength is the same across each network link and is denoted by and the removal rate for each node is . Writing down the equations of the extended state space, there are two types of closure which need to be proved: one for the triples and one for the triples (see (1)). The graph has three triples , but it is sufficient to prove exactness for one of them. Hence we want to prove the following two relations:

(12)

For brevity, we adopt the alternative notation:

We introduce:

By differentiating this, substituting in from the process equations and grouping terms, we obtain

(13)

where:

Differentiating we get

(14)

where:

The derivatives of and can be obtained similarly.

Differentiating we get

(15)

where:

The derivative for can also be obtained. Finally, differentiating and we obtain:

(16)

To conclude the proof of the exactness of the closure, we first assume that the initial state is not mixed; that is, one of the possible configurations has probability 1 at . Then it is easy to see that for all (see Lemma 3.2 in Section 3.5 for a proof in a more general context). Hence the differential equations for and show that and for all . The differential equation for then implies that (and similarly for ). This implies that (and similarly ). The differential equation for shows that which is what we wanted to show. The other triple closure in (12) can be proved similarly.

Remark.

In fact, we have proved several closure relations .

The closure relations each consist of two pairs which are visualised in Figure 5. For reference, we refer to these as the left pair and the right pair referring to their position in this figure.

Figure 5: Each box illustrates the relevant node states for the four parts of the closure relation in the equation above it. The node states on the left and the right correspond to the two terms in the closure relation. The node numbers correspond to the same positions as in Figure 4

Looking at these closure relations, we can form two observations:

  1. For a given node , the number of times it appears as is the same in the left and the right pair, and similarly with the number of times it appears as . For example, with , node 1 (the left node) has one and one for both pairs and node 4 (central node) has two ’s in both pairs. For , the number of ’s at nodes 1,2,3 and 4 is in both pairs and the number of ’s is given by in both pairs.

  2. Any SI pairing on the left appears exactly the same number of times on the right. For example in , appears once on the left and once on the right and appears twice on the left and twice on the right. Observing that only pairs and pairs appear in the closure relations, a consequence is that pairs also have this property.

These observations will be of key importance for developing the general proof in the following two sections.

3.4 General closure relations

In general, to show that the closure relationship (2) is exact for the tree graph, we need to show that where

and . Our proof of this is via induction using a sequence of closures analogous to the proof in the case of the star graph in Section 3.3.

We shall consider many closure relations. In general we specify that they are composed of two pairs of motif states and and that the closure is exact if where

We formalise the observations we made about the closure relations for the star graph at the end of Section 3.3 by defining what we term “compatible pairs”.

Definition 3.1.

For all .

For all