Change points, memory and epidemic spreading in temporal networks
Tiago P. Peixoto Laetitia Gauvin
physics.soc-ph cs.SI physics.data-an





Department of Mathematical Sciences and Centre for Networks
and Collective Behaviour, University of Bath, Claverton Down, Bath BA2
7AY, United Kingdom
ISI Foundation, Via Chisola 5, 10126 Torino, Italy

ISI Foundation, Via Chisola 5, 10126 Torino, Italy




  Dynamic networks exhibit temporal patterns that vary across different
  time scales, all of which can potentially affect processes that take
  place on the network. However, most data-driven approaches used to
  model time-varying networks attempt to capture only a single
  characteristic time scale in isolation — typically associated with
  the short-time memory of a Markov chain or with long-time abrupt
  changes caused by external or systemic events. Here we propose a
  unified approach to model both aspects simultaneously, detecting
  short and long-time behaviors of temporal networks. We do so by
  developing an arbitrary-order mixed Markov model with change points,
  and using a nonparametric Bayesian formulation that allows the Markov
  order and the position of change points to be determined from data
  without overfitting. In addition, we evaluate the quality of the
  multiscale model in its capacity to reproduce the spreading of
  epidemics on the temporal network, and we show that describing
  multiple time scales simultaneously has a synergistic effect, where
  statistically significant features are uncovered that otherwise would
  remain hidden by treating each time scale independently.


Change points, memory and epidemic spreading in temporal networks
    Laetitia Gauvin
    
=================================================================





§ INTRODUCTION


Recent advances in the study of network systems — usually with social,
technological and biological origins — have been moving beyond the
more traditional approach of considering them as static or growing
entities, and instead have been introducing more realistic descriptions
that allow them to change arbitrarily in
time <cit.>. This effort includes
modeling of the time-varying network structure <cit.>, as well as processes
that take place on this dynamic environment, such as epidemic
spreading <cit.>. Further recent
works <cit.>
have highlighted the role of memory, burstiness and time ordering as key
features of empirical temporal networks that affect dynamical processes
taking place on it.

Most approaches, however, rely on a
characteristic time scale on which they describe the dynamics. These can
be divided, roughly, into approaches that model temporal correlations
via Markov chains relating short-time memory with future
behavior <cit.>,
and those that model the dynamics at longer times, usually via network
snapshots  <cit.>
or discrete change points <cit.>. In reality,
however, most systems exhibit both kinds of dynamics, and focusing on a
single aspect comes at the expense of ignoring the other. In this work,
we introduce a data-driven modeling approach that includes both aspects
simultaneously, and is capable of uncovering both the short-time Markov
properties as well a the long-time abrupt changes.

We develop a Bayesian formulation that allows both the change points and
the Markov order to be inferred from data in a principled manner,
prevents overfitting and enables model selection. As an extraneous
evaluation of our approach, we investigate the behavior of epidemic
spreading both in the original data and in artificial ones generated
from our inferred models. We show that the most plausible models tend to
mix both short-time memory and many change points, and those tend to
capture well the nontrivial epidemic behavior observed in the original
data. Importantly, the inferred models with change points typically
uncover higher-order memory than the simpler stationary variants,
demonstrating that the mixed approach is more powerful than considering
individual ones in isolation.

This paper is divided as follows. In Sec. <ref> we present
the epidemic models that will be used for the model comparison. In
Sec. <ref> we describe our modeling and inferring approach,
and apply it to empirical data.  In Sec. <ref> we
finalize with a conclusion.



§ RESULTS




 §.§ Proximity networks and epidemic dynamics


In the interest of simplicity, we will consider a minimal model of
temporal networks and epidemic dynamics that takes place on it. The
most central simplification we will make is that the dynamics takes
place in discrete time, so that the placement of edges forms a temporal
sequence, where only one edge is placed at any given time. Real
dynamical networks and epidemic spreading occur in continuous time, but
our objective here is not to construct a detailed realistic model, but
rather to illustrate how multiple time scales can be described
simultaneously. More realistic features can then be added to the model
at a later stage.

