# Analytic solution of Markovian epidemics without re-infections
on heterogeneous networks

Massimo A. Achterberg<EMAIL_ADDRESS>Faculty of Electrical
Engineering, Mathematics and Computer Science, Delft University of Technology,
P.O. Box 5031, 2600 GA Delft, The Netherlands  Piet Van Mieghem  Faculty of
Electrical Engineering, Mathematics and Computer Science, Delft University of
Technology, P.O. Box 5031, 2600 GA Delft, The Netherlands

###### Abstract

Most epidemic processes on networks can be modelled by a compartmental model,
that specifies the spread of a disease in a population. The corresponding
compartmental graph describes how the viral state of the nodes (individuals)
changes from one compartment to another. If the compartmental graph does not
contain directed cycles (e.g. the famous SIR model satisfies this property),
then we provide an analytic, closed-form solution of the continuous-time
Markovian compartmental model on heterogeneous networks. The eigenvalues of
the Markovian process are related to cut sets in the contact graph between
nodes with different viral states. We illustrate our finding by analytically
solving the continuous-time Markovian SI and SIR processes on heterogeneous
networks. We show that analytic extensions to e.g. non-Markovian dynamics,
temporal networks, simplicial contagion and more advanced compartmental models
are possible. Our exact and explicit formula contains sums over all paths
between two states in the SIR Markov graph, which prevents the computation of
the exact solution for arbitrary large graphs.

##  I Introduction

Most research in network epidemiology focuses on so-called compartmental
models, in which the population is subdivided into several groups based on the
current viral state of the individual [1]. The SIS and SIR models are the most
elementary compartmental models, in which (S) represents a susceptible, but
healthy individual, (I) represents an infected and infectious individual and
(R) represents a resistant or removed individual, who is no longer able to
contract the disease. Both SIS and SIR are composed of two processes: (i)
infected individuals can infect susceptible neighbours and (ii) infected
individuals can recover, either becoming susceptible again (SIS) or receiving
permanent immunity (SIR).

The individuals can be represented by a set $\mathcal{N}$ of $N$ nodes in a
contact graph $G$. The set $\mathcal{L}$ of $L$ links describes the
connections between nodes with a certain weight. The elements $a_{ij}$ of the
$N\times N$ adjacency matrix $A$ denote whether there is a link ($a_{ij}=1$)
between nodes $i$ and $j$ or not ($a_{ij}=0$). We assume that the contact
graph $G$ is connected, such that any pair of nodes is connected by a path in
the contact graph.

General epidemic models are comprised of $c$ compartments. We assume that
interactions only occur between neighbouring nodes and that after such an
interaction, at most one of those nodes changes its viral state. Then the
transitions between compartments can be classified into two types [2]: link-
based transitions, such as the infection of a susceptible node by a
neighbouring infected node and node-based transitions, such as nodes
recovering from the disease. Many results in epidemics have been obtained in
the well-mixed limit, i.e. all individuals can directly contact all other
individuals [3]. In reality, each individual has different characteristics and
cannot contact any other member of the population, leading to heterogeneous
behaviour during an epidemic [1]. Therefore, we focus in this work on epidemic
processes on heterogeneous networks.

Properties of such network-based models are usually computed by mean-field
approximations [4] or Monte Carlo simulations [5]. The downside is that Monte
Carlo simulations are known to converge slowly, which is problematic if the
probability of a certain event is small. Simulating is especially tricky
around the epidemic threshold, which is precisely the regime in which most
real-world epidemics occur. On the other hand, the mean-field approximation
works well for homogeneous, well-mixed contact graphs but is less precise for
heterogeneous networks [6]. Instead of resorting to either Monte Carlo
simulations or mean-field approximations, we focus in this work on determining
the _exact_ time-dependent solution for general Markovian compartmental
models.

For the computation the exact time-dependent solution, we require an efficient
enumeration of all states of the underlying Markov chain. Since each of the
$N$ nodes can attain $c$ viral states, the number of states equals $c^{N}$.
The set of all possible states is denoted by $\mathbf{S}$. Providing explicit
constructions for the states111Throughout this work, we will use _state_ and
_configuration_ when referring to one element of the state space $\mathbf{S}$,
interchangeably. and the corresponding infinitesimal generator on
heterogeneous networks for general Markovian compartmental models is a tedious
task. Several researchers have attempted to construct infinitesimal generators
for specific epidemic models on networks. For example, Van Mieghem _et al_.
[7] provided a construction based on binary numberings for SIS epidemics and
for SIS epidemics with self-infections [8]. Simon _et al_. [9] suggested a
tridiagonal block structure for the infinitesimal generator of the SIS
process. Each block contains configurations with the same number of infected
nodes $k$. Configurations with $k$ infected nodes can only make transitions to
configurations with $k+1$ and $k-1$ infected nodes, such that the
infinitesimal generator $Q$ has a block-tridiagonal structure. Economou _et
al_. [10] developed a similar tridiagonal block structure, but used a
different ordering within each block. For SIR epidemics, López-García [11]
groups configurations based on the number of recovered nodes. Within each
block with the same number of recovered nodes, configurations are grouped
based on the number of infected nodes. Within each block with the same number
of recovered _and_ infected nodes, a lexicographical ordering is applied. The
Generalised Epidemic Mean-Field (GEMF) model by Sahneh _et al_. [2] describes
disease transmission in general compartmental models, where the state space is
constructed using a tensor or Kronecker product formulation. Merbis [12] used
with a similar tensor product construction to derive the time-dependent
equations for general epidemic models. In this tensor formulation, Merbis and
Lodato [13] derived an exact solution for the SI process on unweighted graphs.
Another representation in terms of a tensor-product formulation for SIR
epidemics is provided by Dolgov and Savostyanov [14].

Mappings of a temporal network to a directed, static graph have been
established. An example of such a directed, static graph is the event graph
[15], whose nodes represent events (e.g. corresponding to the creation or
removal of links in the contact graph) and its links represent possible
follow-up events. The space of all states in Markovian epidemics has a similar
structure, although nodes and links have a different meaning. We believe that
our solution method is superior, because it can also include temporal
networks, whereas an extension of the event graph to handle epidemic spreading
is not straightforward.

General compartmental epidemic models, as formalised in the GEMF framework
[2], can be solved by the eigendecomposition of the corresponding
infinitesimal generator. If the compartmental graph is a tree, i.e. the
compartmental graph does not contain any cycles, then the corresponding Markov
graph is also a tree, whose infinitesimal generator $Q$ only has upper-
triangular non-zero elements. Consequently, the eigenvalues and eigenvectors
can be determined efficiently. If cycles occur, then the infinitesimal
generator $Q$ is no longer upper-triangular and our elegant solution method is
not applicable, because the eigenvalues and eigenvectors can only be
determined numerically, which is infeasible for large graphs.

In this work, we present for the first time an analytic closed-form solution
of transient continuous-time Markovian epidemic processes, i.e. without re-
infections, on heterogeneous networks. We demonstrate our method by
considering SI epidemics in Section II. We present an in-depth explanation of
the eigenvalues of the infinitesimal generator in Section II.1 and derive the
full solution in Section II.2. We discuss the benefits and difficulties of the
exact solution in Section II.3 and compare our exact results with simulations
in Section II.4. We then generalise our exact solution on non-Markovian
dynamics, temporal networks, simplicial contagion and self-infections in
Section II.5. We also provide an extensive analysis on SIR epidemics in
Section III: the SIR eigenvalues in Section III.1, SIR eigenvectors in Section
III.2 and the full solution in Section III.3. We demonstrate the power of our
exact method by analytically computing several properties in Section III.4-I.
Finally, we conclude in Section IV.

##  II The Markovian SI process

The Susceptible-Infected (SI) process describes the spread of information,
diseases, innovations or neural activity over a network. Starting from a
single infected node, the infection spreads to all its neighbours, which again
infects their neighbours, until the whole network is infected [16]. The
Markovian SI process is actually a Markovian discovery process, or,
equivalently, a stochastic shortest path process [17, Chapter 16]. The SI
process describes spreading phenomena without damping or curing (for example,
HIV viruses in humans or the spread of innovations) or the spread evolves
extremely fast, effectively flooding the population and prohibiting curing or
recovering (for example, epileptic seizures in the human brain [18] and the
spread of information on social media).

The SI process subdivides the nodes into two compartments: the set
$\mathcal{I}\subseteq\mathcal{N}$ of infected (I) nodes are infected, aware or
activated and the complementary set ${\cal S}={\cal N}\backslash{\cal I}$ of
susceptible (S) nodes are healthy, unaware or idle. We further assume that an
infection over the direct link from node $i$ towards a susceptible node $j$ is
a Poisson process with rate $\beta_{ij}$ and is independent of all other
links. The associated compartmental graph is shown in Figure 1a. The SI
process is a continuous-time Markov process containing $2^{N}$ states, because
each of the $N$ nodes in the contact graph is either susceptible or infected.
The $2^{N}\times 1$ state vector is
$\mathbf{s}(t)=(s_{0}(t),s_{1}(t),\ldots,s_{2^{N}-1}(t))^{T}$, whose element
$s_{i}(t)$ describes the probability that the SI process is in state $i$ at
time $t$. Since the process must certainly be, for any time $t$, in one of the
possible states, $\sum_{i=0}^{2^{N}-1}s_{i}(t)=1$ for all times $t$. The
transition rates between the states are denoted by the $2^{N}\times 2^{N}$
infinitesimal matrix $Q$.

(a)  Compartmental graph (b)  Contact graph (c)  Markov graph

Figure 1: The spread of SI epidemics involves three networks: (a) the
compartmental graph specifies the transitions between the different
compartments, (b) the nodes in the contact graph $G$ represent individuals and
the links specify the connections among the individuals and (c) the
combination of the SI process from (a) and the contact network from (b) leads
to the Markov graph containing all possible configurations and the transitions
between those configurations. Grey nodes indicate absorbing states. The number
of compartments is 2 and the number of nodes $N$, such that the Markov graph
is of size $2^{N}$. In this example, $N=4$.

We describe the SI process using the labeling in [7]. The idea of the
construction in [7, 8] is that the state $i$ represents the viral state of $N$
nodes in the contact graph. We define the binary variables $x_{k}(i)$
indicating whether node $k$ is infected in configuration $i$. Then the _viral
state vector_ $\mathbf{x}=(x_{N},x_{N-1},\ldots,x_{2},x_{1})^{T}$ describes
the viral state of all nodes. By regarding the viral state vector $\mathbf{x}$
as a binary number, the configuration number $i$ is computed as

$i=\sum_{k=1}^{N}x_{k}(i)2^{k-1},$

where $x_{k}=1$ if node $k$ is infected and $x_{k}=0$ if node $k$ is
susceptible. For a complete contact graph with $N=4$ nodes, the possible
transitions for the SIS process are shown in Figure 1. For example, the all-
healthy state $i=0$ has representation $(0000)$, indicating that all nodes in
state $0$ are healthy. Similarly, for state $5$ the binary representation is
$(0101)$, indicating that node 2 and 4 are healthy and node 1 and 3 are
infected.

A fundamental property of Markov processes is that events in continuous time
occur sequentially and simultaneous events do not happen almost surely. This
means that transitions between states in the Markov chain can only occur if
the binary representation of those states differs exactly one bit, which is
illustrated in Figure 1c. However, not all states differing one bit are
connected in Figure 1c. The reason is that the SI process only considers an
infection process: upward-pointing transitions in the Markov graph in Figure
1c cannot occur, because upward transitions correspond to the recovery of a
node. Furthermore, if the SI process is in the all-healthy state, then the
spread never activates, because there is no initial infection. Hence, the all-
healthy state $i=0$ cannot be entered nor left, thus state $i=0$ has no
incoming nor outgoing arrows in the SI process. Finally, we remark that the
contact graph $G$ in Figure 1b describes the existence of links between pairs
of nodes. Naturally, the existence and weight of the links in the Markov graph
are influenced by the existence of links in the contact graph $G$.

Using the binary notation of the state space $\mathbf{S}$, we describe the
transitions between the states using the $2^{N}\times 2^{N}$-dimensional
infinitesimal generator $Q$. To simplify the notation, we introduce the
$N\times N$, possibly asymmetric, weighted and non-negative adjacency matrix
$B$ with elements $\tilde{\beta}_{kl}=a_{kl}\beta_{kl}$. The elements $q_{ij}$
of the infinitesimal generator $Q$ for the SI process are [17, Chapter 17]

$\displaystyle q_{ij}$
$\displaystyle=\sum_{k=1}^{N}\tilde{\beta}_{mk}x_{k}(i),\qquad\text{if
}j=i+2^{m-1}$ (1a) $\displaystyle\qquad\qquad\text{ with }m=1,2,\ldots,N\text{
and }x_{m}(i)=0.$ $\displaystyle q_{ii}$
$\displaystyle=-\sum_{\begin{subarray}{c}j=0\\\ j\neq
i\end{subarray}}^{2^{N}-1}q_{ij},$ (1b)

for $i,j=1,2,\ldots,2^{N}$. The non-zero elements of the infinitesimal
generator $Q$ are visualised in Figure 2. The infinitesimal generator $Q$ is
an upper triangular matrix, because in absence of curing events, there is no
transition from state $i$ to state $j<i$ because the number of infected nodes
is larger for state $i$ than for state $j$. The infinitesimal generator $Q$
for a contact graph with $N$ nodes can also be constructed recursively in
terms of the infinitesimal generator of a graph with $N-1$ nodes and then
adding one node [8].

Figure 2: The structure of the infinitesimal generator $Q$ for the SI process
on a complete graph with $N=7$ nodes.

The governing equation for the continuous-time Markov chain is

$\frac{d\,\mathbf{s}^{T}(t)}{dt}=\mathbf{s}^{T}(t)Q,$ (2)

whose solution is

$\mathbf{s}^{T}(t)=\mathbf{s}^{T}(0)e^{Qt}.$

The solution can be further detailed as [7, p. 4]

$\mathbf{s}(t)=\pi+\sum_{i=1}^{2^{N}-1}e^{\lambda_{i}t}\sum_{m=0}^{n_{i}-1}\mathbf{r}_{i,m}\frac{t^{m}}{m!}.$
(3)

where $n_{i}$ denotes the multiplicity of eigenvalue $\lambda_{i}$,
$\pi=(0,0,\ldots,0,1)^{T}$ is the steady-state vector and the vector
$\mathbf{r}_{i,m}$ is related to the right- and left-eigenvectors of $Q$ and
the initial condition $\mathbf{s}(0)$. Unfortunately, the infinitesimal
generator $Q$ is asymmetric and is not even normal222An $N\times N$ matrix $A$
is _normal_ if it commutes with its conjugate transpose: $AA^{*}=A^{*}A$.,
preventing further simplifications in Eq. (3).

The SI process is completely described by the set (2) of $2^{N}$ linear
differential equations. Frequently, an exponentially large state space is the
fingerprint of a non-linear process. In particular, SIS and SIR epidemics
exhibit a phase transition around the epidemic threshold. Below the threshold,
the epidemic dies out exponentially fast, whereas above threshold, a non-zero
fraction of the population remains infected over a long time. On the contrary,
the SI process does not exhibit a phase transition, because the SI process
always converges – given that the process starts with at least one infected
node – to the all-infected state, irrespective how small, but non-zero, the
infection rates $\beta_{ij}>0$ are. For simple graphs with homogeneous
transition rates333Homogeneous infection rates are $\beta_{kl}=\beta$ for all
nodes $k,l$, such that $B=\beta A$., like the complete graph or the star
graph, an exact analysis is possible [19] by exploiting symmetry of the
contact graph [9]. Otherwise, exactly solving equation (2) or (3) for the SI
process is infeasible for large $N>20$ and one resorts to Monte Carlo
simulations [20, 5].

The general solution (3) contains the eigenvalues and eigenvectors of the
infinitesimal generator $Q$. The eigenvalues $\lambda_{i}$ of the
infinitesimal generator $Q$ are of primary importance in the solution (3),
because the eigenvalues are inversely proportional to the relaxation time of
the corresponding eigenmode. The eigenvalue $\lambda_{i}$ is therefore related
to the average time for the SI process to make a transition from configuration
$i$ to another configuration $j$. The next section further elaborates on the
eigenvalues of the SI process.

###  II.1 Eigenvalues of the infinitesimal generator $Q$

The $2^{N}$ eigenvalues of the $2^{N}\times 2^{N}$ infinitesimal generator $Q$
of a Markov chain have negative real part and at least one eigenvalue is zero
[17]. Since the infinitesimal generator of the SI process is upper triangular
as illustrated in Figure 2, the eigenvalues of the SI process appear on the
diagonal of the infinitesimal generator $Q$, hence,
$\lambda_{i}=q_{ii}=-\sum_{j=i+1}^{2^{N}}q_{ij}$ for all
$i=0,1,\ldots,2^{N}-1$. Due to the triangular structure of the infinitesimal
generator $Q$, the eigenvalue $i$ can be understood as to belong to
configuration $i$. Prior to describing the physical significance of the
eigenvalue $\lambda_{i}$, we introduce some notation.

Each configuration $i$ describes which nodes in the network are infected and
which nodes are susceptible. For a given configuration $i$, we partition the
set of nodes $\mathcal{N}$ into two groups. The group
$\mathcal{I}_{i}\subseteq{\cal{N}}$ defines the set of infected nodes in
configuration $i$. Similarly, $\mathcal{S}_{i}$ contains the susceptible nodes
in configuration $i$. Each node is either infected or susceptible, so
$\mathcal{N}=\mathcal{I}_{i}\cup\mathcal{S}_{i}$ for all configurations $i$.
We define the _cut set_ as the set of links with one node in $\mathcal{S}_{i}$
and one node in $\mathcal{I}_{i}$.

The eigenvalue $\lambda_{i}$ is equal to the sum over all possible transitions
from configuration $i$ to any other configuration $j$, where configuration $j$
differs exactly one bit (one node) from configuration $i$. Starting in
configuration $i$, visualised in Figure 3, the possible configurations $j$ are
those where one susceptible node will be infected by one of the already
infected nodes. Hence, the eigenvalue $\lambda_{i}$ equals (minus) the sum
over all weighted links in the S-I cut set;

$\lambda_{i}=-\sum_{k\in\mathcal{I}_{i}}\sum_{l\in\mathcal{S}_{i}}\tilde{\beta}_{kl},$
(4)

or, following the notation of [17];

$\lambda_{i}=-\sum_{k=1}^{N}\sum_{l=1}^{N}x_{k}(i)(1-x_{l}(i))\tilde{\beta}_{kl}.$
(5)

