Complete hierarchies of SIR models on arbitrary networks with exact and approximate moment closure

0.65 MB·English

by  Kieran J Sharkey

Complete hierarchies of SIR models on arbitrary
5 networks with exact and approximate moment
1
0
closure
2

r
p Kieran J. Sharkeya∗and Robert R. Wilkinsona
A

4
2
(a)DepartmentofMathematicalSciences,UniversityofLiverpool,PeachStreet,

Liverpool, L69 7ZL
]
E
P
Abstract
.
o
i
b We first generalise ideas discussed by Kiss et al. (2015) to prove a theorem for
- generating exact closures (here expressing joint probabilities in terms of their
q
[ constituentmarginalprobabilities)forsusceptible-infectious-removed(SIR)dy-

namics on arbitrary graphs (networks). For Poisson transmission and removal
2
processes, this enables us to obtain a systematic reduction in the number of
v
3 differential equations needed for an exact ‘moment closure’ representation of
5 the underlying stochastic model. We define ‘transmission blocks’ as a possible
3 extensionoftheblockconceptingraphtheoryandshowthattheorderatwhich
6
the exact moment closure representation is curtailed is the size of the largest
0
transmission block. More generally, approximate closures of the hierarchy of
.
1 momentequationsforthese dynamicsaretypicallydefinedforthe firstandsec-
0
ond order yielding mean-field and pairwise models respectively. It is frequently
5
impliedthat,inprinciple,closedmodelscanbewrittendownatarbitraryorder
1
: if only we had the time and patience to do this. However, for epidemic dy-
v
namicsonnetworks,thesehigher-ordermodelshavenotbeendefinedexplicitly.
i
X Here we unambiguously define hierarchies of approximate closed models that
r canutilise subsystemstates ofany order,andshow how well-knownmodels are
a
special cases of these hierarchies.
1 Introduction
A primary method for incorporating spatial structure and other contact struc-
turesintoepidemic modelsis touse anetworkofcontacts[1]. While simulation
of stochastic models is straightforward on these networks, obtaining differen-
tial equation descriptions of the relevant time series is more complex. Here we
∗Correspondingauthor(Email: [e-mail];Tel: [phone]).
1
considerthe constructionofa hierarchyofmomentequationswhich, instatisti-
calphysics,issometimesknownastheBogoliubov-Born-Green-Kirkwood-Yvon
(BBGKY) hierarchy after the names of its originators. The method was ap-
plied to population-level network-based epidemic and ecological models in the
1990swheretruncationofthehierarchywasmadeatsecondorderyieldingpair-
approximation models [2, 3, 4, 5]. Higher-order truncation of this hierarchy at
the level of triples has also been investigated [2, 6, 7]. With increasing compu-
tational resources it has also become numerically viable to consider these hier-
archiesintermsofindividuals,ratherthanpopulation-levelquantities[8,9,10].
Aparticularlyimportantfeatureofthe individual-levelrepresentationis that it
enables us to establish exactness for finite populations in certain circumstances
(see [11] and also [12] and [13] by different methods).
HerewegeneraliseideasdiscussedbyKissetal. [12]andalsonotedbyNew-
man [14], and apply them to arbitrarydirected networks. We also observe that
they apply to non-Markovian as well as Markovian SIR dynamics. Depending
onthenetwork,wefindthatforMarkoviandynamics,exactclosedmodelsexist
at all levels of the hierarchyof moment equations. The exact models and exact
closures considered in [11] and [12] then represent special cases.
While the majority of moment closure models do not go beyond closure at
the levelofpairs (secondorder),it is frequently statedthat, in principle, closed
models at any order can be constructed. However, such higher-order models
are rarely defined explicitly. Here, in the Markovian SIR epidemic context, we
shalldefinehierarchiesofclosedmodelsthatcanbeconstructedunambiguously
at all orders by a systematic truncation method. In fact, we shall define and
investigate several hierarchies of approximate models. All of these converge
to exact representations at truncation orders which depend on the underlying
network structure and all of them have either the pair-level model discussed
in[8]and[9]orthe variantofthis modeldiscussedin[11]asthe lowest(zeroth)
order level of truncation.
Thenextsectiondiscussestherelevantbackgroundconceptsuponwhichour
work builds. Section 3 introduces the exact closure theorem which defines the
conditionsunderwhichsimplificationstothehierarchyofequationscanbemade
for particular networks. Section 4 introduces approximate closures leading to
complete hierarchies of approximate models.
2 Background concepts
Apart from Theorem 3.1 which applies more generally, we shall consider a
Markovian class of SIR models on contact networks. In particular, we con-
sider a directed graph D =(V,A) consisting of N =|V| individuals/nodes and
asetAofarcs. Wealsolabeleachindividualaccordingtosomearbitraryorder-
ingsuchthatifi∈V theni∈{1,2,...,N}. Eachindividualcanbeinonlyone
ofthreestatesS,I,orRatanygiventime. Nodej ∈V,wheninfectious,makes
‘infectious contacts’ to node i∈V via a Poisson process of rate T ≥0, where
ij
T >0⇔(j,i)∈AandwhereweassumethatT =0foralli∈V. Ifnodeiis
ij ii
2
susceptible when it receives an infectious contact then it immediately becomes
infectious. Itwillthenremaininfectiousforanexponentiallydistributedperiod,
withparameterγ ,afterwhichitbecomesrecoveredwhichisanabsorbingstate
i
for the individual. We thus have a continuous-time Markov chain with a state
space of size 3N. Except where otherwise stated, we also assume initial con-
ditions such that the states of all nodes are initially statistically independent.
Thisassumptionencompassesallpure-stateinitialconditions,suchasaspecific
individual being infectious with all others susceptible, and it also incorporates
binomially distributed initial conditions. Uniform initial conditions can also be
exactly represented with additional computation [13].
Definition 2.1. S , I and R denote the indicator random variables for the
i i i
events that node i ∈ V is susceptible, infectious and removed respectively. De-
pending on the context, it will also be convenient to refer to S , I and R as the
i i i
corresponding events themselves.
The hierarchy comprises of a sequence of equations containing the first mo-
ments and mixed moments of the random variables S and I . Using angle
i i
bracketstodenoteexpectationvalues,itcanbeshown[9]thatthemasterequa-
tion (or Kolmogorovforward equations) implies the following rate equations:
N
˙
hS i = − T hS I i,
i ij i j
j=1
X
N
˙
hI i = T hS I i−γ hI i. (1)
i ij i j i i
j=1
X
where S I is a product of the indicator randomvariables which also specifies a
i j
state ofthe subsystemofordertwocomprisingofthe pair ofnodes i andj. For
this pair state we have:
N N
˙
hS I i = T hS S I i− T hS I I i
i j jk i j k ik i j k
k6=i k6=j
X X
−T hS I i−γ hS I i. (2)
ij i j j i j
More generally, for these models, the master equation allows us to write
down a rate equation for the probability of any subsystem state of size n in
termsofsubsystemstatesofsizen andsubsystemstatesofsizen+1. We state
this as a theorem below (Theorem 2.1).
Followingpriorwork[11],butwithanotationalsimplificationbroughtabout
by using the same index for all system and subsystem states, we define an
alternative notation to Definition 2.1 that is useful for keeping track of the
hierarchy of moment equations in this context.
Definition 2.2. We use the following notation to denote subsystem states.
• ψ is a subsystem comprising of the set of nodes W ⊂V.
W
3
• Let A be a mapping from the elements of W to {S,I,R}, and let A be the
i
image of node i∈W under A. Thus, A can also be interpreted as a pure
state for subsystem ψ , i.e. the state where, for all i∈W, individual i is
W
in state A .
i
• ψA denotestheindicatorrandomvariable fortheeventthatψ isinstate
W W
A. Thus the probability of the event that subsystem ψ is in state A is
W
P(ψ = A) = hψAi. As in Definition 2.1, it is also convenient to refer
W W
to ψA as the event that ψ is in state A.
W W
Remark. For the event where node i is in a susceptible state, we can draw the
followingcorrespondencebetween the notations: ψS =S , and similarlyfor the
i i
infectious and removed states.
Definition 2.3. Let k ∈ W ⊆ V and X ∈ {S,I,R} and let A be a state of
subsystem ψ . Then, hX(ψA) denotes the indicator random variable or event
W k W
ψA, but where the state of node k is changed to state X leaving the states of all
W
other nodes unchanged. Note that if A =X then hX(ψA)=ψA.
k k W W
Theorem 2.1. For any subsystem ψ , the probability that it is in state A is
W
governed by the rate equation:
hψ˙Ai = (A =S) − T (A =I)hψAi− T hψAI i
W 1 k " kn1 n W kn W n #
kX∈W nX∈W n∈XV\W
+ (A =I) T (A =I)hhS(ψA)i−γ hψAi
k kn n k W k W
1 " 1 #
k∈W n∈W
X X
+ (A =I) T hhS(ψA)I i
k kn k W n
1 " #
k∈W n∈V\W
X X
+ (A =R) γ hhI(ψA)i , (3)
k k k W
1 " #
k∈W
X
where here, and throughout this paper, the indicator (.) is equal to 1 if its
1
argument is true and is equal to zero otherwise.
This theorem is proved in [11]. Starting with subsystem states that are
only composed of susceptible or infectious individuals, repeated application of
equation 3 to each of these states as well as to any subsystem states that arise
on its right-hand side can never result in subsystem states with a removed
individual. This is due to the absence of hR in equation 3. Hence, for these
k
4
subsystem states, (A =R)=0 for all k ∈W so equation 3 becomes:
k
1
hψ˙Ai = (A =S) − T (A =I)hψAi− T hψAI i
W 1 k " kn1 n W kn W n #
kX∈W nX∈W n∈XV\W
+ (A =I) T (A =I)hhS(ψA)i−γ hψAi
k kn n k W k W
1 " 1 #
k∈W n∈W
X X
+ (A =I) T hhS(ψA)I i . (4)
k kn k W n
1 " #
kX∈W n∈XV\W
Equations 1 and 2 can now be seen to be special cases of this theorem.
By applying equation 4 to every individual in the network for states S and
I and then reapplying to every new subsystem state which emerges, we obtain
a closed set of differential equations for a set M of subsystem states. However,
|M|willgenerallybeverylargeformostsystems,preventingnumericalsolution.
To reduce the number of equations, we need to introduce a mechanism to
curtail the generation of new subsystem states. In the next section, we discuss
scenariosin which this can be done where the emerging system is still an exact
representation of the underlying stochastic process. Following this, we shall
consider hierarchies of approximate closed models.
3 Exact closed models
Here we prove a theorem pertaining to arbitrary SIR dynamics on arbitrary
networks. We then use this to derive a class of exact models for Markovian
SIR dynamics on arbitrary networks. We illustrate this with some examples,
and finally state a theorem specifying the maximum size of subsystem needed
to exactly represent the dynamics on any given network.
3.1 Exact closure theorem
For a given directed graph D =(V,A) with set V of nodes/individuals and set
A of arcs, we make the following definitions:
Definition 3.1. IN(X) is the set of individuals that can reach at least one
member of X ⊆V by following a directed path. Note that X ⊆IN(X).
Definition 3.2. Let X,Y,Z ⊂V be disjoint and non-empty. The set of nodes
Z is ‘dynamically partitioning’ with respect to X and Y if and only if we have
f (X,Y,Z)=1 where:
E
1 if IN(X)∩IN(Y)=∅ (in D−Z)
f (X,Y,Z)= (5)
E
(0 otherwise
and D−Z is the vertex-set deleted subgraph consisting of nodes V \Z. Here,
E is chosen to represent ‘exact’; this is appropriate since we shall now see that
f (X,Y,Z)=1 implies the existence of an exact closure relation.
E
5
Remark. If the network is undirected then f (X,Y,Z)=1 if and only if there
E
is no path between X and Y in D−Z.
Theorem3.1. WeconsiderstochasticSIRdynamicsdefinedonatime-invariant
network where the initial conditions are such that the states of individual nodes
are initially statistically independent. Let ψA,ψB and ψC be indicator random
X Y Z
variables or events where X,Y,Z ⊆ V are disjoint and nonempty. If Z is dy-
namically partitioningwithrespecttoX andY, andallnodes insubsystemstate
C are susceptible (C =S ∀i∈Z), then provided that hψCi=6 0,
i Z
hψAψCihψBψCi
hψAψBψCi= X Z Y Z . (6)
X Y Z hψCi
Z
Proof. If the infection has not passed through Z (which is guaranteed by all
nodes in state C being susceptible), the states of the individuals in X are sta-
tistically independent of the states of the individuals in Y. This is true since
f (X,Y,Z)=1 implies that there are no individuals from which both a mem-
E
berofX anda memberofY canbe reachedwithouttraversingamember ofZ.
We have:
P(ψA,ψB,ψC|ψC)=P(ψA,ψC|ψC)P(ψB,ψC|ψC).
X Y Z Z X Z Z Y Z Z
Given that P(ψC)6=0, we have:
Z
P(ψA,ψB,ψC) P(ψA,ψC)P(ψB,ψC)
X Y Z = X Z Y Z ,
P(ψC) P(ψC) P(ψC)
Z Z Z
from which the result follows.
Remark. For the case of zero denominator, note that P(ψC) = 0 implies that
Z
P(ψA,ψB,ψC)=0.
X Y Z
NoticethatwemadenoassumptionsabouttheSIRdynamicsinprovingthis
theorem and that it is therefore not restricted to Markoviansystems, although
it is the Markoviancase that we shallbe applying it to in the remainder of this
paper.
The theorem is a generalisation of the main result in [12] which is stated
in terms of single dynamically partitioning individuals on undirected networks.
In that context they are referred to simply as partitioning individuals due to
their correspondence to graph partitioning. Some examples of where the ex-
act closure theorem can be applied are shown in Figure 1. In this Figure and
throughout the remainder of the paper, network links without arrowheads de-
note undirected links whereas those with arrowheads denote directed links.
Figure 1a is typical of the dynamical partitioning we shall consider in this
paper. Applying Theorem 3.1, we see that there is dynamical partitioning
about node 2, so we have hI S S I S i = hI S ihS S I S i/hS i. For Fig-
1 2 3 4 5 1 2 2 3 4 5 2
ure 1b we can dynamically partition about a cluster of susceptible nodes. In
6
a) b) c)
I1 I1 I2 S1 I2
S S S S S S
5 2 3 4 5 6 I S
3 4
I4 S3 I7 I8 S5 I6
I
S 7
9
Figure1: Threeexamplesofnetworkstateswherethelocationofthesusceptible
nodes allows the application of the exact closure theorem. Here directed links
have arrowheads and undirected links do not.
fact there are two exact closures we can write down: hI I S S S S I I S i =
1 2 3 4 5 6 7 8 9
hI I S S S ihS S S S I I S i/hS S S i=hI I S S S S ihS S S I I S i/hS S S i.
1 2 3 4 6 3 4 5 6 7 8 9 3 4 6 1 2 3 4 5 6 3 4 5 7 8 9 3 4 5
InFigure1cwecanapplytheexactclosuretheoremtoobtainhS I S S I I i=
1 2 4 5 6 7
hS I S ihS S I I i/hS i. Note that I is not included in this closure.
1 2 4 4 5 6 7 4 3
For our purposes, we are interested in a special case of the exact closure
theorem which is captured by the following corollary.
Corollary 3.1. For subsystem state A of ψ , if A =S where k ∈W, and if
W k
f (n,W \k,k)=1 where n∈V \W, then
E
hψAihS I i
hψAI i= W k n . (7)
W n hS i
k
This corollary is illustrated by the example in Figure 1a. By applying this
7
to equation 4 we obtain:
hψ˙Ai = (A =S) − T (A =I)hψAi
W 1 k " kn1 n W #
k∈W n∈W
X X
+ (A =I) T (A =I)hhS(ψA)i−γ hψAi
k kn n k W k W
1 " 1 #
k∈W n∈W
X X
− (A =S) T 1−f (n,W \k,k) hhS(ψA)I i
k kn E k W n
1 "
kX∈W n∈XV\W (cid:16) (cid:17)
hhS(ψA)ihS I i
+f (n,W \k,k) k W k n
E
hS i
k #
+ (A =I) T 1−f (n,W \k,k) hhS(ψA)I i
k kn E k W n
1 "
kX∈W n∈XV\W (cid:16) (cid:17)
hhS(ψA)ihS I i
+f (n,W \k,k) k W k n . (8)
E
hS i
k #
For an arbitrary network, by applying equation 8 to the indicator random
variables S and I for all i∈{1,2,...,N}, and then reapplying it to every new
i i
subsystemstatethatemerges,aclosedsetofdifferentialequationsforthe exact
time-evolution of the probability of an individual being in a particular state is
obtained for all individuals. The number of equations that will be needed is
limited by the closures that are made possible by the exact closure theorem.
Definition 3.3. For a given network, the induced set M of subsystem states
E
is obtained by applying equation 8 to every individual (for states S and I)in the
network, and then reapplying to every new subsystem state that emerges. M
E
is then the full set of subsystem states that emerge during this process.
Remark. It follows that S and I (∀i ∈ V) and S I (∀i,j ∈ V : T > 0)
i i i j ij
represent members of M for any network.
E
3.2 Examples
Before determining the network structures under which dynamical partitioning
occurs more generally,we consider some examples. For further examples in the
context of undirected networks the reader is directed to [12].
3.2.1 Example 1
Definition 3.4. A network is a tree network if and only if its underlying graph
(all directed edges are replaced by undirected edges) is a tree or forest.
8
a) b) c)
6 5
4
1 4
1 4 5
1 2
2 3
3
2 3
Figure 2: Some example graphs. For dynamics on these graphs, we assume a
generic removal rate g and a transmission rate of 1 across all links.
Theorem3.2. ForMarkovian SIRdynamics onatreenetworkwherethestates
of all individuals are initially statistically independent, the following equations
hold exactly:
N
˙
hS i = − T hS I i,
i ij i j
j=1
X
N
˙
hI i = T hS I i−γ hI i,
i ij i j i i
j=1
X
N N
hS S ihS I i hS I ihS I i
˙ i j j k i k i j
hS I i = T − T ,
i j jk ik
hS i hS i
j i
k=1,k6=i k=1,k6=j
X X
−T hS I i−γ hS I i,
ij i j j i j
N N
hS S ihS I i hS S ihS I i
˙ i j i k i j j k
hS S i = − T − T . (9)
i j ik jk
hS i hS i
i j
k=1,k6=j k=1,k6=i
X X
Proof. For such tree networks, every individual is dynamically partitioning rel-
ative to any two of its neighbours on the underlying graph. Hence, the above
system follows directly from repeated application of equation 8, starting with
hS i, hI i ∀i∈V.
i i
Remark. Thisisthepairwisemodelthatwasshowntobeexactfortreenetworks
in [11].
3.2.2 Example 2
Consider the graph in Figure 2a. Let us suppose that all nodes have the same
removal rate g and that the transmission rate across all network links is unity.
For simplicity we shallalso makethis assumption throughthe remainderof the
explicit examples in this paper. We canapply Corollary3.1 which is embedded
inequation8tobuilduptheinducedsubsystemstatesM . Letusjustconsider
E
the infectious probability of node 1 to see how this works. We have:
˙
hI i=hS I i−ghI i.
1 1 2 1
9
Hereandthroughoutthepaperweordernodesaccordingtothenumericalorder
oftheirlabels;therelevantmotifstructuresneedtobeunderstoodwithreference
to the associated graph. Now, node 2 is dynamically partitioning with respect
to nodes 1 and 3, and it is also dynamically partitioning with respect to nodes
1 and 4. Hence:
˙
hS I i = hS S I i+hS S I i−hS I i−ghS I i
1 2 1 2 3 1 2 4 1 2 1 2
hS S ihS I i hS S ihS I i
1 2 2 3 1 2 2 4
= + −(1+g)hS I i.
1 2
hS i hS i
2 2
Rather than a complete analysis of all induced subsystem states that arise, we
take the single pair state S I from this equation as an example. Here, node
2 3
3 is not dynamically partitioning with respect to nodes 2 and 4 but node 2 is
dynamically partitioning with respect to 1 and 3 so:
hI S ihS I i
˙ 1 2 2 3
hS I i=hS S I i− −hI S I i−(1+g)hS I i.
2 3 2 3 4 4 2 3 2 3
hS i
2
Then for hS S I i, node 2 is dynamically partitioning with respect to node 1
2 3 4
and nodes 3 and 4 so:
hI S ihS S I i
˙ 1 2 2 3 4
hS S I i=− −(2+g)hS S I i.
2 3 4 2 3 4
hS i
2
Weseethathere,M representsasignificantdimensionalreductioninthenum-
E
ber of induced subsystem states compared to the full set of induced subsystem
states M.
3.2.3 Example 3
For the undirected graph in Figure 2b there is dynamical partitioning about
node 1. Starting with (for example) the infectious probability for node 1, we
have:
˙
hI i=hS I i+hS I i+hS I i+hS I i−ghI i,
1 1 2 1 4 1 5 1 6 1
where again we are assuming transmission rates of unity and a removal rate g
for each node. Now, choosing the first of these pairs to develop one part of the
induced set M gives:
E
hS I ihS I i hS I ihS I i
˙ 1 2 1 5 1 2 1 6
hS I i = hS S I i−hS I I i− −
1 2 1 2 3 1 2 4
hS i hS i
1 1
−(1+g)hS I i, (10)
1 2
and then for the first of these triples:
hS S I ihS I i hS S I ihS I i
˙ 1 2 3 1 5 1 2 3 1 6
hS S I i = hS S S I i−hS S I I i− −
1 2 3 1 2 3 4 1 2 3 4
hS i hS i
1 1
−(1+g)hS S I i.
1 2 3
10