More specifically, we consider temporal networks that occur as a
sequence of edges ß={x_t}, where x_t=(u,v)_t is an edge between
nodes u and v observed at time t, with t={1,2,…,E} where
E is the number of edges. Although this formulation is general, we
focus in particular on proximity networks, obtained by tracking
volunteers with wearable sensors over a period of
time <cit.>,
so that an edge (u,v)_t is recorded if the respective people came
closer than a given radius at time t. Data recorded in this manner
possess enough time resolution for our analysis, and also serve as a
plausible scenario for epidemic
spreading <cit.>.

In the above scenario, we assume that an infection can only occur at
time t over the current “active” edge (u,v)_t. If the epidemics
follows the Susceptible-Infected-Recovered (SIR) model, and σ_u(t)
∈{S,I,R} is the state of node u at time t,
we have at each time step t:


  * If (u, v)_t is the current edge, with (σ_u(t-1),
  σ_v(t-1)) = (S, I) or (I, S), the infection spreads with
  probability β, so that (σ_u(t), σ_v(t)) = (I,
  I).

  * For every infected node u with σ_u(t-1)=I, it becomes
      recovered σ_u(t)=R with probability γ. 

The parameters β and γ control the infection and recovery
rates, respectively. We also consider the
Susceptible-Infected-Susceptible (SIS) model, which is a variation of
the above, where in the second step the infected nodes become susceptible,
σ_u(t)=S, instead of recovered. In both cases, we consider the
total number of infected nodes at given time t,

    X(t) = ∑_uδ_σ_u(t),I

where δ is the Kronecker delta.
For any positive recovery rate γ > 0, the long-time behavior of
the SIR model is always lim_t→∞X(t) = 0, as the outbreak
invariably dies out, whereas in the SIS model it can have an indefinite
permanence. In the following, we will use the behavior of X(t) as a
proxy for the comparison between data and model in capturing the
underlying network dynamics.

When considering epidemics on dynamical networks, there are two
properties that are believed to be crucial for the spreading
process <cit.>: 1. The distribution of number of contacts per link,
i.e. the frequency of token x in sequence ß, and 2. The
distribution of waiting (or inter-event) times, i.e. the time between
two occurrences of the same edge. We will have these two aspects in
mind when elaborating our models.



 §.§ Models for temporal networks


Our objective is to construct a generative model for temporal networks
that includes both short-term memories and abrupt change points. We
begin by formulating a stationary version, without change points, and
show how it is insufficient to capture many features in the data. We
then extend the model to include change points, and perform a
comparison.



  §.§.§ Stationary Markov chains

We consider sequences of discrete tokens, i.e. edges, ß={x_t} with
t∈{1,…,E} being the time and x_t∈{1,…,D} the
set of unique edges with cardinality D, which are generated from a
stationary Markov chain of order n, i.e. they occur with probability

    P(ß|,n) = ∏_t p_x_t,_t-1 = ∏_x,p_x,^a_x,,

where  corresponds to the transition matrix and p_x_t,_t-1
is the probability of observing token x_t given the previous n
tokens _t-1={x_t-1,…,x_t-n} in the sequence, and
a_x, is the number of observed transitions from memory  to
token x. This serves a simple model for temporal networks, where each
possible token corresponds to an edge in the network, i.e. x_t ≡
(i,j)_t, as we considered previously. Despite its simplicity, this
model is able to reproduce arbitrary edge frequencies, determined by the
steady-state distribution of the tokens x, as well as causal temporal
correlations between edges. This means that the model should be able to
reproduce properties of the data that can be attributed to the distribution
of number of contacts per link, which are believed to be important for
epidemic spreading <cit.>. However, due to its Markovian nature, the
dynamics will eventually forget past states, and converge to the
limiting distribution (assuming the chain is ergodic and
aperiodic). This latter property means that the model should be able to
capture nontrivial statistics of waiting times only at a short time
scale, comparable to the Markov order.




Given the above model, the simplest way to proceed would be to infer
transition probabilities from data using maximum likelihood,
i.e. maximizing Eq. <ref> under the normalization constraint
∑_xp_x,=1. This yields

    p̂_x, = a_x,/k_.

where k_=∑_xa_x, is the number of transitions originating
from .
However, if we want to determine the most appropriate Markov order
n that fits the data, the maximum likelihood approach cannot be used,
as it will overfit, i.e. the likelihood of Eq. <ref>
will increase monotonically with n, favoring the most complicated
model possible, and thus confounding statistical fluctuations with
actual structure. Instead, the most appropriate way to proceed is to
consider the Bayesian posterior distribution

    P(n|ß) = P(ß|n)P(n)/P(ß),