If the contact graph $G$ has homogeneous link weights such that
$\tilde{\beta}_{kl}=\beta a_{kl}$, where $a_{kl}$ specifies whether a link
exists between node $k$ and $l$, then $\lambda_{i}/\beta$ _equals (minus) the
number of links in the S-I cut set_.

Equation (4) provides the exact eigenvalue $\lambda_{i}$ for every
configuration $i$. The time-dependent solution (3) shows the implication of
the eigenvalues $\lambda_{i}$. If the eigenvalue is large (in modulus), the
corresponding state will converge exponentially fast to zero. For small (in
modulus) eigenvalues, the convergence is much slower. The SI process can be
regarded as a discovery process, that starts with a set of discovered nodes
and discovers adjacent nodes. If $t\to\infty$ on a connected graph, all nodes
will be discovered. Large (in modulus) eigenvalues have a small contribution
to the total discovery time, but small (in modulus) eigenvalues will have a
significantly larger contribution to the total discovery time in the SI
process. In particular, a small (in modulus) eigenvalue $\lambda_{i}$ of a
configuration $i$ corresponds to a small weight of the S-I cut set (see Figure
3) and it takes a long time to transfer the disease between these two groups
of nodes. The minimal eigenvalue (equivalent to finding the minimum weighted
cut in the contact graph) can be obtained efficiently using the Stoer-Wagner
algorithm [21].

Figure 3: The cut set in SI epidemics on a contact graph with 9 nodes. The
eigenvalue $\lambda_{i}$ that belongs to configuration $i$ is (minus) the sum
over all weighted links in the S-I cut set. The links in the cut set are
illustrated by dashed lines.

Finally, equation (4) illustrates that eigenvalue $\lambda=0$ occurs twice;
the all-healthy state (where $\mathcal{I}_{i}=\emptyset$) and the all-infected
state (where $\mathcal{S}_{i}=\emptyset$) have an empty cut set and the
corresponding eigenvalue is zero. Hence, the all-healthy state and the all-
infected are both steady states of the SI process. The all-healthy state is
unstable, because adding a single infected node will lead to more infected
nodes, whereas the all-infected state is stable.

Given a configuration $i$, then the number of infected nodes $k$ in
configuration $i$ can be calculated as

$k=\sum_{l=1}^{N}x_{l}(i).$ (6)

For several contact graphs, the eigenvalues can be computed analytically.

For the complete graph with homogeneous infection rates, the eigenvalues are

$\lambda=-\beta k(N-k),\qquad\text{with multiplicity }\binom{N}{k}$

where $0\leq k\leq N$ is the number of infected nodes in configuration $i$.
The SI process on the complete graph with homogeneous infection rates and
starting with a fixed number $k$ of infected nodes is presumably the easiest
SI process, because in this case, the SI process can be transformed to a birth
and death process, whose time-dependent solution can be calculated exactly
[17, Ch. 16]. Starting from $k$ infected nodes and ending with $l$ infected
nodes, the total spreading time is distributed as the sum of independent
exponential random variables, with parameters $\lambda_{m}$ where $k\leq m\leq
l-1$ represents the number of infected nodes [16, 22].

For the star graph with homogeneous infection rates, the eigenvalues with
multiplicity $\binom{N-1}{k}$ are

$\lambda=\begin{cases}-\beta(N-1-k),\qquad&\text{if hub infected}\\\ -\beta
k,\qquad&\text{if hub healthy}\end{cases}$

where $0\leq k\leq N-1$ is the number of infected leaf nodes in configuration
$i$. In this case, the eigenvalues can be simplified to the form
$\lambda=-\beta k$ with multiplicity $2\binom{N-1}{k}$ for $0\leq k\leq N-1$.

For the cycle graph with homogeneous infection rates, the eigenvalues are

$\lambda_{i}=-2\beta p(k,i),$

where $p(k,i)$ depends on the number of susceptible neighbours of the $k$
infected nodes in configuration $i$. Contrary to the complete graph and the
star graph, the position of the infected nodes is crucial in the cycle graph.
Starting in configuration $i$ with $k$ infected nodes, the $k$ infected nodes
can be grouped into $p\leq k$ connected, infectious components. The total
number of susceptible neighbours is $2p$, because each infected component
connects to two susceptible neighbours. After infecting one of the susceptible
neighbours, the number of connected, infected components either stays constant
or reduces by one, because two components will merge if the in-between
susceptible node has been infected. In that case, the number of components
reduces to $p-1$. The number of connected, infected components can therefore
only decrease over time. Since the eigenvalue $\lambda_{i}$ only decreases (in
modulus) over time, the discovery speed of the SI process also slows down over
time.

In the Erdős-Rényi graph [23] with homogeneous infection rates with weight
$\beta=1$, links exist independent of other links with probability $p$. The
eigenvalue probability distribution $\textnormal{Pr}[\lambda=-l]$ of a random
configuration $i$ follows from the law of total probability;

$\displaystyle\textnormal{Pr}[\lambda=-l]$
$\displaystyle=\sum_{k=0}^{N}\textnormal{Pr}\left[\lambda=-l\ |\ \text{state
}i\text{ has }k\text{ infected nodes}\right]$
$\displaystyle\qquad\quad\times\textnormal{Pr}[\text{state }i\text{ has
}k\text{ infected nodes}].$

Suppose that configuration $i$ consists of $k$ infected nodes. Then the
eigenvalue $\lambda_{i}=-l$ is (minus) the number of links in the cut set of
configuration $i$. The cut set contains maximally $k(N-k)$ links, where each
link exists independently of the other links with probability $p$. Hence, the
probability distribution is a Binomial distribution with $k(N-k)$ possible
links and probability $p$;

$\displaystyle\textnormal{Pr}[\lambda=-l\ |\ \text{state }i\text{ has }k\text{
infected nodes}]$ $\displaystyle\quad=\binom{k(N-k)}{l}p^{l}(1-p)^{k(N-k)-l}.$

The probability for configuration $i$ to have $k$ infected nodes is equal to
the probability for any graph with $N$ nodes to have $k$ infected nodes, which
equals $2^{-N}\binom{N}{k}$. Combining all, we find

$\displaystyle\textnormal{Pr}[\lambda=-l]=$
$\displaystyle\left(\frac{p}{1-p}\right)^{l}2^{-N}\sum_{k=0}^{N}\binom{k(N-k)}{l}\binom{N}{k}(1-p)^{k(N-k)},$
(7)

with the convention that $\binom{k(N-k)}{l}$ is zero if $l<0$ or $l>k(N-k)$.
Hence, the eigenvalue $\lambda$ is bounded between $\lambda=0$ and
$\lambda=-k(N-k)$. We plot the exact solution (II.1) with numerical
simulations in Figure 4, which shows excellent agreement.

Figure 4: The eigenvalue distribution of the infinitesimal generator $Q$ of
the SI process in an Erdős-Rényi graph with $N=10$ nodes and link-existence
probability $p=0.36$ averaged over 100,000 simulations. The exact solution
(II.1) is shown as a red dashed line. The simulation results do not follow the
exact solution for $\lambda\leq-20$ because the number of simulations is too
small compared to the probability of occurrence of these eigenvalues.

###  II.2 The exact SI solution

An explicit form of the SI solution (3) requires the computation of the
eigenvectors of the infinitesimal generator $Q$. We denote by
$\mathbf{v}_{0},\mathbf{v}_{1},\ldots,\mathbf{v}_{2^{N}-1}$ and
$\mathbf{w}_{0},\mathbf{w}_{1},\ldots,\mathbf{w}_{2^{N}-1}$ the right- and
left-eigenvectors of $Q$, respectively. The left-eigenvector
$\mathbf{w}_{2^{N}-1}$, corresponding to the all-infected state $i=2^{N}-1$
with eigenvalue $\lambda_{2^{N}-1}=0$, is equal to the steady-state vector
$\pi$:

$\mathbf{w}_{2^{N}-1}=\pi=(0,0,\ldots,0,1)^{T}.$ (8)

The corresponding right-eigenvector is the all-ones vector
$\mathbf{v}_{2^{N}-1}=(1,1,\ldots,1)^{T}$. Since there are no transitions
possible to and from the all-healthy state $i=0$ in the SI process, _we remove
the all-healthy state from the Markov chain_ , which reduces the Markov chain
to $2^{N}-1$ states.

The SI process for homogeneous transitions rates was recently solved by Merbis
and Lodato [13], thus we focus on heterogeneous infection rates
$\tilde{\beta}_{kl}$. In many practical applications, the weighted infection
rates $\tilde{\beta}_{kl}$ are real numbers and we can safely assume the
following:

###### Assumption 1

The eigenvalues of the infinitesimal generator $Q$ are distinct.444The
eigenvalues may also be degenerate, as long as the algebraic multiplicity
equals the geometric multiplicity for all eigenvalues. Then the matrix is
diagonalisable and Eq. (9) holds. It is, however, unclear when the algebraic
and geometric multiplicities are equal for general contact graphs.,555Suppose
that the weighted infection rates $\tilde{\beta}_{kl}$ are real and
independently distributed, the probability that two eigenvalues are equal is
almost surely zero on a finite graph.,666All contact graphs with homogeneous
transition rates do not satisfy Assumption 1. We know that in homogeneous
contact graphs, the eigenvalue equals the number of links in the cut set.
Consider one infected node and all other nodes susceptible, then the size of
the cut set exactly equals the infected node’s degree. Since every connected
graph has at least two nodes with the same degree [24], at least one
eigenvalue is degenerate.

Following Assumption 1, the infinitesimal generator $Q$ is diagonalisable and
the left- and right eigenvectors from different eigenvalues are orthogonal:
$\mathbf{w}_{i}^{T}\mathbf{v}_{j}=\delta_{ij}$ for all $1\leq i,j\leq
2^{N}-1$, where $\delta$ indicates the Kronecker delta function [24]. Then the
solution $\mathbf{s}(t)$ in Eq. (3) simplifies to

$\mathbf{s}(t)=\pi+\sum_{i=1}^{2^{N}-1}c_{i}e^{\lambda_{i}t}\mathbf{w}_{i}$
(9)

where the constant $c_{i}=\mathbf{s}(0)^{T}\mathbf{v}_{i}$. Moreover, any
infinitesimal generator $Q$ has row sum zero and thus a zero eigenvalue with
corresponding all-one right-eigenvector $\mathbf{v}_{2^{N}-1}=\mathbf{u}$.
Hence, all left-eigenvectors $\mathbf{w}$ are orthogonal to the right-
eigenvector $\mathbf{v}_{2^{N}-1}=\mathbf{u}$, which implies that each left-
eigenvector (except $\mathbf{w}_{2^{N}-1}$) sums element-wise to zero [24,
art. 140–142].

The set of nodes that are initially infected at time $t=0$ remain infected for
all times $t>0$, because nodes cannot cure. Consider a configuration $i$ in
the Markov graph. If one of the initially infected nodes is susceptible in
configuration $i$, then configuration $i$ cannot be reached from the initial
state and configuration $i$ can be removed from the Markov graph.

Using Assumption 1 and exploiting the upper-triangular structure of the
infinitesimal generator $Q$, we can explicitly compute the eigenvectors. We
say that state $i$ can _reach_ state $j$ if there is a directed path in the
Markov graph from state $i$ to state $j$. If state $i$ can reach state $j$,
then we denote $i\to j$.

Let us first consider the right-eigenvectors $\mathbf{v}_{i}$. For
configurations $i$ that correspond to one infected node, its right-eigenvector
is the basis vector $\mathbf{e}_{i}$. For configurations $i$ with two infected
nodes, the right-eigenvector will have non-zero elements at the positions that
correspond to all states that can reach state $i$ (including state $i$
itself). Let us consider an example with $N=3$ nodes in Figure 5.

(a)  (b)

Figure 5: (a) Example of a weighted contact graph with $N=3$ nodes and (b) the
corresponding Markov graph. The grey nodes can be removed from the Markov
graph, because the SI process starts with node 1 infected (indicated by the
black node) and all grey nodes indicate configurations where node 1 is
susceptible.

###### Example 1

The contact graph in Figure 5a consists of 3 nodes. Node 1 is initially
infected whereas node 2 and 3 are susceptible. The Markov graph in Figure 5b
is comprised of only four states: (00), (01), (10) and (11), where
$(x_{3}x_{2})$ represent the viral state of nodes 3 and 2, respectively. We
have omitted the viral state of node 1, as node 1 is always infected. The
infinitesimal generator $Q$ equals

$Q=\begin{pmatrix}q_{00}&q_{01}&q_{02}&0\\\ 0&q_{11}&0&q_{13}\\\
0&0&q_{22}&q_{23}\\\ 0&0&0&0\end{pmatrix}$

where $q_{00}=-3,q_{01}=1,q_{02}=2,q_{11}=-6,q_{13}=6,q_{22}=-5,q_{23}=5$. For
configuration 0=(00) with eigenvalue $\lambda_{0}=q_{00}=-3$, the right-
eigenvector $\mathbf{v}_{0}$ equals the basis vector
$\mathbf{e}_{0}=(1,0,0,0)^{T}$. For configuration $1=(01)$, the eigenvector
can be computed using the definition;

$\displaystyle(Q-q_{11}I)\mathbf{v}_{1}=0\quad\Longleftrightarrow$
$\displaystyle\quad\left(\begin{array}[]{cccc|c}q_{00}-q_{11}&q_{01}&q_{02}&0&0\\\
0&0&0&q_{13}&0\\\ 0&0&q_{22}-q_{11}&q_{23}&0\\\ 0&0&0&-q_{11}&0\\\
\end{array}\right)$

The last three rows reveal that $(\mathbf{v}_{1})_{2}=(\mathbf{v}_{1})_{3}=0$.
Now we choose $(\mathbf{v}_{1})_{1}=1$. Then the solution follows as
$(\mathbf{v}_{1})_{0}=\frac{q_{01}}{q_{11}-q_{00}}=-\frac{1}{3}$. We conclude
that $\mathbf{v}_{1}=(-1/3,1,0,0)^{T}$ and the only non-zero elements in
$\mathbf{v}_{1}$ correspond to configurations that can reach configuration 1.
Analogously, we can find $\mathbf{v}_{2}=(-1,0,1,0)^{T}$. For
$\mathbf{v}_{3}$, we obtain

$(Q-q_{33}I)\mathbf{v}_{3}=0\quad\Longleftrightarrow\quad\left(\begin{array}[]{cccc|c}q_{00}&q_{01}&q_{02}&0&0\\\
0&q_{11}&0&q_{13}&0\\\ 0&0&q_{22}&q_{23}&0\\\ 0&0&0&0&0\\\ \end{array}\right)$

We choose $(\mathbf{v}_{3})_{3}=1$. Then the elements
$(\mathbf{v}_{3})_{2}=-\frac{q_{23}}{q_{22}}$ and
$(\mathbf{v}_{3})_{1}=-\frac{q_{13}}{q_{11}}$. The element
$(\mathbf{v}_{3})_{0}$ can be computed in the same manner and equals
$(\mathbf{v}_{3})_{0}=-\frac{q_{01}(\mathbf{v}_{1})_{1}+q_{02}(\mathbf{v}_{1})_{2}}{q_{00}}$.
Here, the iterative nature of the construction is clear;
$(\mathbf{v}_{3})_{0}$ depends on $(\mathbf{v}_{3})_{1}$ and
$(\mathbf{v}_{3})_{2}$. $\hfill\square$

In general, the right-eigenvector $\mathbf{v}_{i}$ of a certain configuration
$i$ can be constructed iteratively;

$\displaystyle(\mathbf{v}_{i})_{j}$ $\displaystyle=0,\qquad\qquad\text{if
state }j\text{ cannot reach state }i$ $\displaystyle(\mathbf{v}_{i})_{i}$
$\displaystyle=1,\qquad\qquad\text{by construction}$
$\displaystyle(\mathbf{v}_{i})_{h}$
$\displaystyle=\frac{q_{ih}}{\lambda_{i}-\lambda_{h}}=\frac{\displaystyle\sum_{l=1}^{N}\tilde{\beta}_{ml}x_{l}(h)}{\lambda_{i}-\lambda_{h}},$
$\displaystyle\qquad\qquad\qquad\text{if }i-h=2^{m},m=0,1,\ldots,N-1$
$\displaystyle(\mathbf{v}_{i})_{g}$
$\displaystyle=\sum_{h=0}^{2^{N}-1}\frac{q_{hg}}{\lambda_{i}-\lambda_{g}}(\mathbf{v}_{i})_{h}$
$\displaystyle\qquad\qquad\qquad\text{if }h-g=2^{n},n=0,1,\ldots,N-1$
$\displaystyle\vdots$

Configurations $h$ and $i$ differ only at position $m$, which corresponds to
node $m+1$ being infected in configuration $i$ and susceptible in
configuration $h$. Similarly, configuration $g$ and $h$ differ only at
position $n$ corresponding to node $n+1$, etc. The iterative formulation
consists of at most $N-1$ steps, because at each iteration, one node turns
from susceptible into infected, there are $N$ nodes in the contact graph and
we start with at least one infected node. We emphasise that this construction
is infeasible for degenerate eigenvalues.

The left-eigenvectors $\mathbf{w}_{i}$ can be derived similarly, from which we
find

$\displaystyle(\mathbf{w}_{i})_{h}$ $\displaystyle=0,\qquad\qquad\text{if
state }h\text{ cannot be reached,}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\text{starting from state
}i.$ $\displaystyle(\mathbf{w}_{i})_{i}$ $\displaystyle=1,\qquad\qquad\text{by
construction}$ $\displaystyle(\mathbf{w}_{i})_{j}$
$\displaystyle=\frac{q_{ij}}{\lambda_{i}-\lambda_{j}}=\frac{\displaystyle\sum_{l=1}^{N}\tilde{\beta}_{ml}x_{l}(j)}{\lambda_{i}-\lambda_{j}},$
$\displaystyle\qquad\qquad\qquad\text{if }j-i=2^{m},m=0,1,\ldots,N-1$
$\displaystyle(\mathbf{w}_{i})_{k}$
$\displaystyle=\sum_{j=0}^{2^{N}-1}\frac{q_{jk}}{\lambda_{i}-\lambda_{k}}(\mathbf{w}_{i})_{j}$
$\displaystyle\qquad\qquad\qquad\text{if }k-j=2^{n},n=0,1,\ldots,N-1$
$\displaystyle\vdots$

By construction, the eigenvectors are orthonormal
$\mathbf{w}_{i}^{T}\mathbf{v}_{i}=1$ for all configurations $i$, because the
eigenvectors have only one non-zero value in common, which is the 1 at
position $i$. Using the right-eigenvectors $\mathbf{v}_{i}$ and the left-
eigenvectors $\mathbf{w}_{i}$, the solution (9) can be computed explicitly.

The element $(\mathbf{v}_{i})_{h}$ of the right eigenvector $\mathbf{v}_{i}$
is non-zero if state $h$ can reach state $i$, in other words, if there is a
directed path in the Markov graph from state $h$ to state $i$. The opposite
holds for the left-eigenvector $\mathbf{w}_{i}$, whose element $j$ is non-zero
if there is a directed path in the Markov graph from state $i$ to state $j$.

Thus, starting a certain state $i$, the element $h$ of the right-eigenvector
$\mathbf{v}_{i}$ specifies whether the process can start in configuration $h$
and converge to configuration $i$. Similarly, the non-zero element $j$ of the
left-eigenvector $\mathbf{w}_{i}$ specifies whether the process can start in
configuration $i$ and converge to configuration $j$.

###  II.3 Computational feasibility

Irrespective of the application domain, the quantity of interest of the SI
process should be computable in a reasonable time. The $2^{N}\times 2^{N}$
infinitesimal generator $Q$ cannot be computed nor stored if the number of
nodes $N\geq 25$. Fortunately, the exact solution $\mathbf{s}(t)$ under
Assumption 1 can be computed without explicitly constructing the infinitesimal
generator $Q$.

Only the eigenvalues and eigenvectors of the infinitesimal generator $Q$ are
required in the solution (9). For a given configuration $i$, the eigenvalue
$\lambda_{i}$ can be computed based on equation (5) in $\mathcal{O}(N^{2})$
operations. To compute all $2^{N}$ eigenvalues, we need
$\mathcal{O}(N^{2}2^{N})$ operations. The iterative procedure to compute a
single left- and right-eigenvector takes $\mathcal{O}(N2^{N})$ operations if
the eigenvalues are known. Since there are $2^{N}$ eigenvectors, the total
number of required operations to compute all eigenvectors is
$\mathcal{O}(N4^{N})$. Hence, for networks with $N>20$ nodes, computing all
eigenvalues and eigenvectors remains infeasible.

Even though the computation of all eigenvalues and eigenvectors is infeasible
for large networks, for certain quantities of interest, not all eigenvalues or
eigenvectors are required. Here, we consider two exemplary scenarios. First,
we consider the probability $\textnormal{Pr}[X_{2}(t)=1]$ that node 2 is
infected at time $t$. To compute $\textnormal{Pr}[X_{2}(t)=1]$, we sum over
all states in $\mathbf{S}$ in which node 2 is infected. We introduce the
$2^{N}\times 1$ vector $\mathbf{m}=(0,1,1,0,0,1,\ldots,1)^{T}$, which is $1$
at position $j$ if the binary representation of configuration $j$ has a 1 at
position 2, and zero otherwise. Then
$\textnormal{Pr}[X_{2}(t)=1]=\mathbf{m}^{T}s(t)$ and using (9), we find

$\textnormal{Pr}[X_{2}(t)=1]=\mathbf{m}^{T}\pi+\sum_{i=1}^{2^{N}-2}c_{i}e^{\lambda_{i}t}\mathbf{m}^{T}\mathbf{w}_{i}.$
(10)

The inner product $\mathbf{m}^{T}\pi=1$, because node 2 will be infected after
infinitely long time. We know that each left-eigenvector $\mathbf{w}_{i}$ sums
element-wise to zero. Then for all configurations $i$ for which the non-zero
elements of the left-eigenvector $\mathbf{w}_{i}$ overlap with the non-zero
elements of $\mathbf{m}$, the inner product is zero. Only the states $i$ that
can be reached from state $2$ will have a zero inner product, which are
$2^{N-1}$ states. Since $c_{i}$ is non-zero for $2^{N-1}$ states, we conclude
that the sum in (10) can be simplified from $2^{N}-1$ states to $2^{N-2}$
states. Unfortunately, summing over $2^{N-2}$ states is still exponentially
large and solving for large networks remains impossible.

As a second example, we consider the probability that all nodes are infected
at time $t$. We multiply the solution $\mathbf{s}(t)$ with the vector
$\mathbf{m}=(0,0,\ldots,0,1)^{T}$, which is the all-zeros vector, except the
last element is one. Regarding the product $m^{T}\mathbf{w}_{i}$ in (10), we
conclude that only the last element $(\mathbf{w}_{i})_{2^{N}-1}$ is of
importance. Unfortunately, computing the last element of the left-eigenvector
$\mathbf{w}_{i}$ is as difficult as computing the whole vector, because of the
iterative construction of the eigenvector $\mathbf{w}_{i}$.

Another key property of epidemics involves the average fraction of infected
nodes, also known as the prevalence. Computation of the prevalence requires
the computation of the probability of infection of each node individually,
which is definitely computationally harder than the previous two examples.

The provided examples consider limit cases, namely the infection probability
of a single node and all nodes, respectively. For both examples, exponentially
many eigenvalues and eigenvectors are required to build up the solution (9).
The examples illustrate that calculating the exact, time-varying solution of
the infection probability in SI epidemics on large contact graphs is generally
infeasible.

###  II.4 Numerical simulations

We illustrate the accuracy of our exact solution method in Figure 6 for two
contact graphs: an Erdős-Rényi graph with $N=10$ nodes and link-connectivity
$p=0.33$ and a cycle graph with $N=12$ nodes. The simulations start with node
1 infected. The results are averaged over 10,000 Monte Carlo simulations. We
compare our exact solution and the Monte Carlo simulations with the
N-Intertwined Mean-Field Approximation (NIMFA), which assumes every pair of
random variables is uncorrelated [26]. The exact solution coincides with the
Monte Carlo simulations, whereas the mean-field approximation vastly
overestimates the time-dependent prevalence [27].

(a)  (b)

Figure 6: The exact solution (solid line), Monte Carlo simulations (asterisks)
and the mean-field approximation (dashed line) of the Markovian SI process on
(a) an Erdős-Rényi graph with $N=10$ nodes and $p=0.33$ and (b) a cycle graph
with $N=12$ nodes. The infection rates $\beta_{ij}$ are chosen uniformly at
random between $\beta_{\min}=0.1$ and $\beta_{\max}=1.1$. Initially, node 1 is
infected. The mean-field approximation deviates in both cases significantly
from the exact solution.

As we concluded earlier, the large (in modulus) eigenvalues have a small
effect on the dynamics at large times. Presumably, the exact solution (9) can
be accurately approximated using a truncation method. The number of modes is
truncated at index $m$, such that the approximated solution
$\tilde{\mathbf{s}}$ equals

$\tilde{\mathbf{s}}(t)\approx\pi+\sum_{i=1}^{m}c_{i}e^{\lambda_{i}t}\mathbf{w}_{i}$
(11)

The approximation of Eq. (9) by Eq. (11) introduces the error

$e(t)=\|\mathbf{s}(t)-\tilde{s}(t)\|=\left\|\sum_{i=m+1}^{2^{N-1}}c_{i}e^{\lambda_{i}t}\mathbf{w}_{i}\right\|\leq\sum_{i=k+1}^{2^{N-1}}e^{\lambda_{i}t}\|c_{i}\mathbf{w}_{i}\|$

where $\|\cdot\|$ is a vector norm and we used the triangle inequality
$\|\mathbf{a}+\mathbf{b}\|\leq\|\mathbf{a}\|+\|\mathbf{b}\|$. Using the fact
that the eigenvectors are normalised, we find

$e(t)\leq\sum_{i=m+1}^{2^{N-1}}e^{\lambda_{i}t}\|c_{i}\mathbf{w}_{i}\|<e^{\lambda_{m+1}t}\sum_{i=m+1}^{2^{N-1}}|c_{i}|$

such that the error scales as $\mathcal{O}(e^{\lambda_{m+1}t})$. Similarly to
$\varepsilon$-SIS dynamics [25, Appendix D], the truncation method performs
poor for small times, see Figure 7. The main reason is as follows. Suppose the
SI epidemic starts at a single node that has only one link with a very small
weight. Under the truncation approximation, the eigenmode corresponding to the
infection from the seed node to the second node will be disregarded, whereas
this link has a profound impact on the dynamics.

Figure 7: The exact solution (solid line) and the truncated solution (dotted
line) with $m=N^{2}$ out of $2^{N}$ states of the Markovian SI process on an
Erdős-Rényi graph with $N=8$ nodes. The infection rates $\beta_{ij}$ are
chosen uniformly at random between $\beta_{\min}=0.9$ and $\beta_{\max}=1.1$.
Initially, node 1 is infected.

###  II.5 Extensions of Markovian SI epidemics

For Markovian SI processes on heterogeneous networks, we have derived the
explicit solution (9). To demonstrate the versatility of our method, we extend
our results on the SI process to non-Markovian epidemics, temporal networks,
simplicial contagion and the inclusion of self-infections.

####  II.5.1 Non-Markovian dynamics

Most epidemic models assume a memoryless Markov process, i.e. the probability
to infect a neighbour is exponentially distributed over time. Non-Markovian
effects in Markov chains can be taken into account using fractional calculus,
as was recently derived by Van Mieghem [28]. The governing equations (2) of
the Markov chain change to

$D^{\alpha}\mathbf{s}_{\alpha}^{T}(t)=\mathbf{s}_{\alpha}^{T}(t)Q^{\alpha},$
(12)

where $D^{\alpha}$ is the Caputo fractional derivative and $0<\alpha<1$. As
derived by Van Mieghem [28], the solution (9) becomes

$\mathbf{s}(t)=\pi+\sum_{i=1}^{2^{N}-1}c_{i}E_{\alpha,1}\left((\lambda_{i}t)^{\alpha}\right)\mathbf{w}_{i}$
(13)

where $E_{\alpha,1}(z)$ is the Mittag-Leffler function [29] of the complex
number $z$. Up to our knowledge, Eq. (13) is the first exact solution of a
non-Markovian epidemic process on an heterogeneous network. The exact solution
(13) is plotted in Figure 8 for $0.5\leq\alpha\leq 1$. Decreasing $\alpha$
increases the prevalence at short times but slows down the disease spread for
large times. Compared to the Markovian case $\alpha=1$, decreasing $\alpha$
increases the variance of the infection time distribution, causing a large
probability of both very slow and very fast infections. At short times, this
results in a larger infection probability, but it also takes much longer for
all nodes to become infected.

Figure 8: The exact solution of the non-Markovian SI process with Caputo
fractional derivative $\alpha$ on an Erdős-Rényi graph with $N=10$ nodes. The
infection rates $\beta_{ij}$ are chosen uniformly at random between
$\beta_{\min}=0.9$ and $\beta_{\max}=1.1$. Initially, node 1 is infected. The
case $\alpha=1$ coincides with the Markovian case.

####  II.5.2 Temporal networks

So far, we have assumed that the contact network is fixed. In reality, contact
networks are varying over time due to movements of individuals. Given a
temporal contact graph $G$ that changes finitely many times in the interval
$[0,T]$, each connected contact graph $G_{n}$ in an interval $[t_{n},t_{n+1}]$
is represented by the adjacency matrix $A^{(n)}$. We assume that the contact
graph remains connected at all times. Then we may compute the eigenvectors and
eigenvalues for each of those adjacency matrices $A^{(n)}$ and the solution
becomes

$\mathbf{s}(t)=\pi+\sum_{i=1}^{2^{N}-1}c_{i}^{(n)}e^{\lambda_{i}^{(n)}t}\mathbf{w}_{i}^{(n)},\qquad
t_{n}\leq t\leq t_{n+1}$ (14)

where $c_{i}^{(n)}=\mathbf{s}(t_{n})\mathbf{v}_{i}^{(n)}$. Within each
interval $t_{n}\leq t\leq t_{n+1}$, the dynamics are equivalent to the case of
the fixed network. Going to the next time interval requires different
eigenvalues and eigenvectors, and additionally, the initial condition of the
new interval must equal the final state of the previous interval. Although the
eigenvectors and eigenvalues may be different for each time interval, each
contact graph is assumed to be connected, thus the steady-state vector
$\pi=(0,0,\ldots,0,1)^{T}$, representing the state in which all nodes are
infected.

####  II.5.3 Simplicial contagion

In the study of opinion dynamics on networks, it is observed that multiple
neighbours of a node may be required to persuade a node to adopt a particular
strategy or opinion. Such higher-order spreading phenomena on networks are
known as simplicial contagion [30], which was first applied to the SIS model
by Iacopini _et al._ [31]. The conceptual idea is that, besides node-node
interactions, connected infected triangles may increase the probability of
converting a particular node to the other opinion, even larger than the sum of
the individual infection rates. Adding simplicial contagion to the SI process
does not alter the _structure_ of the infinitesimal generator $Q$, because
simplicial infections require the existence of triangles in the original
network structure, whose individual links already had contributions to the
infinitesimal generator. Instead, the _link weights_ corresponding to
triangles are increased for certain configurations in the infinitesimal
generator $Q$. The eigenvalues are no longer the weighted sum over all links
in the S-I cut set, but additionally include the weighted sum over all
2-infected-1-susceptible-node triangles in the contact graph. The eigenvectors
can be computed in the usual way, such that the time-dependent solution (9)
can be obtained.

Figure 9 illustrates the impact of simplicial contagion on the classical SI
process. Adding the homogeneous simplicial infection rate $\beta_{\Delta}$
increases the prevalence $y$, but only up to a certain point, because some
parts of the contact graph do not contain any triangles and therefore nodes
cannot be infected by higher-order simplexes.

Figure 9: The exact solution of the SI process with simplicial contagion on an
Erdős-Rényi graph with $N=8$ nodes. The infection rates $\beta_{ij}$ are
chosen uniformly at random between $\beta_{\min}=0.2$ and $\beta_{\max}=1.0$
and the triangle interaction rate $\beta_{\Delta}$ is taken equal to $0,3$ and
$100$. Initially, node 1 is infected. The Monte Carlo simulations are shown as
red asterisks.

####  II.5.4 Self-infections

Contrary to the standard SI process, where the spread or discovery of items is
solely related to the contact graph, some spreading phenomena can be triggered
by external processes, which are unrelated to the spread on the contact graph.
For example, individuals or companies may adopt a certain product or
technology without any interference with others. The situation in which a node
becomes infected spontaneously without inference with other nodes, is
described as a self-infection with rate $\varepsilon$. The Bass model
describes the spread of direct infections and self-infections of new products
used in companies [32] and is actually equivalent to the SI process with self-
infections [8].

The infinitesimal generator $Q$ of the $\varepsilon$-SI process equals

$\displaystyle q_{ij}$
$\displaystyle=\varepsilon_{m}+\sum_{k=1}^{N}\tilde{\beta}_{mk}x_{k}(i),\qquad\text{if
}j=i+2^{m-1}$ (15) $\displaystyle\qquad\qquad\text{ with }m=1,2,\ldots,N\text{
and }x_{m}(i)=0.$ $\displaystyle q_{ii}$
$\displaystyle=-\sum_{\begin{subarray}{c}k=0\\\ k\neq
i\end{subarray}}^{2^{N}-1}q_{ki},$ (16)

Thus, adding self-infections can both increase the link weight and add more
links in the Markov graph. The resulting infinitesimal generator $Q$ will
remain upper-triangular, because the total number of infected nodes can only
increase. In this case, the eigenvalue $\lambda_{i}$ equals the weighted sum
of the links in the S-I cut set plus the total (possibly weighted) sum of the
self-infection rate of all susceptible nodes in configuration $i$.

Adding the self-infection process to the SI process increases the prevalence
significantly, as shown in Figure 10. In the limit where the self-infection
process is much larger than the infection process, each node adopts the
innovation independent of all other nodes and the prevalence $y$ is just the
sum of $N$ independent exponential random variables with rate $\varepsilon$.

Figure 10: The exact solution of the Markovian SI process with self-infections
on an Erdős-Rényi graph with $N=10$ nodes. The infection rates $\beta_{ij}$
are chosen uniformly at random between $\beta_{\min}=0.9$ and
$\beta_{\max}=1.1$. The self-infection rate $\varepsilon$ is taken equal to
$0,0.3$ and $0.7$. Initially, node 1 is infected. Figure 11: The Markov graph
of the SIR process on the complete graph with $N=3$ nodes. Each state
$(x_{3}x_{2}x_{1})$ encodes the configuration (viral state) of the nodes;
susceptible (0), infected (1) and recovered (2). The configuration number $i$
is shown in bold.

##  III The SIR process

Contrary to the SI process, whose applicability is limited, network-based SIR
epidemics has a very rich and long history in modelling network epidemiology
[1]. The SIR process is a transient process, because after being infected,
nodes will recover and can never be re-infected. We solve here the SIR process
on a heterogeneous network by proposing a labelling of the state space based
on trinary numerals.

Again, to simplify the notation, we denote the $N\times N$, possibly
asymmetric, weighted and non-negative adjacency matrix $B$ with elements
$\tilde{\beta}_{kl}=a_{kl}\beta_{kl}$.

We introduce $x_{k}$ as the viral state of node $k$:

$\displaystyle x_{k}$ $\displaystyle=0,\qquad\text{if node }k\text{ is
susceptible}$ $\displaystyle x_{k}$ $\displaystyle=1,\qquad\text{if node
}k\text{ is infected}$ $\displaystyle x_{k}$ $\displaystyle=2,\qquad\text{if
node }k\text{ is recovered}$

The viral state vector $\mathbf{x}=(x_{N},x_{N-1},\ldots,x_{1})^{T}$ describes
the viral state of all nodes. Given a particular viral state vector
$\mathbf{x}$, we may compute the configuration number $i$ using the trinary
numerals;

$i=\sum_{k=1}^{N}x_{k}(i)\cdot 3^{k-1}$

The trinary numbering ensures that each viral state vector $\mathbf{x}$
corresponds to a unique configuration number $i$, which ranges between $i=0$
and $i=3^{N}-1$. Then the infinitesimal generator $Q$ of the SIR process
equals

$\displaystyle q_{ij}$
$\displaystyle=\sum_{k=1}^{N}\tilde{\beta}_{mk}\cdot\mathbf{1}_{\\{x_{k}(i)=1\\}},\qquad\text{if
}j=i+3^{m-1}$ $\displaystyle\qquad\qquad\text{ with }m=1,2,\ldots,N\text{ and
}x_{m}(i)=0$ $\displaystyle q_{ij}$
$\displaystyle=\delta_{m},\qquad\qquad\text{if }j=i+3^{m-1}$
$\displaystyle\qquad\qquad\text{ with }m=1,2,\ldots,N\text{ and }x_{m}(i)=1$
$\displaystyle q_{ii}$ $\displaystyle=-\sum_{\begin{subarray}{c}j=0\\\ j\neq
i\end{subarray}}^{3^{N}-1}q_{ji},$