which involves the integrated marginal
likelihood <cit.>

    P(ß|n) = ∫ P(ß|,n)P(|n) ,

where the prior probability P(|n) encodes the amount of knowledge we
have on the transitions  before we observe the data. If we possess no
information, we can be agnostic by choosing a uniform prior

    P(|n) = ∏_(D-1)!δ(1-∑_xp_x,),

which assumes that all transition probabilities are equally likely.
Inserting Eq. <ref> and <ref> in
Eq <ref>, and calculating the integral we obtain

    P(ß|n) = ∏_(D-1)!/(k_+D-1)!∏_x a_x,!.

The remaining prior, P(n), that represents our a
priori preference to the Markov order, can also be chosen in an
agnostic fashion in a range [0,N], i.e.

    P(n) = 1/N+1.

Since this prior is a constant, the upper bound N has no effect on the
posterior of Eq. <ref>, provided it is
sufficiently large to include most of the distribution.


Differently from the maximum-likelihood approach described previously,
the posterior distribution of Eq. <ref> will select
the size of the model to match the statistical significance available,
and will favor a more complicated model only if the data cannot be
suitably explained by a simpler one, i.e. it corresponds to an
implementation of Occam's razor that prevents overfitting.

When applying this approach to empirical data, we observe that it favors
n=0 for all datasets we considered (not shown), indicating that a
higher-order model is not statistically justified, as can be seen in
Fig. <ref>.  However, if we generate temporal networks
from the fitted models,
i.e. sequence of edges using the transition probabilities
p̂_x, = a_x,/k_, they exhibit epidemic dynamics that
are very different from what we observe on the empirical time-series, as
can be seen in Fig. <ref>: for the original data,
the epidemic spreading is marked by abrupt changes in the infection
rate, which are not reproduced by the model for any value of Markov
order n — even those that overfit. Therefore, these patterns in the
epidemic dynamics seem to stem from changes in the underlying structure
of the temporal network that are not captured by the above Markov
model. Among other things, this means that the behavior cannot be
explained by a heterogeneous distribution of edge frequencies, as this
is well described by the model. As we show in the next section, the
situation changes considerably once we generalize the model to
incorporate heterogeneous Markov chains with change points.




  §.§.§ Markov chains with change points




We attempt to model the abrupt changes observed in the previous section
by non-stationary transition probabilities p_x,x⃗ that change
abruptly at a given “change point,” but otherwise remain constant
between change points. The occurrence of change points is governed by
the probability q that one is inserted at any given time. The
existence of M change points divide the time series into M+1
temporal segments indexed by l∈{0,…,M}. The variable l_t
indicates to which temporal segment a given time t belongs among the
M segments. Thus, the conditional probability of observing a token x
at time t in segment l_t is given by

    P(x_t, l_t|_t-1,l_t-1) = p_x,^l_t[q(1-δ_l_t,l_t-1)+(1-q)δ_l_t,l_t-1],

where p_x,^l_t is the transition probability inside segment l_t and
q is the probability to transit from segment l to l+1.  The
probability of a whole sequence ß={x_t} and ={l_t} being
generated is then

    P(ß, |,q) = q^M(1-q)^E-M∏_l,x,(p_x,^l)^a_x,^l

where a_x,^l is the number of transitions from
memory  to token x in the segment l.
Note that we recover the stationary model of Eq. <ref> by
setting q=0. The maximum-likelihood estimates of the parameters are

    p̂_x,^l = a_x,^l/k_^l,   q̂ = M/E

where k_^l=∑_xa_x,^l is the number of transitions originating
from  in a segment l.
But once more, we want to infer the model the segments  in a
Bayesian way, via the posterior distribution

    P(|ß,n) = P(ß, |n)/P(ß|n),

where the numerator is the integrated likelihood

    P(ß, |n) = ∫ P(ß, |,q,n)P(|n)P(q)   q

using uniform priors P(q) = 1, and

    P(|n) = ∏_lP(_l|d_l,n)P(d_l),

with the uniform prior

    P(_l|d_l,n) = ∏_(D_l-1)!δ(∑_xp_x,^l-1).

and

    P(d_l) = 2^-D