where $\mathbf{1}_{\\{x\\}}$ is the indicator function, which is one if
condition $x$ is satisfied and zero otherwise. The Markov graph specifies all
possible transitions between the $3^{N}$ states and the rates are given by the
infinitesimal generator $Q$. The Markov graph for a complete graph with $N=3$
nodes is shown in Figure 11. The absorbing states are indicated by shaded
circles. We observe that the Markov graph is a directed tree, i.e. a bipartite
graph. Moreover, we can define the _layer number_ $l(i)=I(i)+2R(i)$, where
$I(i)$ represents the number of infected nodes and $R(i)$ the number of
recovered nodes in configuration $i$. Transitions can only take place from
states in layer $l$ to states in layer $l+1$. Layers with odd number $l$ do
not contain any absorbing states, because an odd layer number $l$ implies that
at least one node must be infected. If a configuration contains one infected
node, that node can always recover, ruling out the existence of absorbing
states in that layer.

Figure 12 shows all non-zero elements of the infinitesimal generator $Q$ with
$N=4$ nodes. Some rows of $Q$ are empty, because these rows correspond to
absorbing states in the SIR process (the grey nodes in Figure 11). Although
the infinitesimal generator $Q$ has dimensions $3^{N}\times 3^{N}$, the total
number of transitions is much lower. For a given configuration $i$, the number
of possible transitions is given by $N$, because each node can change its
viral state to at most one other state. Taking into account the possibility of
not making any change, the number of non-zero elements in the infinitesimal
generator $Q$ is upper bounded by $(N+1)3^{N}$, which is much less than a
dense matrix structure with $3^{2N}$ elements.

Figure 12: The structure of the infinitesimal generator $Q$ for the SIR
process on a complete graph with $N=4$ nodes.

###  III.1 SIR eigenvalues

Similar to the SI process, the eigenvalue $\lambda_{i}$ of the SIR process can
be understood to belong to configuration $i$ because the infinitesimal
generator $Q$ is triangular. The eigenvalue $\lambda_{i}$ of configuration $i$
equals:

$\lambda_{i}=-\sum_{k=1}^{N}\sum_{l=1}^{N}\tilde{\beta}_{kl}\cdot\mathbf{1}_{\\{x_{k}(i)=1\cap\,x_{l}(i)=0\\}}-\sum_{k=1}^{N}\delta_{k}\cdot\mathbf{1}_{\\{x_{k}(i)=1\\}}$
(17)

The first component specifies the sum over the weighted links in the S-I cut
set as in (5), whereas the second part contains the total weighted curing rate
of the infected nodes in configuration $i$. Each configuration $i$ that does
not contain any infected nodes (i.e. $\mathbf{1}_{\\{x_{k}(i)=1\\}}=0$ for all
$k$) is an absorbing state and cannot be left. The total number of states
without infected nodes, thus only containing susceptible and recovered nodes,
equals $2^{N}$. Hence, there are $2^{N}$ absorbing states. The eigenvalue of
an absorbing state equals $\lambda=0$ and eigenvalue $\lambda=0$ has algebraic
multiplicity $2^{N}$. To determine the geometric multiplicity of eigenvalue
$\lambda=0$, i.e. solving $Q\mathbf{v}=\mathbf{0}$, we investigate the
structure of the infinitesimal generator $Q$ in Figure 12. There are $2^{N}$
rows with only zeros, allowing for $2^{N}$ free variables in the right-
eigenvector $\mathbf{v}$. Hence, the algebraic multiplicity is equal to the
geometric multiplicity for $\lambda=0$ and moreover, the eigenvectors with
eigenvalue zero can be chosen orthogonally.

It remains to verify the uniqueness of the non-zero eigenvalues. Similar as
for SI epidemics, we require that all infection rates $\tilde{\beta}_{kl}$ are
different. For SIR epidemics, we additionally require that the contact graph
is complete. Consider the non-complete contact graph $G$ from Figure 13 with a
missing link between the orange and black node. Suppose the black node is
infected and the orange node is either susceptible or recovered. All other
nodes are considered susceptible (represented by the cloud). The eigenvalue
$\lambda$ in (17) equals the SI cut set plus the curing rate of the black
node. The fact that the orange node is susceptible or recovered does not
influence the eigenvalue $\lambda$, because the SI cut set and the sum over
the curing rates are equivalent in both cases. Thus, any non-complete contact
graph has one or more degenerate non-zero eigenvalues.

From simulations it appears that even in the case of missing links in the
contact graph (and subsequently degenerate eigenvalues), the corresponding
eigenspace is of full rank and orthogonal eigenvectors can be found. Explicit
constructions for those eigenvectors are cumbersome and are omitted here.
Instead, we can set infection rates $\tilde{\beta}_{kl}>0$ arbitrarily small
for every non-existing link. The resulting graph is the complete graph with
heterogeneous infection rates, which can accurately approximate any realistic
situation at the cost of a small error and avoids the computation of the
complicated case of degenerate eigenvalues.

Figure 13: The contact graph $G$ is represented by a green cloud containing
susceptible nodes, an infected black node and the orange node. The link
between the black and orange node is missing. Whether the orange node is
susceptible or recovered does not influence the computation of the eigenvalue
$\lambda$ from Eq. (17), thus the eigenvalue $\lambda$ is degenerate.

###  III.2 SIR eigenvectors

The right- and left-eigenvectors $\mathbf{v}_{i}$ and $\mathbf{w}_{i}$ of the
SIR process can be constructed in a similar fashion as for the SI process in
Section II.2. We focus on the right-eigenvector $\mathbf{v}_{i}$, but the
left-eigenvector $\mathbf{w}_{i}$ can be constructed similarly.

For a configuration $i$ that corresponds to one infected node, the right-
eigenvector $\mathbf{v}_{i}$ is the basis vector $\mathbf{e}_{i}$. For a
configuration $i$ with two infected nodes, the right-eigenvector
$\mathbf{v}_{i}$ will have non-zero elements at the positions that correspond
to all states that can reach state $i$ (including state $i$ itself). In fact,
we follow the same procedure as for the SI process. The only exception occurs
when $\lambda_{i}=0$, which is a degenerate eigenvalue. Fortunately, for
$\lambda=0$ the algebraic and geometric multiplicity are equal, such that the
eigenspace is of full rank and we may choose the eigenvectors of the zero
eigenvalues as standard basis vectors: for configuration $i$ with eigenvalue
$\lambda_{i}=0$, we choose $\mathbf{v}_{i}=\mathbf{e}_{i}$.

In general, the right-eigenvector $\mathbf{v}_{i}$ of a certain configuration
$i$ can be constructed iteratively;

If $\displaystyle\mathbf{\lambda_{i}=0;}\textbf{ then
}\mathbf{v}_{i}=\mathbf{e}_{i}$ If $\displaystyle\mathbf{\lambda_{i}\neq
0;}\textbf{ then}$ $\displaystyle(\mathbf{v}_{i})_{j}=0,\qquad\qquad\text{if
state }j\text{ cannot reach state }i$
$\displaystyle(\mathbf{v}_{i})_{i}=1,\qquad\qquad\text{by construction}$
$\displaystyle(\mathbf{v}_{i})_{h}=\frac{q_{ih}}{\lambda_{i}-\lambda_{h}}$
$\displaystyle\phantom{(\mathbf{v}_{i})_{h}}=\frac{\displaystyle\sum_{l=1}^{N}\tilde{\beta}_{ml}\cdot\mathbf{1}_{\\{x_{l}(h)=1\cap\,x_{m}(h)=0\\}}+\delta_{m}\cdot\mathbf{1}_{\\{x_{m}(h)=1\\}}}{\lambda_{i}-\lambda_{h}}$
$\displaystyle\qquad\qquad\qquad\text{if }i-h=3^{m},m=0,1,\ldots,N-1$
$\displaystyle(\mathbf{v}_{i})_{g}=\sum_{h=0}^{3^{N}-1}\frac{q_{hg}}{\lambda_{i}-\lambda_{g}}(\mathbf{v}_{i})_{h}$
$\displaystyle\qquad\qquad\qquad\text{if }h-g=3^{n},n=0,1,\ldots,N-1$
$\displaystyle\ldots$

Configurations $h$ and $i$ differ only at position $m$ (corresponding to node
$m+1$), where node $m+1$ is infected in configuration $i$ and susceptible in
configuration $h$ or node $m+1$ is recovered in $i$ and infected in $h$.
Similarly, configuration $g$ and $h$ differ only at position $n$,
corresponding to node $n$, etc.

###  III.3 SIR solution

After the derivation of the orthonormal eigenvectors for the zero and non-zero
eigenvalues, the infinitesimal generator $Q$ is diagonalisable and the time-
dependent solution $\mathbf{s}(t)$ becomes

$\mathbf{s}(t)=\sum_{i=0}^{3^{N}-1}c_{i}e^{\lambda_{i}t}\mathbf{w}_{i}$ (18)

where $c_{i}=\mathbf{s}(0)^{T}\mathbf{v}_{i}$ and
$\mathbf{v}_{i},\mathbf{w}_{i}$ are the right and left-eigenvector of $Q$,
respectively. Defining the $3^{N}\times 1$ vector $\mathbf{m}$, whose elements
$m_{i}$ equal the number of infected nodes in configuration $i$, the time-
varying prevalence $y(t)$ follows as

$y(t)=\frac{1}{N}\sum_{i=0}^{3^{N}-1}c_{i}e^{\lambda_{i}t}\mathbf{m}^{T}\mathbf{w}_{i}$

which can be simplified to

$y(t)=\sum_{i=0}^{2^{N}-1}\tilde{c}_{i}+\sum_{i=0}^{3^{N}-2^{N}-1}\tilde{c}_{i}e^{\lambda_{i}t}$
(19)

where $\tilde{c}_{i}=\frac{1}{N}c_{i}\mathbf{m}^{T}\mathbf{w}_{i}$. The first
term in Eq. (19) contains all zero eigenvalues and the second term contains
all non-zero (negative) eigenvalues. Surprisingly, the time-dependent
prevalence $y(t)$ in (19) is just a function of the eigenvalues $\lambda_{i}$
and the (complicated) variables $\tilde{c}_{i}$.

We can further simplify (19). The contribution of all non-zero eigenvalues in
the second term of (19) converges exponentially fast to zero for $t\to\infty$,
whereas the first term in (19) corresponding to the zero eigenvalues, is
fixed. Since the prevalence $y\to 0$ for $t\to\infty$, it follows that the
first term must be zero. Hence, the solution (19) reduces to

$y(t)=\sum_{i=0}^{3^{N}-2^{N}-1}\tilde{c}_{i}e^{\lambda_{i}t}$ (20)

Figure 14 shows the exact solution (20) and Monte Carlo simulations, which
match perfectly. For SI epidemics, we considered non-Markovian dynamics using
fractional derivatives, included time-varying contact networks, added higher-
order simplicial contagion and added self-infections. The same extensions can
straightforwardly be applied to SIR epidemics.

Figure 14: The exact solution of the Markovian SIR process on a complete graph
with $N=7$ nodes. The infection rates $\tilde{\beta}_{kl}$ are chosen
uniformly at random between $\beta_{\min}=0.2$ and $\beta_{\max}=1.0$ and the
curing rates $\delta_{k}$ are chosen uniformly between $\delta_{\min}=0.1$ and
$\delta_{\max}=0.6$. Initially, node 1 is infected.

####  III.3.1 Explicit SIR solution

Since the Markov graph only contains directed links and re-infections are
impossible, the SIR process is transient. Written as a system of $3^{N}$
linear differential equations, the equations (2) can be solved exactly [33].
Assuming that the initial condition $\mathbf{s}(0)=0$, except for some state
$i$ where $s_{i}(0)=1$, then the exact solution follows as

$\displaystyle\begin{split}s_{j}(t)&=0,\qquad\text{if }j<i\\\
s_{i}(t)&=e^{\lambda_{i}t},\\\
s_{j}(t)&=q_{jl}\left(\frac{e^{\lambda_{j}t}-e^{\lambda_{l}t}}{\lambda_{l}-\lambda_{j}}\right)-\sum_{h=2}^{l-j}(-1)^{h}\sum_{k_{1}=j+1}^{l-1}\sum_{k_{2}=k_{1}+1}^{l-1}\cdots\\\
&\quad\sum_{k_{h-1}=k_{h-2}+1}^{l-1}q_{j,k_{1}}\left({\prod_{r=1}^{h-2}}q_{k_{r},k_{r+1}}\right)q_{k_{h-1},l}\\\
&\hskip
14.22636pt\times\sum_{n=0}^{h}\frac{e^{\lambda_{k_{n}}t}}{{\prod\nolimits_{r=0;r\neq
n}^{h}}\left(\lambda_{k_{r}}-\lambda_{k_{n}}\right)},\qquad\text{if
}j>i\end{split}$ (21)

where $h$ denotes the number of hops of a path between starting state $i$ and
final state $j$. In fact, Eq. (21) is just the non-vectorised form of (18),
directly providing intuition for the coefficients $\tilde{c}_{i}$ and left-
eigenvectors $\mathbf{w}_{i}$.

The maximum hopcount $h=2N-1$ in Eq. (21), because of the following reason.
Consider the most extreme case, where the process starts with a single
infected node and all other nodes are susceptible and the process ends with
all nodes recovered. During the process, $N-1$ infections and $N$ recoveries
must take place, leading to a total of $2N-1$ transitions. The path is denoted
by $\mathcal{P}=(i,k_{1},k_{2},\ldots,k_{h-1},j)$. Naturally, not all paths
are valid paths, because of the physical constraints of the SIR process (i.e.
nodes must be infected before they can recover). In Appendix E, we demonstrate
that the maximal number of paths $\mathcal{P}$ is at least $N!$ for any
connected contact graph $G$. Thus, in spite of having the closed-form
analytical solution (21), it is infeasible to compute the exact solution for
any moderately sized network.

Surprisingly, the Markovian SIR process with homogeneous and heterogeneous
infection and curing rates can be solved exactly on tree graphs with initially
one infected node [34, 35]. The paradoxical discrepancy is due to our
computation of the exact solution for all states, whereas [34, 35] merely
focus on global properties of the SIR dynamics, such as the prevalence and the
average fraction of recovered nodes.

###  III.4 Epidemic peak time in SIR

A key property of SIR epidemics is the _epidemic peak time_ , which is the
time at which on average most nodes are infected. Knowing when the number of
infections reaches the peak allows decision-makers to employ timely counter-
measures. If the infection rate in the SIR process is above the epidemic
threshold, then an outbreak occurs with large probability. After potentially a
long time, the disease dies out, exhibiting a peak in the number of cases at
the epidemic peak time. Due to the heterogeneous infection and curing rates,
multiple local extrema may be observed [36].

The epidemic peak can be derived explicitly for the Markovian SIR process. The
derivative of $y(t)$ in (19) equals

$y^{\prime}(t)=\sum_{i=0}^{3^{N}-2^{N}-1}\tilde{c}_{i}\lambda_{i}e^{\lambda_{i}t}.$
(22)

The epidemic peak time $t_{\text{peak}}$ obeys
$y^{\prime}(t_{\text{peak}})=0$. We determine the epidemic peak time
$t_{\text{peak}}$ using the Newton-Raphson and second-order Newton-Raphson
method (see Appendix A for the derivation). Methods like Newton-Raphson are
suitable for the determination of the epidemic peak time using our solution,
since all derivatives of the prevalence $y(t)$ can be straightforwardly
computed.

(a)  (b)

(c)  (d)

Figure 15: (a) The time-varying prevalence $y(t)$ in SIR epidemics on the
complete graph with $N=7$ nodes for 100 trials. Each trial randomly generated
uniform random variables $\tilde{\beta}_{kl}\in[0.05,0.85]$ and
$\delta_{k}\in[0.1,0.9]$. (b) The derivative of the prevalence $y$ with
respect to time $t$ for each of the 100 trials. (c-d) The error of finding the
true epidemic peak using the (c) Newton-Raphson method and the (d) second-
order Newton-Raphson method with 0–5 iterations. Initially, $\tilde{t}_{0}=0$.
The minimal error for finding the peak time is $10^{-10}$, because the “exact”
peak is also determined numerically using the bisection method. Newton-Raphson
in (c) does not converge, whereas second-order Newton-Raphson in (d) almost
always converges quickly.

Starting with initial guess $\tilde{t}_{0}$, the Newton-Raphson method finds a
root $t^{*}$ of $y^{\prime}(t)$ iteratively as

$\tilde{t}_{k+1}=\tilde{t}_{k}-\frac{y^{\prime}(\tilde{t}_{k})}{y^{\prime\prime}(\tilde{t}_{k})},$
(23)

and converges, provided that the initial point $\tilde{t}_{0}$ is sufficiently
close, to the root $t^{*}$. Performing Newton-Raphson on Eq. (22) with
starting point $\tilde{t}_{0}=0$ yields for $k\geq 0$:

$\tilde{t}_{k+1}=\tilde{t}_{k}-\frac{\displaystyle\sum_{i=0}^{3^{N}-2^{N}-1}\tilde{c}_{i}\lambda_{i}e^{\lambda_{i}\tilde{t}_{k}}}{\displaystyle\sum_{i=0}^{3^{N}-2^{N}-1}\tilde{c}_{i}\lambda_{i}^{2}e^{\lambda_{i}\tilde{t}_{k}}}.$

The Newton-Raphson method converges quadratically777Quadratic convergence
means that the number of correct decimal digits roughly doubles with each
iteration. to a local root, which is not necessarily equal to a global maximum
or minimum. Similarly, the second-order Newton-Raphson method is given by (see
Appendix A)

$\tilde{t}_{k+1}=\tilde{t}_{k}+\frac{-y^{\prime\prime}(\tilde{t}_{k})\pm\sqrt{y^{\prime\prime}(\tilde{t}_{k})^{2}-2y^{\prime}(\tilde{t}_{k})y^{\prime\prime\prime}(\tilde{t}_{k})}}{y^{\prime\prime\prime}(\tilde{t}_{k})}.$
(24)

We verify the accuracy of both approaches by comparing the peak time estimate
with the true epidemic peak time $t_{\text{peak}}$ as follows: On the complete
graph with $N=7$ nodes, we repeat 100 times: Draw uniformly distributed
$\tilde{\beta}_{kl}\in[0.05,0.25]$ and $\delta_{k}\in[0.3,0.5]$ and compute
for each set of parameters the exact peak time using Eq. (22) and the
bisection method. Figure 15a shows the time-varying prevalence $y(t)$ for each
of the 100 trials. The epidemic peak is approximately located at
$t_{\text{peak}}\approx 2$. Figure 15c shows the absolute difference between
the real peak time $t_{\text{peak}}$ and the peak time estimated by the
Newton-Raphson method for $k=0,\ldots,5$ iterations. For almost all trials,
unfortunately, the error did not converge to zero. Contrary, the error of the
second-order Newton-Raphson (2nd-NR) method in Figure 15d rapidly converges to
zero for almost888 It might happen that divergence occurs if the third
derivative $y^{\prime\prime\prime}(t)\approx 0$, in which case a higher-order
Newton-Raphson method should be implemented. all trials. Due to the cubic
convergence of 2nd-NR, only three iterations suffice to accurately approximate
the true epidemic peak time $t_{\text{peak}}$.

There are multiple reasons why the Newton-Raphson method does not converge to
the true epidemic peak time $t_{\text{peak}}$. For example, for trials that
satisfy $y^{\prime}(0)>0$ and $y^{\prime\prime}(0)>0$ and starting at
$\tilde{t}_{0}=0$, Eq. (23) tells us that $\tilde{t}_{1}<0$. Subsequent
estimates $\tilde{t}_{k}$ for $k\geq 1$ will converge to even smaller
(negative) values. For other trials, it holds that $y^{\prime\prime}(0)>0$ but
$y^{\prime\prime}(0)$ is very small. The Newton-Raphson estimate
$\tilde{t}_{1}$ is then so large, that subsequent estimates $\tilde{t}_{k}$
will diverge to $+\infty$, as $t^{*}=\infty$ is also a valid solution of
$y^{\prime}(t^{*})=0$. Additionally, the prevalence $y(t)$ can be non-
monotonic and exhibits multiple peaks. Figure 16 shows that such an effect can
already appear in small graphs (although it is not very pronounced).

Figure 16 additionally demonstrates that the mean-field approximation performs
very poorly. It is known that mean-field approximations are poor for small
contact graphs [37], but allowing for heterogeneous transition rates leads to
even worse mean-field estimates.

Figure 16: The exact solution (solid line) and mean-field approximation
(dashed line) for SIR epidemics on the complete graph with $N=7$ nodes. The
heterogeneous mean-field approximation from Eq. (30) poorly matches the exact
solution. Additionally in this example, the Markovian prevalence $y(t)$ is not
monotonic.

###  III.5 Probability of $k$ infected nodes

One of the advantages of our trinary numeral system for constructing the state
space $\mathbf{S}$ of the $3^{N}$-sized Markov chain, is that many properties
of the SIR process can be computed with ease. For example, the probability
that $k$ nodes are simultaneously infected at time $t$ characterises an
epidemic outbreak and can be used to assess the maximum impact of the disease
on the population. Given the solution $\mathbf{s}(t)$, the probability of $k$
infected nodes can be computed by projecting on the vector $\mathbf{m}$, whose
elements $m_{i}=1$ if configuration $i$ has $k$ infected nodes. In practice,
we count the number of ones in $\mathbf{x}_{i}$ to compute $m_{i}$. Figure 17
shows the probability of $k$ infected nodes in a contact graph with 7 nodes.
At time $t=0$, the probability of 1 infected node is 1, whereafter it
converges quickly to smaller values. In this example, the prevalence is rather
high, implying that many nodes can be infected simultaneously. Figure 17
illustrates the probability that all nodes are infected simultaneously is
$0.1$ at its maximum value, which is very high.

Figure 17: The probability of $k$ infected nodes in SIR epidemics on a
complete graph with $N=7$ nodes and uniformly generated
$\tilde{\beta}_{kl}\in[0.4,0.9]$ and $\delta_{k}\in[0.1,0.6]$. The probability
of $k=1$ infected node equals 1 at $t=0$, whereas all other curves start at 0.

###  III.6 Probability of group-level infections

Many stochastic and deterministic epidemic models on networks focus on the
prevalence or the probability of $k$ simultaneously infected nodes. The
benefit of our exact method is the ability to exactly determine the
probability that any group of $k$ nodes is simultaneously infected:

$y^{(k)}(t)=\frac{1}{\binom{N}{k}}\sum_{\mathcal{S}\subseteq\mathcal{N},|\mathcal{S}|=k}\textnormal{Pr}[X_{i_{1}}=1,X_{i_{2}}=1,\ldots,X_{i_{k}}=1].$

For $k=1$, we recover the prevalence $y(t)$. We know that $y^{(1)}(t)\geq
y^{(2)}(t)\geq\ldots\geq y^{(N)}(t)$, because
$\textnormal{Pr}[X=1,Y=1]\leq\textnormal{Pr}[X=1]$ for any two random
variables $X$ and $Y$. Figure 18 shows the joint infection probability of $k$
groups of nodes on a contact graph with $N=7$ nodes. All curves roughly
exhibit their peak at the same epidemic peak time, but the height of the peak
decreases with increasing $k$ as expected.

Figure 18: The joint infection probability of $k$ groups of nodes in SIR
epidemics on a complete graph with $N=7$ nodes and uniformly generated
$\tilde{\beta}_{kl}\in[0.4,0.9]$ and $\delta_{k}\in[0.1,0.6]$.

The SIR epidemic model exhibits a phase transition around the epidemic
threshold. Below the epidemic threshold, the disease dies out exponentially
fast and above the threshold, the disease persists and infects a significant
part of the population. It is conjectured that around the epidemic threshold,
the probability of $k$ groups being infected is roughly equal for several
values of $k$. The epidemic threshold $\tau_{c}$ for Markovian SIR dynamics is
not known, but lower-bounded by the mean-field threshold $\tau_{c}^{(1)}$,
which is given in Eq. (31) in Appendix B. Figure 19 shows the peak time in (a)
and height of the peak in (b). Around the epidemic threshold, which is
approximated located at the normalised effective infection rate
$x=\tau/\tau_{c}^{(1)}=1$, we do not observe that all curves are approximately
of the same order of magnitude. Most likely, a contact graph with $N=7$ nodes
is too small to exhibit a clear epidemic threshold.

(a)  (b)

Figure 19: Finding the (a) time and (b) size of the epidemic peak for $k$
groups of infected nodes in SIR epidemics on a complete graph with $N=7$ nodes
and uniformly generated $\tilde{\beta}_{kl}\in[0.4,0.9]$ and
$\delta_{k}\in[0.1,0.6]$. Around the $x=1$, which lower-bounds the epidemic
threshold, all curves do not have the same magnitude.

###  III.7 Final epidemic size

Of primary interest is the determination of the total number of infected nodes
over the course of the disease outbreak, which is also known as the final
epidemic size. The final epidemic size is equivalent to the average final
fraction of recovered nodes. We split up the $3^{N}$ states in Eq. (18) into
the set of $2^{N}$ absorbing states with eigenvalue $\lambda_{i}=0$ and the
remaining states;

$\mathbf{s}(t)=\sum_{i=0}^{2^{N}-1}c_{i}\mathbf{w}_{i}+\sum_{i=2^{N}}^{3^{N}-1}c_{i}e^{\lambda_{i}t}\mathbf{w}_{i}.$

In the limit $t\to\infty$, the second term becomes zero. Then the steady-state
distribution $\mathbf{s}_{\infty}$ is given by

$\mathbf{s}_{\infty}=\sum_{i=0}^{2^{N}-1}c_{i}\mathbf{w}_{i}.$

Denote the vector $\mathbf{m}_{R}=[0,0,1,0,0,1,1,1,2,0,...]^{T}$, whose
elements $(m_{R})_{i}$ describe the number of recovered nodes in configuration
$i$. Then the final fraction $R_{\infty}$ of recovered nodes equals

$R_{\infty}=\sum_{i=0}^{2^{N}-1}c_{i}\mathbf{m}_{R}^{T}\mathbf{w}_{i}.$ (25)

Equation (25) is a closed-form solution for the final size of the epidemic.
The mean-field approximation from Appendix B does not allow for a closed-form
solution, due to the non-linearity of the mean-field approximation. Figure 20
shows the average fraction of infected and recovered nodes over time. Although
the discrepancy between the curves looks small, the final epidemic size
$R_{\infty}\approx 0.9$, whereas the mean-field estimates
$R_{\infty}^{(1)}\approx 0.9999$, which is far off.

As a side note, we have chosen to compute the _average_ final size, but the
final size _distribution_ can be obtained in a similar manner, by constructing
a matrix $M_{R}$, whose element $(M_{R})_{i,j}=1$ if in configuration $i$
contains exactly $j$ recovered nodes.

Figure 20: The average fraction of infected nodes (blue) and average fraction
of recovered nodes (green) in a complete, weighted graph with $N=7$ nodes for
the Markovian solution (solid), mean-field solution (dashed) and Monte Carlo
simulations (asterisks). Parameters are the same as in Figure 14.

###  III.8 Exceeding hospital capacity

Hospitals have a limited capacity to treat patients simultaneously. Thus, the
number of simultaneously infected people must not exceed the hospital
capacity. If the hospital capacity is forecasted to be insufficient, counter-
measures must be taken, e.g. taking social distancing measures or scaling up
the hospital capacity. Given the disease parameters and the hospital capacity
$c_{h}$, we can compute the probability of exceeding the hospital capacity as

$p_{h}=\max_{t}\textnormal{Pr}[I(t)>c_{h}]=\max_{t}\mathbf{m}_{H}^{T}\mathbf{s}(t),$

where the element $(m_{H})_{i}$ of the $3^{N}\times 1$ vector $\mathbf{m}_{H}$
describes if configuration $i$ contains more than $c_{h}$ infected nodes and
$I(t)$ is the random variable describing the number of simultaneously infected
nodes. Using (18), we obtain

$p_{h}=\max_{t}\sum_{i=0}^{3^{N}-1}c_{i}e^{\lambda_{i}t}\mathbf{m}_{H}^{T}\mathbf{w}_{i}.$
(26)

If the probability to exceed the hospital capacity $p_{h}$ is too large,
policymakers may suggest to close facilities, shops or schools. Using
estimates of the impact of those closings on the infection rate allows for the
computation of the new exceedance probability $p_{h}$. Since the exceedance
probability is a rare event, standard Monte Carlo simulations are unsuitable
for the determination of this probability. For example, if it is required that
the hospital capacity is not exceeded with probability $p_{h}=10^{-6}$, at
least $10^{7}$ Monte Carlo simulations are required to accurately estimate the
probability, demonstrating the necessity of our exact solution.

###  III.9 The Laplace transform

In the limit of large networks, the solution (20) is composed of exponentially
many terms and cannot be computed. Fortunately, we can rewrite (20) using Abel
summation in the limit $N\to\infty$ to find (see Appendix C for the
derivation):

$y(t)=t\int_{0}^{\infty}e^{-xt}g(x)dx\hskip 28.45274pt(t>0)$ (27)

where

$g(x)=\sum_{l=0}^{[\lambda^{-1}(x)]}\tilde{c}_{l}.$ (28)

Here, $[x]$ indicates the integer part of $x$ and $\lambda^{-1}(x)$ is the
inverse eigenvalue function. The function $g(x)$ is the inverse Laplace
transform of $\frac{y(t)}{t}$ and contains all information about the
prevalence, but is expressed in the frequency $x$ domain, rather than the time
$t$ domain. For the complete graph with $N=7$ nodes, Figure 21a shows that the
contribution $g(x)$ of each frequency $x$ exhibits a complicated structure.
While (27) was derived in the limit $N\to\infty$, Figure 21b shows that the
prevalence $y(t)$ from (27) for a contact graph with only 7 nodes
approximately agrees with the exact solution (19). Equation (27) is derived
under the assumption $N\to\infty$, which is accurate everywhere, except around
$t=0$. For a finite graph, the initial condition $y(0)=1/N$ is a finite
number, whereas for $N\to\infty$, the initial prevalence converges to zero,
causing a discrepancy between the two solutions. Nevertheless, for most
networks of a realistic size, Eq. (27) seems to be an accurate approximation
of the solution (19).

(a)  (b)

Figure 21: (a) The function $g(x)$ from Eq. (28) and (b) the prevalence (27)
of the SIR process on a contact graph with $N=7$ nodes and uniformly generated
$\tilde{\beta}_{kl}\in[0.4,0.9]$ and $\delta_{k}\in[0.1,0.6]$. Additionally,
(b) shows the high-amplitude filter with cut-off amplitude $\epsilon=1$ and
low-pass filter with frequency $x\in[0,2]$. The low-pass filter uses
approximately 330 out of 2059 non-zero eigenvalues.

Instead of the whole function $g(x)$, partial information on $g(x)$ may be
sufficient to reconstruct the time-varying prevalence $y(t)$. We apply two
techniques from signal processing to $g(x)$. Method 1 performs a high-
amplitude filter, keeping only frequencies with amplitude larger than some
bound $\epsilon$ and all other contributions are set to zero;

$g_{1}(x)=\begin{cases}g(x),&\qquad g(x)\geq\epsilon\\\
0,&\qquad\text{otherwise}\end{cases}$

Second, we employ a low-pass filter, only keeping the contributions of
frequencies in the interval $[0,f_{0}]$;

$g_{2}(x)=\begin{cases}g(x),&\qquad x\leq f_{0}\\\
0,&\qquad\text{otherwise}\end{cases}$

The main reason behind $g_{2}(x)$ is that small (in modulus) frequencies
correspond to small (in modulus) eigenvalues, which have the largest
contributions in the Laplace transform (27). Figure 21b shows the result of
applying both filters to $g(x)$. The low-pass filter is superior to the high-
amplitude filter, because small values of $x$ have a more significant
contribution in the Laplace transform (27). Even for an amplitude cut-off
$\epsilon=1$ in Figure 21b, which is relatively small, the low-pass filter
better approximates the true solution.

Interestingly, eigenmode truncation (as explained in Section II.3) is only
effective for large times, as shown in Figure 22. Eigenmode truncation also
maintains the smallest (in modulus) eigenvalues, but apparently performs much
worse compared to the low-pass filter of the Laplace transform. The primary
reason for this discrepancy is that eigenmode truncation considers finite
networks whereas the Laplace transform was derived under the assumption of an
infinitely large network. Even for a network with $N=7$ nodes, the Laplace
transform (27) performs better using the same number of eigenmodes. The
advantage of the latter is that the solution can no longer explode around
$t=0$ because the $t$-term in (27) guarantees $y(0)=0$.

Figure 22: Eigenmode truncation in SIR epidemics on the complete graph with
$N=7$ nodes and uniformly generated $\tilde{\beta}_{kl}\in[0.4,0.9]$ and
$\delta_{k}\in[0.1,0.6]$. Approximately 330 out of 2059 non-zero eigenvalues
are used.

##  IV Conclusion

We investigated continuous-time Markovian compartmental models on networks. We
focussed on transient epidemic compartmental models such as SI and SIR, where
re-infections do not appear. Then the compartmental graph describing the
transitions between compartments does not contain directed cycles. We showed
that such a continuous-time Markovian compartmental model admits an analytic
solution using eigendecomposition of the infinitesimal generator. We
demonstrated our method on the SI and SIR processes. We derived the
eigenvalues of the infinitesimal generators and showed that the eigenvalues
are related to S-I cut sets in the contact graph. Assuming that all non-zero
eigenvalues are distinct, we constructed the exact solution of the SI and SIR
process on heterogeneous networks. Additionally, we showed that the transient
compartmental model can be extended by adding self-infections, simplicial
contagion, temporal networks and non-Markovian dynamics while maintaining an
exact, analytic solution for the time-varying infection probabilities of all
nodes.

The advantage of our exact approach is twofold. First, the analytic solution
on small networks can serve as a calibration tool for (non)-Markovian
compartmental simulators on networks and subsequently allows researchers to
determine the number of simulation events to achieve a certain accuracy.
Second, we showed that the underlying state space scales exponentially with
the number of nodes. Much effort is devoted to simplify the state space, e.g.
by aggregating states in the Markov graph. Our exact time-dependent solution
can assess the quality of the proposed aggregation methods in an exact manner.

Finally, the advance in quantum computers and quantum algorithms may help to
exactly determine the time-varying prevalence for much larger networks. There
is an obvious analogy between the viral state of a node in a network and a
quantum state of a particle in a quantum device. Also, the entire state space
of a quantum system is increasing exponentially in the number of qubits.
Hence, we expect a computational breakthrough and an unraveling of the SIR
epidemic phase transition with the further development of quantum computers.

Acknowledgements We thank Robin Persoons and Brian Chang for useful
discussions on the infinitesimal generator for the GEMF framework.

This research has been funded by the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No 101019718).

##  Appendix A Newton’s root-finding method

We shortly explain Newton’s method of finding the root $t^{*}$ of a function
$f(t)$. In our case, we are interested in the root of $y^{\prime}(t)$ at
$t_{\text{peak}}$, but we write $f(t)$ and $t^{*}$ for generality.

###  A.1 First-order Newton-Raphson

The Taylor expansion of the function $f(t^{*})$ around $t^{*}=t$ equals

$\begin{split}f(t^{*})&=f_{0}(t)+(t^{*}-t)f_{1}(t)\\\
&\quad+(t^{*}-t)f_{2}(t)+\mathcal{O}((t^{*}-t)^{3})\end{split}$ (29)

where we defined

$f_{k}(t)=\frac{f^{(k)}(t)}{k!}.$

Since $t^{*}$ is a zero of $f(t)$, we find $f(t^{*})=0$. Using only terms up
to first order in (29), we find

$0=f_{0}(t)+(t^{*}-t)f_{1}(t)+\mathcal{O}((t^{*}-t)^{2})$

which can be rewritten as

$t^{*}=t-\frac{f_{0}(t)}{f_{1}(t)}+\mathcal{O}((t-t^{*})^{2})$

Taking $t^{*}=\tilde{t}_{k+1},t=\tilde{t}_{k}$ and neglecting the second-order
terms, we find the Newton-Raphson iterative scheme:

$\tilde{t}_{k+1}=\tilde{t}_{k}-\frac{f_{0}(\tilde{t}_{k})}{f_{1}(\tilde{t}_{k})}.$

###  A.2 Second-order Newton-Raphson

Instead of using only first-order terms, we can also use second-order terms in
(29). Again using $f(t^{*})=0$, we find