being the prior for the alphabet d_l of size D_l inside segment l,
sampled uniformly from all possible subsets of the overall alphabet of
size D. Performing the above integral, we obtain

    P(, |n) =  2^-D(M+1)M!(E-M)!/(E+1)!×
    ∏_l∏_(D_l-1)!/(k_^l+D_l-1)!∏_x a_x,^l!.


Like with the previous stationary model, both the order and the
positions of the change points can be inferred from the joint posterior
distribution

    P(,n|) = P(, |n)P(n)/P(),

in a manner that intrinsically prevents overfitting. This constitutes a
robust and elegant way of extracting this information from data, that
contrasts with non-Bayesian methods of detecting change points using
Markov chains that tend to be more
cumbersome <cit.>, and is more versatile than
approaches that have a fixed Markov order <cit.>.

The exact computation of the posterior of Eq. <ref>
would require the marginalization of the above distribution for all
possible segments , yielding the denominator P(|n), which is
unfeasible for all but the smallest time series. However, it is not
necessary to compute this value if we sample  from the posterior
using Monte Carlo. We do so by making move proposals →' with
probability P('|), and accepting it with probability a
according to the Metropolis-Hastings criterion <cit.>

    a = min(1,P('|,n)P(|')/P(|,n)P('|)),

which does not require the computation of P(|n) as it cancels out in
the ratio. If the move proposals are ergodic, i.e. they allow every
possible partition  to be visited eventually, this algorithm will
asymptotically sample from the desired posterior. Here we use the
following move proposal scheme, choosing between one the following
actions with equal probability:

  
  * We select a segment randomly and split it
        in a random point in the middle.
  
  * We merge two adjacent segments.
  
  * We move a randomly chosen boundary to a
        random position between the two enclosing ones.

We perform this algorithm many times, starting from a single segment,
and waiting sufficiently long for equilibration — determined by
observing if the likelihood value no longer changes significantly —
and we choose the partition with the largest probability across
runs. For all datasets we investigated, we observed a fast convergence
of this algorithm, which typically shows very little variation between
runs.


Note that it is also possible to change the Markov order during the
algorithm, by proposing moves n→ n', and using the
Metropolis-Hastings criterion to accept or reject them. However, we
found that Markov order typically settles very early in the algorithm,
and no longer changes during the remaining run, as it incurs a
macroscopic change in the likelihood. Since changing the Markov order is
an expensive operation of order O(E), we have found it is best to
leave it fixed during the MCMC, and select it later according to the
likelihood value.

Once a fit is obtained, we can compare the above model with the
stationary one by computing the posterior odds ratio

    Λ = P(,n|)/P(_0,n_0|) = P(,|n)/P(,_0|n_0),

where _0 is the partition into a single interval (which is
equivalent to the stationary model). A value Λ>1
[i.e. P(,|n) > P(,_0|n_0)] indicates a larger evidence for
the nonstationary model.  As can be seen in Fig. <ref>,
we observe indeed a larger evidence for the nonstationary model for all
Markov orders. In addition to this, using this general model we identify
n=1 as the most plausible Markov order, in contrast to the n=0
obtained with the stationary model. Therefore, identifying change points
allows us not only to uncover patterns at longer time scales, but the
separation into temporal segments enables the identification of
statistically significant patterns at short time scales as well, which
would otherwise remain obscured with the stationary model — even
though it is designed to capture only these kinds of correlations.





The improved quality of this model is also evident when we investigate
the epidemic dynamics, as shown in
Fig. <ref>. In order to obtain an estimate of
the number of infected based on the model, we generated different
sequences of edges using the fitted segments and transition
probabilities p̂_x,^l = a_x,^l / k_^l in each of the
segments estimated with Markov orders going from 0 to 3. We
simulated SIR and SIS processes on top of the networks generated and
averaged the number of infected over many instances. Looking at
Fig. <ref>, we see that the inferred positions
of the change-points tend to coincide with the abrupt changes in
infection rates, which show very good agreement between the empirical
and generated time-series. For higher Markov order, the agreement
improves, although the improvement seen for n>1 is probably due to
overfitting, given the results of Fig. <ref>. We note
also that the fact that n=0 provides the worse fit and agreement with
epidemic dynamics shows that it is not only the existence of change
points, but also the inferred Markov dynamics that contribute to the
quality of the model in reproducing the epidemic spreading.


In order to examine the link between the structure of the network and
the change points, we fitted a layered hierarchical degree-corrected
stochastic block
model <cit.> to the
data, considering each segment as a separate edge layer. From the figure
Fig. <ref>) we can see that the density of
connections between node groups vary in a substantial manner, suggesting
that change point marks an abrupt transition in the typical kind of
encounters between students — representing breaks between classes,
meal time, etc (see Fig. <ref>). This yields an
insight as to why these changes in pattern may slow down or speed up an
epidemic spreading: if students are confined to their classrooms,
contagion across classrooms is inhibited, but as soon they are free to
move around the school grounds, so can the epidemic.

We explore further the match between data and model by measuring the
distribution of waiting times between temporal edges, i.e. the time
interval between the occurrence in the time series of the same edge in
the network, shown in Fig. <ref> for both Markov models. For
the empirical dataset, the waiting time distribution shows a
characteristic peak at short times, and a broad decay for longer ones.
For the stationary model, the distributions obtained with the fitted
models show significant discrepancy — for both long and short times
— except when the Markov order is increased to n=3, which, according
to our Bayesian analysis cannot be used as an explanation for the data,
as it represents an overfit. However, for the nonstationary model with
change points, we observe a fair agreement between data and model for
the most-likely model with n=1, across all time scales. The
nonstationary model also provides an explanation to the shape of the
distribution at longer times, which shows a separation of time scales
inside individual stationary segments, from larger ones across change
points (marked as vertical line in Fig. <ref>). In addition
to this, the fact that the n=0 model does not reproduce the short time
behavior of the distribution shows that the Markov property inside each
stationary segment is indeed a necessary ingredient of the model. The
model that best fits the data is able to reproduce with a quite good
degree of approximation the distribution of waiting times, across all
time scales. This point is in agreement with previous results
highlighting the importance of the heterogeneity of inter-event times
for dynamical processes <cit.>, but here we see how two different time
scales are sufficient to reproduce a large fraction of the observed
behavior.

In appendix <ref> we show that the same behavior is
obtained for a variety of different datasets.



§ DISCUSSION


In this work we presented a data-driven approach to model temporal
networks that is based on the simultaneous description of the network
dynamics in two time scales: 1. The occurrence of the edges according to
an arbitrary-order Markov chain, 2. The abrupt transition of the Markov
transition probabilities at specific change-points. We developed a
Bayesian framework that allows the inference of the change points and
Markov order from data in manner that prevents overfitting, and enables
the selection of competing models.

We have applied our approach to a variety of empirical proximity
networks, and we have evaluated the inferred models based on their
capacity to reproduce the epidemic spreading observed with the original
data. We have seen that the nonstationary model accurately reproduces
the highly-variable nature of the infection rate, with changes
correlating strongly with the inferred change points. Furthermore, we
showed that the inferred model also accurately reproduces the waiting
time statistics in the empirical data, both at small and large time
scales, neither of which are accurately captured if the different time
scales are analyzed in isolation.

We argue that, ultimately, the incorporation of such temporal
heterogeneity is indispensable for the evaluation of the speeding up or
slowing down of processes taking place on dynamic
networks <cit.>, and
the development of mitigating strategies against
epidemics <cit.>.

Although our model successfully captures key properties of real dynamic
networks, it can still be made more realistic in a variety of ways. For
instance, it can be extended to continuous time via the incorporation of
waiting time distributions between events, as done in
Ref. <cit.>. Furthermore, it remains also to be
seen how the approach presented here can be extended to scenarios where
edges are allowed both to appear and disappear from the network, so that
its dynamics can no longer be represented simply by a sequence of
edges. And lastly, it would be desirable to provide a more direct
connection between the edge probabilities and change points with
large-scale network descriptors, such as community structure.


  The authors thank Ciro Cattuto, André Panisson and Anna Sapienza for useful
  conversations.








§ OTHER DATASETS


Here we show that very similar results to those described in the main
text are also encountered for other proximity datasets. In
Fig. <ref> (I) we show the analysis for the temporal behavior
of students in a primary school <cit.>,
which shows a very clear correlation of the change in infection rate and
the inferred change points. If we inspect the network structure inside
each temporal segment, we see that amounts to periods of time where the
students are either confined into classes, or mingling in larger
groups. A similar behavior is seen if Fig. <ref> (II) for
people (staff and patients) in a hospital
ward <cit.>.