$0=f(t)+(t^{*}-t)f_{1}(t)+(t^{*}-t)^{2}f_{2}(t)+\mathcal{O}((t^{*}-t)^{3})$

Neglecting the third-order terms, we can solve for $t^{*}$;

$t^{*}=t+\frac{-f_{1}(t)\pm\sqrt{f_{1}(t)^{2}-4f_{0}(t)f_{2}(t)}}{2f_{2}(t)}$

Using $t^{*}=\tilde{t}_{k+1}$ and $t=\tilde{t}_{k}$, we find the second-order
Newton-Raphson iterative scheme:

$\tilde{t}_{k+1}=\tilde{t}_{k}+\frac{-f_{1}(\tilde{t}_{k})\pm\sqrt{f_{1}(\tilde{t}_{k})^{2}-4f_{0}(\tilde{t}_{k})f_{2}(\tilde{t}_{k})}}{2f_{2}(\tilde{t}_{k})}$

Substituting $f(t)=y^{\prime}(t)$ gives Eq. (24). The major advantage of the
second-order Newton-Raphson method is that the radius of convergence seems
larger than Newton-Raphson and the speed of convergence is faster (cubic
versus quadratic convergence).

##  Appendix B Mean-field equations

The heterogeneous mean-field approximation of the stochastic SI(S) process is
known as the N-Intertwined Mean-Field Approximation (NIMFA) [26] and equals

$\frac{d\,v_{k}}{dt}=-\delta_{k}v_{k}+(1-v_{k})\sum_{j=1}^{N}\tilde{\beta}_{kl}v_{l}$

where $v_{k}(t)$ is the probability that node $k$ is infected at time $t$ and
$\tilde{\beta}_{kl}=\beta_{kl}a_{kl}$.

Similarly, the heterogeneous individual mean-field approximation for the SIR
process is given by [38]

$\displaystyle\begin{split}\frac{d\,s_{k}}{dt}&=-\sum_{l=1}^{N}\tilde{\beta}_{kl}s_{k}v_{l}\\\
\frac{d\,v_{k}}{dt}&=\sum_{l=1}^{N}\tilde{\beta}_{kl}s_{k}v_{l}-\delta_{k}v_{k}\end{split}$
(30)

where $s_{k}(t)$ and $v_{k}(t)$ represent the probabilities that node $k$ is
susceptible or infected at time $t$, respectively.

The mean-field SIS and SIR threshold is [39]

$\tau_{c}^{(1)}=\lambda_{\max}(S^{-1}B)$ (31)

where $S=\text{diag}(\delta_{1},\ldots,\delta_{N})$ is the $N\times N$ curing
rate matrix and $B$ is the $N\times N$ infection rate matrix, whose elements
are $\tilde{\beta}_{kl}$.

##  Appendix C Abel summation

A finite-$n$-analysis on the sum $\sum_{k=m}^{n}a_{k}\,b_{k}$, where $n\geq
m$, illustrates interesting formal manipulations dealing with partial sums,
$s_{k}=\sum_{l=m}^{k}a_{l}$ where $k\geq m$, attributed to Niels Abel.

The basic observation is that $a_{k}=s_{k}-s_{k-1}$ for $k>m$. Thus,

$\displaystyle\sum_{k=m}^{n}a_{k}\,b_{k}$
$\displaystyle=a_{m}\,b_{m}+\sum_{k=m+1}^{n}(s_{k}-s_{k-1})\,b_{k}$
$\displaystyle=a_{m}\,b_{m}+\sum_{k=m+1}^{n}s_{k}\,b_{k}-\sum_{k=m}^{n-1}s_{k}\,b_{k+1}$
$\displaystyle=a_{m}\,b_{m}+\sum_{k=m+1}^{n-1}s_{k}\,(b_{k}-b_{k+1})$
$\displaystyle\qquad+s_{n}\,b_{n}-s_{m}\,b_{m+1}$

Since $s_{m}=a_{m}$, we arrive at Abel’s partial summation valid for any
integer $m\leq n$,

$\sum_{k=m}^{n}a_{k}\,b_{k}=\sum_{k=m}^{n-1}\left(\sum_{l=m}^{k}a_{l}\right)\,(b_{k}-b_{k+1})+b_{n}\left(\sum_{l=m}^{n}a_{l}\right)$
(32)

A reversal of the $k$\- and $l$-sum again returns to the sum at left-hand
side.

As an application of Abel summation, we consider the sum
$h_{n}\left(t\right)=\sum_{k=0}^{n}a_{k}\,e^{-\lambda_{k}t}$, where
$0\leq\operatorname{Re}\left(\lambda_{0}\right)\leq\operatorname{Re}\left(\lambda_{1}\right)\leq\cdots\leq\operatorname{Re}\left(\lambda_{n}\right)$,
that can appear as a solution of $n+1$-th order linear differential equation.
Abel summation (32) yields

$h_{n}\left(t\right)=\sum_{k=0}^{n-1}\left(\sum_{l=0}^{k}a_{l}\right)\,(e^{-\lambda_{k}t}-e^{-\lambda_{k+1}t})+e^{-\lambda_{n}t}\left(\sum_{l=0}^{n}a_{l}\right)$

Invoking
$-t\int_{\lambda_{k+1}}^{\lambda_{k}}e^{-xt}dx=e^{-\lambda_{k}t}-e^{-\lambda_{k+1}t}$,
we obtain

$h_{n}\left(t\right)=t\sum_{k=0}^{n-1}\int_{\lambda_{k}}^{\lambda_{k+1}}\left(\sum_{l=0}^{k}a_{l}\right)e^{-xt}dx+e^{-\lambda_{n}t}\left(\sum_{l=0}^{n}a_{l}\right)$

Let the function $y=\lambda\left(x\right)$ define the sequence
$\left\\{\lambda_{k}\right\\}_{0\leq k\leq n}$ by
$\lambda\left(k\right)=\lambda_{k}$ at integer values of $x$. The inverse
function $x=\lambda^{-1}\left(y\right)$ maps the index to
$k=\lambda^{-1}\left(\lambda_{k}\right)$. Since
$k\leq\lambda^{-1}\left(x\right)<k+1$ for $\lambda_{k}\leq x<\lambda_{k+1}$,
the integer part $\left[.\right]$ of
$\left[\lambda^{-1}\left(x\right)\right]=k$ for $\lambda_{k}\leq
x<\lambda_{k+1}$ and we have that

$\displaystyle\sum_{k=0}^{n-1}\int_{\lambda_{k}}^{\lambda_{k+1}}\left(\sum_{l=0}^{k}a_{l}\right)e^{-xt}dx$
$\displaystyle=\sum_{k=0}^{n-1}\int_{\lambda_{k}}^{\lambda_{k+1}}e^{-xt}\sum_{l=0}^{\left[\lambda^{-1}\left(x\right)\right]}a_{l}\;dx$
$\displaystyle=\int_{\lambda_{0}}^{\lambda_{n}}e^{-xt}\sum_{l=0}^{\left[\lambda^{-1}\left(x\right)\right]}a_{l}\;dx$

which leads to

$h_{n}\left(t\right)=t\int_{\lambda_{0}}^{\lambda_{n}}e^{-xt}\sum_{l=0}^{\left[\lambda^{-1}\left(x\right)\right]}a_{l}\;dx+e^{-\lambda_{n}t}\left(\sum_{l=0}^{n}a_{l}\right)$
(33)

Relation (33) shows for $t=0$ that $h_{n}\left(0\right)=\sum_{l=0}^{n}a_{l}$
for all $n$ (also when $n\rightarrow\infty$).

However, if $\lambda_{n}\rightarrow\infty$ for $n\rightarrow\infty$ and
$\lambda_{0}=0$, then the last sum in (33) vanishes for all positive $t>0$ and

$\frac{h_{\infty}\left(t\right)}{t}=\frac{1}{t}\sum_{k=0}^{\infty}a_{k}\,e^{-\lambda_{k}t}=\int_{0}^{\infty}e^{-xt}g\left(x\right)dx\hskip
28.45274pt\left(t>0\right)$ (34)

is the Laplace transform of the sum
$g\left(x\right)=\sum_{l=0}^{\left[\lambda^{-1}\left(x\right)\right]}a_{l}$.
When $t$ tends to zero as $t=\frac{c}{\lambda_{n}}$ in the limit
$n\rightarrow\infty$, then $e^{-\lambda_{n}t}=e^{-c}$, which is a constant,
and the first term in (33) vanishes, but the second term remains and
$e^{-c}=1$ so that $h_{\infty}\left(0\right)=\sum_{l=0}^{\infty}a_{l}$, which
is a prerequisite for continuity for $t\geq 0$ as well. Hence, the point $t=0$
needs care999If $h_{\infty}\left(0\right)=0$, then the limit
$\lim_{t\rightarrow
0}\frac{h_{\infty}\left(t\right)}{t}=h_{\infty}^{\prime}\left(0\right)$ and
(34) shows that
$\int_{0}^{\infty}g\left(x\right)dx=h_{\infty}^{\prime}\left(0\right)$..

Observe from
$g\left(x\right)=\sum_{l=0}^{\left[\lambda^{-1}\left(x\right)\right]}a_{l}$ in
the limit $n\rightarrow\infty$ where $\lambda_{n}\rightarrow\infty$, that
$g\left(0\right)=a_{0}=\lim_{t\rightarrow\infty}h_{\infty}\left(t\right)$,
while
$\lim_{x\rightarrow\infty}g\left(x\right)=\sum_{l=0}^{\infty}a_{l}=h_{\infty}\left(0\right)$.
The inverse Laplace transform is

$g\left(x\right)=\sum_{l=0}^{\left[\lambda^{-1}\left(x\right)\right]}a_{l}=\frac{1}{2\pi
i}\int_{c-i\infty}^{c+i\infty}\frac{h_{\infty}\left(t\right)}{t}e^{xt}dt\hskip
28.45274pt\left(c>0\right)$

The infinitesimal generator $Q$ is minus a weighted Laplacian matrix and the
smallest eigenvalue of any Laplacian matrix is zero. Hence, $\lambda_{0}=0$ is
satisfied, but $\lambda_{n}$ is always finite for any finite graph with $n$
nodes. Only if the graph is sufficiently large, then the Laplace transform
(34) is exact for $t>0$, while $h_{n}\left(0\right)=\sum_{l=0}^{n}a_{l}$ for
all $n$.

##  Appendix D General compartments

The Generalised Epidemic Mean-Field (GEMF) framework describes general
compartmental epidemic models on networks [2]. Each node $k$ at time $t$ is in
one of the $c$ possible compartments. The compartmental graph describes all
possible transitions between the different compartments and additionally
specifies whether the transition is a link-based or a node-based transition.
Figure 23 shows several examples of well-known compartmental graphs in
epidemiology. As shown in Figure 23, multi-links can exist in the
compartmental graph if multiple processes cause the same transition. For
example, infections may happen based on contact between infected and
susceptible nodes, but susceptible nodes can also be infected indirectly,
which is modelled as a self-infection process. Both the infection and self-
infection process (in $\varepsilon$-SIR) change the viral state of a node from
susceptible to infected, thus allowing for multi-links in the compartmental
graph.

Figure 23: Several well-known compartmental epidemic models. Left: SIS and
SIR. Middle: SEIR and SIRS. Right: SEIARDS. Absorbing states are shown as grey
nodes. Full links described link-based transitions (and the letter(s) on the
link describe the compartment(s) that triggers the transition) whereas dashed
links indicate node-based transitions. Here, the exposed component E includes
nodes that have already been infected but are not yet infectious. Furthermore,
A represents asymptomatic nodes that do not show any symptoms but can still
infect susceptible nodes (although at a lower rate) and D represents the
deceased nodes.

We use a base-$c$ numeral system to denote each of the $c$ possible viral
states of each individual. We denote $x_{k}(i)$ as the viral state of node $k$
in configuration $i$, which can be any integer between $0$ and $c-1$. The
viral state vector $\mathbf{x}=(x_{N},x_{N-1},\ldots,x_{1})^{T}$ represents
the viral state of all nodes. Given a particular viral state vector
$\mathbf{x}$, we compute the configuration number $i$ (which is a decimal
number) based on the base-$c$ numerals as

$i=\sum_{k=1}^{N}x_{k}(i)\cdot c^{k-1}$ (35)

Representation (35) guarantees that each viral state vector $\mathbf{x}$
receives a unique configuration number $i$. The total number of configurations
equals $c^{N}$, because each node has $c$ possible viral states and there are
$N$ nodes in the contact graph.

The transitions in the Markov chain with $c^{N}$ states are described by the
$c^{N}\times c^{N}$-dimensional infinitesimal generator $Q$. We assume that
the transitions can be subdivided into a set $\Theta$ of link-based
transitions and a set $\Xi$ of node-based transitions. Each link-based
transition $\theta$ describes the transition of node $l$ to change from
compartment $\theta_{1}$ to $\theta_{2}$ with rate $\beta_{\theta,kl}a_{kl}$
under the influence of an adjacent node $k$ whose current state is
$\theta_{3}$. Additionally, a node-based transition $\xi$ changes node $l$
from state $\xi_{1}$ to $\xi_{2}$ with rate $\delta_{\xi,l}$.

Due to the Markov property, any transition in the Markov graph will change the
viral state of exactly one node. Thus the viral state vector $\mathbf{x}$
before and after the transition will differ by exactly one bit. The _Hamming
distance_ $d_{H}(\mathbf{x},\mathbf{y})$ between two vectors $\mathbf{x}$ and
$\mathbf{y}$ is the number of elements that is different in $\mathbf{x}$ and
$\mathbf{y}$. Transitions in the Markov graph occur if a single node changes
its state, which corresponds to states $\mathbf{x}$ and $\mathbf{y}$ having
Hamming distance $d_{H}(\mathbf{x},\mathbf{y})=1$. Note that the reverse does
not hold, i.e. having $d_{H}(\mathbf{x},\mathbf{y})=1$ does not necessarily
imply that a transition is possible, because transitions are dependent on the
existence of directed links in the compartmental graph (recall Figure 23).

For general Markovian compartmental models with $c$ compartments, the set
$\Theta$ of link-based transitions and the set $\Xi$ of node-based
transitions, the infinitesimal generator $Q$ can be constructed as follows:

$\displaystyle\begin{split}q_{ij}&=\sum_{\theta\in\Theta}\sum_{k=1}^{N}\beta_{\theta,kl}a_{kl}\mathbf{1}_{\\{x_{l}(i)=\theta_{1}\cap\,x_{l}(j)=\theta_{2}\cap\,x_{k}(i)=\theta_{3}\\}}\\\
&\qquad\qquad+\sum_{\xi\in\Xi}\delta_{\xi,l}\mathbf{1}_{\\{x_{l}(i)=\xi_{1}\cap\,x_{l}(j)=\xi_{2}\\}}\\\
&\qquad\text{if }d_{H}(\mathbf{x}(i),\mathbf{x}(j))=1\text{ and
}l=\lfloor\log_{c}(|i-j|)\rfloor+1\\\ \end{split}$ (36a) $\displaystyle
q_{ii}$ $\displaystyle=-\sum_{\begin{subarray}{c}j=0\\\ j\neq
i\end{subarray}}^{c^{N}-1}q_{ij}$ (36b)

where $\log_{c}(\ldots)$ is the base-$c$ logarithm and $\lfloor x\rfloor$ is
$x$ rounded down to the nearest integer. Transitions from configuration $i$ to
$j$ only occur by changing the viral state of exactly one node, thus exactly
one bit in the viral state vector $\mathbf{x}$. For a given configuration $i$,
the maximum possible number of transitions is given by the maximum out-degree
$d_{\max}$ of the compartmental graph (recall Figure 23), because the out-
degree of a viral state specifies the number of possible next viral states.
For a given configuration $i$, the number of possible transitions is given by
$d_{\max}N$, because each node can change its viral state. Taking into account
the possibility of not making any change, the number of non-zero elements in
the infinitesimal generator $Q$ is upper bounded by $(d_{\max}N+1)c^{N}$,
which is significantly less than the dense $c^{N}\times c^{N}$ matrix.

Using the infinitesimal generator $Q$ for any GEMF epidemic model, one can
perform the same analysis as in the main body in this paper, i.e. analysing
the eigenvalues and eigenvectors and computing the exact time-dependent
solution.

##  Appendix E Paths in the SIR Markov graph

Each state in the SIR Markov graph has maximally $N$ out-going links, i.e. the
maximum out-degree is $N$. The reason is that every state describes the viral
state of $N$ nodes. During the Markovian SIR process, only one event can occur
simultaneously, which results in one node changing its state. Since the graph
consists of $N$ nodes, at most $N$ states can be reached from the starting
state. Similarly, the maximum in-degree equals $N$.

We are interested in the _number of paths between state $i$ and $j$_. The
maximum number of paths is attained if (i) state $i$ corresponds to the
initial state, in which one node is infected and the other nodes are
susceptible, (ii) state $j$ corresponds to the all-recovered state and (iii)
the contact graph $G$ equals the complete graph. Then we obtain the following
result:

###### Lemma 2

The maximum number of paths $|\mathcal{P}|$ in the SIR Markov graph for the
complete graph with $N\geq 3$ is bounded by

$(N-1)!\cdot
N!<(2N-1)((N-1)!)^{2}<|\mathcal{P}|<\frac{(2N-1)!}{2^{N-1}}<\frac{(2N)!}{2^{N}}<N^{2N-1}.$
(37)

Proof. Upper bound 1: The amount of paths between state $i$ and $j$ is
maximally equal to $N^{2N-1}$, because the maximal hopcount equals $2N-1$ and
the maximum out-degree equals $N$. The upper bound is never attained, because
many states have out-degree smaller than $N$.

Upper bound 2: Consider adding a positive self-infection rate $\varepsilon>0$
to the SIR process. Then the Markov graph simplifies to only one absorbing
state: the all-recovered state. Not starting from the initial state with a
single infected node, but from the all-healthy state, we definitely have an
upper bound on the number of paths in the SIR process. The number of paths in
this case is equivalent to asking in how many ways the $N$ nodes can change
from S to I, and from I to R. In general, there are $(2N)!$ ways. However, for
a given node, it is required that the node gets infected before the node
recovers. In exactly half of the cases, the ordering is wrong, thus we obtain
$(2N)!/2$ paths. Repeating this argument for all $N$ nodes leads to
$(2N)!/(2^{N})$ possible paths.

Upper bound 3: In our case, the process is initiated with a single infected
node. Using the argument from upper bound 2 for the $N-1$ initially non-
infected nodes, we find $(2N-2)!/(2^{N-1})$ paths. However, besides the
transition from the $N-1$ initially susceptible nodes from $S\to I\to R$, the
initially infected node must also recover. There are exactly $2N-1$ ways to
place the recovery of the initial node in the sequence of the total $2N-1$
events. This leads to $(2N-1)!/(2^{N-1})$ paths. This number of paths is still
an upper bound, because we allow for transitions from absorbing states to non-
absorbing states, which is not possible in the SIR process.

Lower bound 1: Computing all paths is complicated, so we restrict ourselves
here to a subset of all possible paths. We only consider paths that pass
through the all-infected state. Starting from a single infected node, the
other $N-1$ nodes must become infected and the only thing that matters is the
order in which the nodes are infected. This equals to $(N-1)!$ possible paths.
From the all-infected state to the all-recovered state, all nodes must
recover, and the only thing that matters is the order in which the nodes
recover. This yields another $N!$ possible paths. In total, this equals
$(N-1)!\cdot N!$ paths that pass through the all-infected state. Naturally,
there are other paths too, causing the lower bound to be non-tight.

Lower bound 2: The other paths might pass through a state with all nodes
infected, except one node that remains susceptible and another is already
recovered (this only happens if $N\geq 3$). There are in total $(N-1)^{2}$ of
such paths, because the susceptible node can be chosen among all that are not
initially susceptible and the recovered node can be chosen among all nodes not
remaining susceptible. We split up the computation in two parts:

Part 1: from initial state to $(111111102)$ state.
In total $N-2$ nodes must be infected and 1 node must recover. To prevent
issues, we assume that the recovery event is the last event. This is
definitely a lower bound. Then there are $(N-2)!$ possible paths.

Part 2: from $(1111111102)$ state to all-recovered state $(2222222222)$.
In total $1$ node must be infected and $N-1$ nodes must recover. To prevent
issues, we assume that the infection event is the first event. Again, this is
definitely a lower bound. Then there are $(N-1)!$ possible paths.

In total, we conclude that the number of paths through state of the form
$(11111111102)$ is lower-bounded by $(N-1)((N-1)!)^{2}$. Hence, the lower
bound can be improved to $(2N-1)((N-1)!)^{2}$ for $N\geq 3$. $\square$

###### Corollary 3

For any connected contact graph $G$, the maximum number of paths
$|\mathcal{P}|$ in the SIR Markov graph is bounded by

$N!\leq|\mathcal{P}|<\frac{(2N-1)!}{2^{N-1}}$

Proof. The upper bound follows from Lemma 2. The lower bound follows the same
reasoning as the first lower bound in Lemma 2: we only consider the subset of
paths that pass through the all-infected state. The number of paths from the
initial state to the all-infected state is at least 1 in a connected graph,
but the number of paths from the all-infected state to the all-recovered state
equals $N!$, as the recovery of each node may happen independently of each
other. $\square$

Corollary 3 implies that the maximum number of paths grows more than
exponentially in the number of nodes $N$. For the complete graph, Table 1
shows a numerical evaluation of the number of paths $|\mathcal{P}|$, which is
computed by enumerating all paths in the $3^{N}$-sized Markov graph. The
combinatorial explosion of the number of paths prevents the computation of the
exact SIR solution on any network of realistic size.

Table 1: The maximal number of paths in the SIR Markov graph in the complete graph between the initial state with a single infected node and all other nodes susceptible, towards the all-recovered state on a complete contact graph $G$. The lower bounds (LB) and upper bounds (UB) are from Eq. (37). $N$ | LB1 | LB2 | $|\mathcal{P}|$ | UB1 | UB2 | UB3
---|---|---|---|---|---|---
2 | 2 | - | 2 | 3 | 6 | 8
3 | 12 | 20 | 20 | 30 | 90 | 243
4 | 144 | 252 | 444 | 630 | 2,520 | 16,384
5 | 2,880 | 5,184 | 16,944 | 22,680 | 113,400 | 1,952,125
6 | 86,400 | 158,400 | 979,440 | 1,247,400 | 7,484,400 | $3.6\cdot 10^{8}$

###  E.1 Recursive formula for the complete graph

We compute the maximum number of paths in the complete graph exactly using a
recursive formula. As before, the initial state is the state with one infected
node and the other nodes are susceptible. The final state corresponds to the
all-recovered state. Let $r_{kl}$ be the number of paths, starting from the
initial state to a state that has $k$ infected nodes and $l$ recovered nodes,
which is denoted by $(k,l)$. The number of paths $r_{kl}$ can be computed
iteratively, by looking at the previous transition, which is either related to
the infection or recovery of a node. Starting with $k-1$ infected nodes and
$l$ recovered nodes, one of the susceptible nodes can get infected. Since we
assume that the contact graph is the complete graph, all of the infected nodes
can infect the susceptible nodes, leading to $N-k-l+1$ possible paths from
state $(k-1,l)$ to $(k,l)$. Similarly, starting in the state with $k+1$
infected nodes and $l-1$ recovered nodes, we make a transition to $(k,l)$ if
one of the infected nodes recovers. Since there are $k+1$ infected nodes, each
of those nodes can recover, leading to $k+1$ possible paths from $(k+1,l-1)$
to $(k,l)$. Then the following recursion for $r_{kl}$ is obtained:

$\displaystyle r_{kl}$ $\displaystyle=$
$\displaystyle(N-k-l+1)\,r_{k-1,l}+(k+1)\,r_{k+1,l-1}$ (38) $\displaystyle
r_{10}$ $\displaystyle=$ $\displaystyle 1.\qquad\qquad\qquad\textnormal{by
definition}$

Since the process is initiated at $k=1$ infected node and $l=0$ recovered
nodes, we define $r_{10}=1$. The number of paths from the initial state to the
all-recovered state equals $r_{0N}$.

Figure 24: The simplified SIR Markov graph on a complete contact graph with
$N$ nodes, here shown with $N=3$. The node labels $kl$ denote the number of
infected nodes $k$ and the number of recovered nodes $l$.

Numerical evaluation of the recursion (38) for $r_{kl}$ indicates that the
number of paths equals

$|\mathcal{P}|=(N-1)!\cdot f(N)$

where $f(N)$ is given by https://oeis.org/A000698. Unfortunately, $f(N)$ does
not seem to possess a closed-form solution. Also, the generating function
$R(x,y)=\sum_{k=0}^{\infty}\sum_{l=0}^{\infty}r_{kl}x^{k}y^{l}$, applied to
the recursion (38) and resulting in the partial differential equation

$(x^{2}-y)\dfrac{\partial R}{\partial x}+xy\dfrac{\partial R}{\partial
y}-(Nx-1)\,R(x,y)=0$

does not seem to admit a closed-form solution.

(a)  Compartmental graph (b)  Markov graph

Figure 25: The spread of SIR epidemics on a network with two connected nodes.
(a) the SIR compartmental graph and (b) the SIR Markov graph with $3^{2}$
states.

##  Appendix F SIR exact solution for small graphs

The SIR equations can be solved explicitly for small graphs. Here, we provide
those exact solutions on the complete and path graph with $N=2$ and $N=3$
nodes.

###  F.1 Solution on $N=2$ nodes

For $N=2$ nodes, the Markov graph contains 9 states and is shown in Figure 25.
Starting in the initial state $i=1$, that is, node 1 being infected and node 2
susceptible, the exact solution $\mathbf{s}(t)$ equals:

$\displaystyle s_{0}(t)$ $\displaystyle=0$ $\displaystyle s_{1}(t)$
$\displaystyle=e^{-(\beta_{21}+\delta_{1})t}$ $\displaystyle s_{2}(t)$
$\displaystyle=\frac{\delta_{1}}{\beta_{21}+\delta_{1}}\left(1-e^{-(\beta_{21}+\delta_{1})t}\right)$
$\displaystyle s_{3}(t)$ $\displaystyle=0$ $\displaystyle s_{4}(t)$
$\displaystyle=\frac{\beta_{21}}{\beta_{21}-\delta_{2}}\left(e^{-(\delta_{1}+\delta_{2})t}-e^{-(\beta_{21}+\delta_{1})t}\right)$
$\displaystyle s_{5}(t)$
$\displaystyle=\frac{\beta_{21}}{\beta_{21}+\delta_{1}-\delta_{2}}e^{-\delta_{2}t}+\frac{\beta_{21}}{\delta_{2}-\beta_{21}}e^{-(\delta_{1}+\delta_{2})t}+\frac{\beta_{21}\delta_{1}}{(\beta_{21}-\delta_{2})(\beta_{21}+\delta_{1}-\delta_{2})}e^{-(\beta_{21}+\delta_{1})t}$
$\displaystyle s_{6}(t)$ $\displaystyle=0$ $\displaystyle s_{7}(t)$
$\displaystyle=e^{-\delta_{1}t}+\frac{\beta_{21}}{\delta_{2}-\beta_{21}}e^{-(\delta_{1}+\delta_{2})t}+\frac{\delta_{2}}{\beta_{21}-\delta_{2}}e^{-(\beta_{21}+\delta_{1})t}$
$\displaystyle s_{8}(t)$
$\displaystyle=\frac{\beta_{21}}{\beta_{21}+\delta_{1}}-e^{-\delta_{1}t}-\frac{\beta_{21}}{\beta_{21}+\delta_{1}-\delta_{2}}e^{-\delta_{2}t}+\frac{\beta_{21}}{\beta_{21}-\delta_{2}}e^{-(\delta_{1}+\delta_{2})t}$
$\displaystyle\qquad+\left(\frac{\delta_{1}}{\beta_{21}+\delta_{1}}-\frac{\delta_{2}}{\beta_{21}-\delta_{2}}+\frac{\delta_{2}-\delta_{1}}{\beta_{21}+\delta_{1}-\delta_{2}}\right)e^{-(\beta_{21}+\delta_{1})t}$

The prevalence follows by a weighted sum over the states $\mathbf{s}(t)$,
where the weight $w_{i}$ equals the number of infected nodes in state $i$.
Then we find

$y(t)=\frac{1}{2}\left(e^{-\delta_{1}t}+\frac{\beta_{21}}{\beta_{21}+\delta_{1}-\delta_{2}}e^{-\delta_{2}t}-\frac{\beta_{21}}{\beta_{21}+\delta_{1}-\delta_{2}}e^{-(\beta_{21}+\delta_{1})t}\right)$
(39)

and the average fraction of recovered nodes $r(t)$ equals

$\displaystyle r(t)$
$\displaystyle=\frac{1}{2}\left(1+\frac{\beta_{21}}{\beta_{21}+\delta_{1}}-e^{-\delta_{1}t}-\frac{\beta_{21}}{\beta_{21}+\delta_{1}-\delta_{2}}e^{-\delta_{2}t}+\frac{\beta_{21}\delta_{2}}{(\beta_{21}+\delta_{1})(\beta_{21}+\delta_{1}-\delta_{2})}e^{-(\beta_{21}+\delta_{1})t}\right)$
(40)

Figure 26 shows the solution of the exact SIR model for two cases: (left)
below the epidemic threshold and (right) above the epidemic threshold. Below
the threshold, the probability to arrive at the state $s_{8}$ is at most 0.33
whereas it is over 0.6 above the threshold.

(a)  (b)

(c)  (d)

Figure 26: Simulation of the SIR model on a graph with $N=2$ nodes. The
recovery rates are $\delta_{1}=2.1$ and $\delta_{2}=2.2$. It shows (a,b) the
probability to be in a particular state, (c,d) the prevalence $y(t)$ and the
average fraction of recovered nodes $r(t)$ for (a,c) infection rate
$\beta_{21}=1$ and (b,d) $\beta_{21}=4$.

###  F.2 Solution on a complete graph with $N=3$ nodes

For $N=3$ nodes, the contact graph is a triangle. For simplicity, we assume
that the contact graph is symmetric: $\beta_{23}=\beta_{32}$. Its Markov graph
contains 27 states and is visualised in Figure 11 in the main text. Starting
in the initial state, that is, node 1 being infected and node 2 and 3 being
susceptible, the exact solution equals

$\displaystyle s_{0}(t)$ $\displaystyle=0$ $\displaystyle s_{1}(t)$
$\displaystyle=e^{-(\beta_{21}+\beta_{31}+\delta_{1})t}$ $\displaystyle
s_{2}(t)$
$\displaystyle=\frac{\delta_{1}}{\beta_{21}+\beta_{31}+\delta_{1}}\left(1-e^{-(\beta_{21}+\beta_{31}+\delta_{1})t}\right)$
$\displaystyle s_{3}(t)$ $\displaystyle=0$ $\displaystyle s_{4}(t)$
$\displaystyle=\frac{\beta_{21}}{\beta_{32}-\beta_{21}+\delta_{2}}\left(e^{-(\beta_{21}+\beta_{31}+\delta_{1})t}-e^{-(\beta_{31}+\beta_{32}+\delta_{1}+\delta_{2})t}\right)$
$\displaystyle s_{5}(t)$
$\displaystyle=\frac{\beta_{21}\delta_{1}}{(\beta_{31}+\delta_{1})(\beta_{21}+\beta_{31}-\beta_{32}+\delta_{1}-\delta_{2})}e^{-(\beta_{32}+\delta_{2})t}$
$\displaystyle\qquad+\frac{\beta_{21}\delta_{1}}{(\beta_{21}-\beta_{32}-\delta_{2})(\beta_{21}+\beta_{31}-\beta_{32}+\delta_{1}-\delta_{2})}e^{-(\beta_{21}+\beta_{31}+\delta_{1})t}$
$\displaystyle\qquad+\frac{\beta_{21}\delta_{1}}{(\beta_{31}+\delta_{1})(-\beta_{21}+\beta_{32}+\delta_{2})}e^{-(\beta_{31}+\beta_{32}+\delta_{1}+\delta_{2})t}$
$\displaystyle s_{6}(t)$ $\displaystyle=0$ $\displaystyle s_{7}(t)$
$\displaystyle=\frac{\delta_{2}}{\beta_{32}+\delta_{2}}e^{-(\beta_{31}+\delta_{1})t}+\frac{\delta_{2}}{\beta_{21}-\beta_{32}-\delta_{2}}e^{-(\beta_{21}+\beta_{31}+\delta_{1})t}+\frac{\beta_{21}\delta_{2}}{(\beta_{32}+\delta_{2})(-\beta_{21}+\beta_{32}+\delta_{2})}e^{-(\beta_{31}+\beta_{32}+\delta_{1}+\delta_{2})t}$
$\displaystyle s_{8}(t)$
$\displaystyle=\frac{\beta_{21}\delta_{1}\delta_{2}}{(\beta_{31}+\delta_{1})(\beta_{21}+\beta_{31}+\delta_{1})(\beta_{32}+\delta_{2})}-\frac{\delta_{1}\delta_{2}}{(\beta_{31}+\delta_{1})(\beta_{32}+\delta_{2})}e^{-(\beta_{31}+\delta_{1})t}$
$\displaystyle\qquad-\frac{\beta_{21}\delta_{1}\delta_{2}}{(\beta_{31}+\delta_{1})(\beta_{21}+\beta_{31}-\beta_{32}+\delta_{1}-\delta_{2})(\beta_{32}+\delta_{2})}e^{-(\beta_{32}+\delta_{2})t}$
$\displaystyle\qquad-\frac{\delta_{1}\delta_{2}(2\beta_{21}+\beta_{31}-\beta_{32}+\delta_{1}-\delta_{2})}{(\beta_{21}+\beta_{31}+\delta_{1})(\beta_{21}-\beta_{32}-\delta_{2})(\beta_{21}+\beta_{31}-\beta_{32}+\delta_{1}-\delta_{2})}e^{-(\beta_{21}+\beta_{31}+\delta_{1})t}$
$\displaystyle\qquad+\frac{\beta_{21}\delta_{1}\delta_{2}}{(\beta_{31}+\delta_{1})(\beta_{21}-\beta_{32}-\delta_{2})(\beta_{32}+\delta_{2})}e^{-(\beta_{31}+\beta_{32}+\delta_{1}+\delta_{2})t}$
$\displaystyle s_{9}(t)$ $\displaystyle=0$ $\displaystyle s_{10}(t)$
$\displaystyle=\frac{\beta_{31}}{\beta_{31}-\beta_{32}-\delta_{3}}\left(-e^{-(\beta_{21}+\beta_{31}+\delta_{1})t}+e^{-(\beta_{21}+\beta_{32}+\delta_{1}+\delta_{2})t}\right)$
$\displaystyle s_{11}(t)$
$\displaystyle=\frac{\beta_{31}\delta_{1}}{(\beta_{21}+\delta_{1})(\beta_{21}+\beta_{31}-\beta_{32}+\delta_{1}-\delta_{3})}e^{-(\beta_{32}+\delta_{3})t}$
$\displaystyle\qquad+\frac{\beta_{31}\delta_{1}}{(\beta_{31}-\beta_{32}-\delta_{3})(\beta_{21}+\beta_{31}-\beta_{32}+\delta_{1}-\delta_{3})}e^{-(\beta_{21}+\beta_{31}+\delta_{1})t}$
$\displaystyle\qquad-\frac{\beta_{31}\delta_{1}}{(\beta_{21}+\delta_{1})(\beta_{31}-\beta_{32}-\delta_{3})}e^{-(\beta_{21}+\beta_{32}+\delta_{1}+\delta_{3})t}$
$\displaystyle s_{12}(t)$ $\displaystyle=0$ $\displaystyle s_{13}(t)$
$\displaystyle=\frac{-\beta_{21}^{2}\beta_{31}-\beta_{21}\beta_{31}^{2}+\beta_{21}\beta_{32}^{2}+\beta_{31}\beta_{32}^{2}+\beta_{21}\beta_{31}\delta_{2}+\beta_{31}\beta_{32}\delta_{2}+\beta_{21}\beta_{31}\delta_{3}+\beta_{21}\beta_{32}\delta_{3}}{(\beta_{21}-\beta_{32}-\delta_{2})(\beta_{21}+\beta_{31}-\delta_{2}-\delta_{3})(-\beta_{31}+\beta_{32}+\delta_{3})}e^{-(\beta_{21}+\beta_{31}+\delta_{1})t}$
$\displaystyle\qquad+\frac{\beta_{31}(-\beta_{21}-\beta_{32})}{(\beta_{21}+\beta_{32}-\delta_{2})(\beta_{31}-\beta_{32}-\delta_{3})}e^{-(\beta_{21}+\beta_{32}+\delta_{1}+\delta_{3})t}$
$\displaystyle\qquad+\frac{\beta_{21}(-\beta_{31}-\beta_{32})}{(\beta_{21}-\beta_{32}-\delta_{2})(\beta_{31}+\beta_{32}-\delta_{3})}e^{-(\beta_{31}+\beta_{32}+\delta_{1}+\delta_{2})t}$
$\displaystyle\qquad+\Bigg{(}\beta_{21}^{2}\beta_{31}+\beta_{21}\beta_{31}^{2}+\beta_{21}^{2}\beta_{32}+2\beta_{21}\beta_{31}\beta_{32}+\beta_{31}^{2}\beta_{32}+\beta_{21}\beta_{32}^{2}+\beta_{31}\beta_{32}^{2}-\beta_{21}\beta_{31}\delta_{2}-\beta_{21}\beta_{32}\delta_{2}$
$\displaystyle\qquad-\beta_{21}\beta_{31}\delta_{3}-\beta_{31}\beta_{32}\delta_{3}\Bigg{)}\Bigg{/}\Bigg{(}(\beta_{21}+\beta_{32}-\delta_{2})(\beta_{31}+\beta_{32}-\delta_{3})(\beta_{21}+\beta_{31}-\delta_{2}-\delta_{3})\Bigg{)}e^{-(\delta_{1}+\delta_{2}+\delta_{3})t}$
$\displaystyle s_{14}(t)$
$\displaystyle=\Bigg{(}\beta_{21}^{4}\beta_{31}\delta_{1}+3\beta_{21}^{3}\beta_{31}^{2}\delta_{1}+3\beta_{21}^{2}\beta_{31}^{3}\delta_{1}+\beta_{21}\beta_{31}^{4}\delta_{1}-2\beta_{21}^{3}\beta_{32}^{2}\delta_{1}-7\beta_{21}^{2}\beta_{31}\beta_{32}^{2}\delta_{1}-7\beta_{21}\beta_{31}^{2}\beta_{32}^{2}\delta_{1}$
$\displaystyle\qquad-2\beta_{31}^{3}\beta_{32}^{2}\delta_{1}+3\beta_{21}^{2}\beta_{32}^{3}\delta_{1}+6\beta_{21}\beta_{31}\beta_{32}^{3}\delta_{1}+3\beta_{31}^{2}\beta_{32}^{3}\delta_{1}-\beta_{21}\beta_{32}^{4}\delta_{1}-\beta_{31}\beta_{32}^{4}\delta_{1}+2\beta_{21}^{3}\beta_{31}\delta_{1}^{2}$
$\displaystyle\qquad+4\beta_{21}^{2}\beta_{31}^{2}\delta_{1}^{2}+2\beta_{21}\beta_{31}^{3}\delta_{1}^{2}-3\beta_{21}^{2}\beta_{32}^{2}\delta_{1}^{2}-6\beta_{21}\beta_{31}\beta_{32}^{2}\delta_{1}^{2}-3\beta_{31}^{2}\beta_{32}^{2}\delta_{1}^{2}+2\beta_{21}\beta_{32}^{3}\delta_{1}^{2}+2\beta_{31}\beta_{32}^{3}\delta_{1}^{2}$
$\displaystyle\qquad+\beta_{21}^{2}\beta_{31}\delta_{1}^{3}+\beta_{21}\beta_{31}^{2}\delta_{1}^{3}-\beta_{21}\beta_{32}^{2}\delta_{1}^{3}-\beta_{31}\beta_{32}^{2}\delta_{1}^{3}-2\beta_{21}^{3}\beta_{31}\delta_{1}\delta_{2}-4\beta_{21}^{2}\beta_{31}^{2}\delta_{1}\delta_{2}-2\beta_{21}\beta_{31}^{3}\delta_{1}\delta_{2}$
$\displaystyle\qquad-2\beta_{21}^{2}\beta_{31}\beta_{32}\delta_{1}\delta_{2}-4\beta_{21}\beta_{31}^{2}\beta_{32}\delta_{1}\delta_{2}-2\beta_{31}^{3}\beta_{32}\delta_{1}\delta_{2}+2\beta_{21}^{2}\beta_{32}^{2}\delta_{1}\delta_{2}+9\beta_{21}\beta_{31}\beta_{32}^{2}\delta_{1}\delta_{2}+6\beta_{31}^{2}\beta_{32}^{2}\delta_{1}\delta_{2}$
$\displaystyle\qquad-2\beta_{21}\beta_{32}^{3}\delta_{1}\delta_{2}-3\beta_{31}\beta_{32}^{3}\delta_{1}\delta_{2}-3\beta_{21}^{2}\beta_{31}\delta_{1}^{2}\delta_{2}-3\beta_{21}\beta_{31}^{2}\delta_{1}^{2}\delta_{2}-3\beta_{21}\beta_{31}\beta_{32}\delta_{1}^{2}\delta_{2}-3\beta_{31}^{2}\beta_{32}\delta_{1}^{2}\delta_{2}$
$\displaystyle\qquad+2\beta_{21}\beta_{32}^{2}\delta_{1}^{2}\delta_{2}+4\beta_{31}\beta_{32}^{2}\delta_{1}^{2}\delta_{2}-\beta_{21}\beta_{31}\delta_{1}^{3}\delta_{2}-\beta_{31}\beta_{32}\delta_{1}^{3}\delta_{2}+\beta_{21}^{2}\beta_{31}\delta_{1}\delta_{2}^{2}+\beta_{21}\beta_{31}^{2}\delta_{1}\delta_{2}^{2}$
$\displaystyle\qquad+3\beta_{21}\beta_{31}\beta_{32}\delta_{1}\delta_{2}^{2}+3\beta_{31}^{2}\beta_{32}\delta_{1}\delta_{2}^{2}-3\beta_{31}\beta_{32}^{2}\delta_{1}\delta_{2}^{2}+\beta_{21}\beta_{31}\delta_{1}^{2}\delta_{2}^{2}+2\beta_{31}\beta_{32}\delta_{1}^{2}\delta_{2}^{2}-\beta_{31}\beta_{32}\delta_{1}\delta_{2}^{3}$
$\displaystyle\qquad-2\beta_{21}^{3}\beta_{31}\delta_{1}\delta_{3}-4\beta_{21}^{2}\beta_{31}^{2}\delta_{1}\delta_{3}-2\beta_{21}\beta_{31}^{3}\delta_{1}\delta_{3}-2\beta_{21}^{3}\beta_{32}\delta_{1}\delta_{3}-4\beta_{21}^{2}\beta_{31}\beta_{32}\delta_{1}\delta_{3}-2\beta_{21}\beta_{31}^{2}\beta_{32}\delta_{1}\delta_{3}$
$\displaystyle\qquad+6\beta_{21}^{2}\beta_{32}^{2}\delta_{1}\delta_{3}+9\beta_{21}\beta_{31}\beta_{32}^{2}\delta_{1}\delta_{3}+2\beta_{31}^{2}\beta_{32}^{2}\delta_{1}\delta_{3}-3\beta_{21}\beta_{32}^{3}\delta_{1}\delta_{3}-2\beta_{31}\beta_{32}^{3}\delta_{1}\delta_{3}-3\beta_{21}^{2}\beta_{31}\delta_{1}^{2}\delta_{3}$
$\displaystyle\qquad-3\beta_{21}\beta_{31}^{2}\delta_{1}^{2}\delta_{3}-3\beta_{21}^{2}\beta_{32}\delta_{1}^{2}\delta_{3}-3\beta_{21}\beta_{31}\beta_{32}\delta_{1}^{2}\delta_{3}+4\beta_{21}\beta_{32}^{2}\delta_{1}^{2}\delta_{3}+2\beta_{31}\beta_{32}^{2}\delta_{1}^{2}\delta_{3}-\beta_{21}\beta_{31}\delta_{1}^{3}\delta_{3}$
$\displaystyle\qquad-\beta_{21}\beta_{32}\delta_{1}^{3}\delta_{3}+3\beta_{21}^{2}\beta_{31}\delta_{1}\delta_{2}\delta_{3}+3\beta_{21}\beta_{31}^{2}\delta_{1}\delta_{2}\delta_{3}+2\beta_{21}^{2}\beta_{32}\delta_{1}\delta_{2}\delta_{3}+4\beta_{21}\beta_{31}\beta_{32}\delta_{1}\delta_{2}\delta_{3}$
$\displaystyle\qquad+2\beta_{31}^{2}\beta_{32}\delta_{1}\delta_{2}\delta_{3}-4\beta_{21}\beta_{32}^{2}\delta_{1}\delta_{2}\delta_{3}-4\beta_{31}\beta_{32}^{2}\delta_{1}\delta_{2}\delta_{3}+2\beta_{21}\beta_{31}\delta_{1}^{2}\delta_{2}\delta_{3}+2\beta_{21}\beta_{32}\delta_{1}^{2}\delta_{2}\delta_{3}$
$\displaystyle\qquad+2\beta_{31}\beta_{32}\delta_{1}^{2}\delta_{2}\delta_{3}-\beta_{21}\beta_{31}\delta_{1}\delta_{2}^{2}\delta_{3}-2\beta_{31}\beta_{32}\delta_{1}\delta_{2}^{2}\delta_{3}+\beta_{21}^{2}\beta_{31}\delta_{1}\delta_{3}^{2}+\beta_{21}\beta_{31}^{2}\delta_{1}\delta_{3}^{2}+3\beta_{21}^{2}\beta_{32}\delta_{1}\delta_{3}^{2}$
$\displaystyle\qquad+3\beta_{21}\beta_{31}\beta_{32}\delta_{1}\delta_{3}^{2}-3\beta_{21}\beta_{32}^{2}\delta_{1}\delta_{3}^{2}+\beta_{21}\beta_{31}\delta_{1}^{2}\delta_{3}^{2}+2\beta_{21}\beta_{32}\delta_{1}^{2}\delta_{3}^{2}-\beta_{21}\beta_{31}\delta_{1}\delta_{2}\delta_{3}^{2}-2\beta_{21}\beta_{32}\delta_{1}\delta_{2}\delta_{3}^{2}$
$\displaystyle\qquad-\beta_{21}\beta_{32}\delta_{1}\delta_{3}^{3}\Bigg{)}\Bigg{/}\Bigg{(}(\beta_{21}+\beta_{31}-\beta_{32}+\delta_{1}-\delta_{2})(-\beta_{21}+\beta_{32}+\delta_{2})(\beta_{31}-\beta_{32}-\delta_{3})$
$\displaystyle\qquad\times(\beta_{21}+\beta_{31}-\beta_{32}+\delta_{1}-\delta_{3})(\beta_{21}+\beta_{31}-\delta_{2}-\delta_{3})(\beta_{21}+\beta_{31}+\delta_{1}-\delta_{2}-\delta_{3})\Bigg{)}e^{-(\beta_{21}+\beta_{31}+\delta_{1})t}$
$\displaystyle\qquad-\frac{\beta_{21}\beta_{32}\delta_{1}}{(\beta_{31}+\delta_{1})(\beta_{21}+\beta_{31}-\beta_{32}+\delta_{1}-\delta_{2})(\beta_{32}-\delta_{3})}e^{-(\beta_{32}+\delta_{2})t}$
$\displaystyle\qquad-\frac{\beta_{31}\beta_{32}\delta_{1}}{(\beta_{21}+\delta_{1})(\beta_{32}-\delta_{2})(\beta_{21}+\beta_{31}-\beta_{32}+\delta_{1}-\delta_{3})}e^{-(\beta_{32}+\delta_{3})t}$
$\displaystyle\qquad+\frac{\beta_{31}\left(\beta_{21}^{2}\delta_{1}+2\beta_{21}\beta_{32}\delta_{1}+\beta_{32}^{2}\delta_{1}+\beta_{21}\delta_{1}^{2}+\beta_{32}\delta_{1}^{2}-\beta_{32}\delta_{1}\delta_{2}\right)}{(\beta_{21}+\delta_{1})(\beta_{21}+\beta_{32}-\delta_{2})(\beta_{21}+\beta_{32}+\delta_{1}-\delta_{2})(\beta_{31}-\beta_{32}-\delta_{3})}e^{-(\beta_{21}+\beta_{32}+\delta_{1}+\delta_{3})t}$
$\displaystyle\qquad+\frac{\beta_{21}\left(\beta_{31}^{2}\delta_{1}+2\beta_{31}\beta_{32}\delta_{1}+\beta_{32}^{2}\delta_{1}+\beta_{31}\delta_{1}^{2}+\beta_{32}\delta_{1}^{2}-\beta_{32}\delta_{1}\delta_{3}\right)}{(\beta_{31}+\delta_{1})(\beta_{21}-\beta_{32}-\delta_{2})(\beta_{31}+\beta_{32}-\delta_{3})(\beta_{31}+\beta_{32}+\delta_{1}-\delta_{3})}e^{-(\beta_{31}+\beta_{32}+\delta_{1}+\delta_{2})t}$
$\displaystyle\qquad+\Bigg{(}\frac{\beta_{31}(\beta_{21}+\beta_{32})}{(\beta_{21}+\beta_{32}+\delta_{1}-\delta_{2})(\beta_{21}+\beta_{31}+\delta_{1}-\delta_{2}-\delta_{3})}$
$\displaystyle\qquad+\frac{\beta_{31}\beta_{32}\delta_{1}}{(\beta_{32}-\delta_{2})(\beta_{21}+\beta_{32}+\delta_{1}-\delta_{2})(\beta_{21}+\beta_{31}+\delta_{1}-\delta_{2}-\delta_{3})}$
$\displaystyle\qquad+\frac{\beta_{21}(\beta_{31}+\beta_{32})}{(\beta_{31}+\beta_{32}+\delta_{1}-\delta_{3})(\beta_{21}+\beta_{31}+\delta_{1}-\delta_{2}-\delta_{3})}$
$\displaystyle\qquad-\frac{\beta_{21}\beta_{32}\delta_{1}}{(\beta_{32}-\delta_{3})(\beta_{31}+\beta_{32}+\delta_{1}-\delta_{3})(-\beta_{21}-\beta_{31}-\delta_{1}+\delta_{2}+\delta_{3})}\Bigg{)}e^{-(\delta_{2}+\delta_{3})t}$
$\displaystyle\qquad-\Bigg{(}\beta_{21}^{2}\beta_{31}+\beta_{21}\beta_{31}^{2}+\beta_{21}^{2}\beta_{32}+2\beta_{21}\beta_{31}\beta_{32}+\beta_{31}^{2}\beta_{32}+\beta_{21}\beta_{32}^{2}+\beta_{31}\beta_{32}^{2}-\beta_{21}\beta_{31}\delta_{2}-\beta_{21}\beta_{32}\delta_{2}$
$\displaystyle\qquad-\beta_{21}\beta_{31}\delta_{3}-\beta_{31}\beta_{32}\delta_{3}\Bigg{)}\Bigg{/}\Bigg{(}(\beta_{21}+\beta_{32}-\delta_{2})(\beta_{31}+\beta_{32}-\delta_{3})(\beta_{21}+\beta_{31}-\delta_{2}-\delta_{3})\Bigg{)}e^{-(\delta_{1}+\delta_{2}+\delta_{3})t}$

We have omitted the remaining states $s_{15},\ldots,s_{26}$, since they are
rather involved. Instead, the formula for the prevalence $y(t)$ is shown
below, which is, unfortunately, also rather involved.

$\displaystyle 3y(t)$
$\displaystyle=e^{-\delta_{1}t}-\frac{\beta_{21}\delta_{3}}{(\beta_{21}+\delta_{1}-\delta_{2})(\beta_{32}+\delta_{3})}e^{-(\beta_{21}+\delta_{1})t}-\frac{\beta_{31}\delta_{2}}{(\beta_{32}+\delta_{2})(\beta_{31}+\delta_{1}-\delta_{3})}e^{-(\beta_{31}+\delta_{1})t}$
$\displaystyle\qquad+\Bigg{(}\beta_{32}(-\beta_{21}^{5}+\beta_{21}^{4}(-5\beta_{31}+4\beta_{32}-3\delta_{1}+2\delta_{2}+3\delta_{3})-\beta_{31}(\beta_{31}-\beta_{32}+\delta_{1}-\delta_{2})$
$\displaystyle\qquad\times(\beta_{31}-\beta_{32}-\delta_{3})\left(\beta_{31}^{2}+(\delta_{1}-2\delta_{2})(\delta_{1}-\delta_{3})+\beta_{32}(-2\delta_{1}+\delta_{2}+\delta_{3})-\beta_{31}(2\beta_{32}-2\delta_{1}+2\delta_{2}+\delta_{3})\right)$
$\displaystyle\qquad-\beta_{21}^{3}(10\beta_{31}^{2}+5\beta_{32}^{2}+3\delta_{1}^{2}-5\delta_{1}\delta_{2}+\delta_{2}^{2}-6\delta_{1}\delta_{3}+6\delta_{2}\delta_{3}+2\delta_{3}^{2}+\beta_{32}(-9\delta_{1}+6\delta_{2}+8\delta_{3})$
$\displaystyle\qquad-\beta_{31}(16\beta_{32}-10\delta_{1}+9\delta_{2}+11\delta_{3}))+\beta_{21}^{2}(-10\beta_{31}^{3}+2\beta_{32}^{3}-\delta_{1}^{3}+4\delta_{1}^{2}\delta_{2}-2\delta_{1}\delta_{2}^{2}+3\delta_{1}^{2}\delta_{3}$
$\displaystyle\qquad-9\delta_{1}\delta_{2}\delta_{3}+3\delta_{2}^{2}\delta_{3}-2\delta_{1}\delta_{3}^{2}+4\delta_{2}\delta_{3}^{2}+\beta_{32}^{2}(-8\delta_{1}+5\delta_{2}+6\delta_{3})+\beta_{32}(6\delta_{1}^{2}-10\delta_{1}\delta_{2}+2\delta_{2}^{2}$
$\displaystyle\qquad-11\delta_{1}\delta_{3}+11\delta_{2}\delta_{3}+3\delta_{3}^{2})+\beta_{31}^{2}(24\beta_{32}-14\delta_{1}+15(\delta_{2}+\delta_{3}))-\beta_{31}(15\beta_{32}^{2}+5\delta_{1}^{2}-14\delta_{1}\delta_{2}$
$\displaystyle\qquad+4\delta_{2}^{2}-15\delta_{1}\delta_{3}+18\delta_{2}\delta_{3}+5\delta_{3}^{2}+\beta_{32}(-25\delta_{1}+20\delta_{2}+22\delta_{3})))-\beta_{21}(5\beta_{31}^{4}-\beta_{31}^{3}(16\beta_{32}-10\delta_{1}$