Modeling the Spread of Multiple Contagions on Multilayer Networks

TL;DR

This paper introduces a novel SIS model for multiple contagions on multilayer networks, incorporating various spreading channels and mutations, with analytical estimates of the basic reproduction number and insights into strain persistence without epidemic thresholds.

Contribution

It develops a new multilayer network SIS model for multiple contagions, analyzing mutation-driven strain persistence and epidemic thresholds, applicable to diverse spreading phenomena.

Findings

Analytical estimation of the basic reproduction number.

Identification of strain persistence without epidemic thresholds.

Model applicability to various spreading phenomena.

Abstract

A susceptible-infected-susceptible (SIS) model of multiple contagions on multilayer networks is developed to incorporate different spreading channels and disease mutations. The basic reproduction number for this model is estimated analytically. In a special case when considering only compartmental models, we analytically analyze an example of a model with a mutation driven strain persistence characterized by the absence of an epidemic threshold. This model is not related to the network topology and can be observed in both compartmental models and models on networks. The novel multiple-contagion SIS model on a multilayer network could help in the understanding of other spreading phenomena including communicable diseases, cultural characteristics, addictions, or information spread through e-mail messages, web blogs, and computer networks.

Figures

Tables

Table 1. Table 1: Summary of notation that is used throughout the paper.

$R$	Set of relation types.
$E_l$	Edge set of layer $l$.
$A_l$	Adjacency matrix of layer $l$.
$H$	Set of contagions.
${\beta }_l^c$	Infection rate of contagion $c$ on layer $l$.
${\mu }_{cb}$	Mutation rate of contagion $c$ into contagion $b$.
${\gamma }_c$	Healing rate of contagion $c$.
$p_i^L\left(t\right)$	Probability node $i$ is infected by all contagions
	in $L$ at time $t$.
$p^L={\sum }_{i=1}^N{p}_i^L$	Population infected by $L$.
$L_b^c$	Set $L$ including $c$, excluding $b$: $L\setminus b\cup c$.
$R_0^c$	Basic reproduction number for contagion $c$
${\beta }_{\text{c}r}^c$	Critical value for $\beta$ of contagion $c$.

Equations

$$R_0=\frac{{\lambda }_{\text{max}}(A)\beta }{\gamma }\equiv \frac{{\beta }_{\mathrm{eff}}}{\gamma },$$

$${\beta }_{cr}=\frac{\gamma }{{\lambda }_{\text{max}}(A)}.$$

$${\dot{p}}_i^c=\sum _{l\in H_c}{p}_i^{cl}{\gamma }_l+p_i^l{\mu }_{lc}+\left(1-\sum _{l\in \mathcal{P}(H)\setminus \varnothing }{p}_i^l\right)f_i^c-p_i^c\left({\gamma }_c+\sum _{l\in H}{\mu }_{cl}+f_i^l\right),$$

$${\dot{p}}_i^L=\sum _{c\in H_L}\left(p_i^{L^c}{\gamma }_c+\sum _{b\in L}{p}_i^{L_b^c}{\mu }_{cb}\right)+\sum _{c\in L}{p}_i^{L_c}{f}_i^c-p_i^L{\gamma }_c-p_i^L\sum _{c\in H_L}\left(f_i^c+\sum _{b\in L}{\mu }_{bc}\right).$$

$${\dot{p}}_i^H=\sum _{c\in H}{p}_i^{H_c}{f}_i^c-p_i^H\sum _{c\in H}{\gamma }_c.$$

$$f_i^c=\sum _{\begin{array}{c}j=1\end{array}}^n\sum _{l\in R}\sum _{\begin{array}{c}A\in P(H)\\ [c\in A]\end{array}}{A}_{ij}^l{\beta }_l^c{p}_j^A.$$

$$\frac{d{ϵ}_i^c}{dt}$$

$$=\sum _{l\in H}{ϵ}_i^l{\mu }_{lc}+\sum _{j=1}^n{ϵ}_j^c{t}_{i,j}^c-{ϵ}_i^c\left({\gamma }_c+\sum _{l\in H}{\mu }_{cl}\right)$$

$$t_{i,j}^c=\sum _{l\in R}{A}_{ij}^l{\beta }_l^c$$

$$\frac{dϵ}{dt}=\mathbf{D}ϵ$$

$$Re\{{\lambda }_1(\mathbf{D})\}<0,$$

$${\lambda }_1(\mathbf{D})<0,$$

$$\frac{d{p}_i^c}{dt}$$

$$\frac{d\mathbf{p}}{dt}\le \mathbf{Dp},$$

$$\frac{d{\mathbf{p}}^{\prime }}{dt}=\mathbf{D}{\mathbf{p}}^{\prime },$$

$$0\le \underset{t\to \infty }{\lim }{p}_i^L\le \underset{t\to \infty }{\lim }\prod _{c_r\in L}{p}_{,i}^{c_r}\prod _{c_r\in H\backslash L}(1-p_i^{c_r})=0\times 1=0$$

$$\sum _c\frac{d{ϵ}_i^c}{dt}=\sum _{j=1}^n\sum _{l\in R}{A}_{ij}^l\sum _c{ϵ}_j^c{\beta }_l^c-\sum _c{ϵ}_i^c{\gamma }_c+\sum _c{ϵ}_i^c\sum _{l\in H^c}\left({\mu }_{lc}-{\mu }_{cl}\right)$$

$$\frac{d{ϵ}_i}{dt}=\sum _{j=1}^n\sum _{l\in R}{A}_{ij}^l\sum _c{ϵ}_j^c{\beta }_l^c-\sum _c{\gamma }_c{ϵ}_i^c$$

$$\frac{d{ϵ}_i}{dt}\le \sum _{j=1}^n{ϵ}_j\sum _{l\in R}{A}_{ij}^l\underset{c}{\max }{\beta }_l^c-\underset{c}{\min }{\gamma }_c{ϵ}_i,$$

$$\frac{dϵ}{dt}$$

$$\frac{d{ϵ}^{\sup }}{dt}$$

$${\lambda }_1({\mathbf{Q}}^{\sup })={\lambda }_1({\mathbf{T}}^{\sup })-{\gamma }^{\sup }$$

$${\lambda }_1({\mathbf{T}}^{\sup })<{\gamma }^{\sup }$$

$$R^{+}=\frac{{\lambda }_1({\mathbf{T}}^{\sup })}{{\gamma }^{\sup }},$$

$$\frac{dϵ}{dt}$$

$$\frac{d{ϵ}^{\inf }}{dt}$$

$${\lambda }_1({\mathbf{T}}^{\inf })>{\gamma }^{\inf }$$

$$R^{-}=\frac{{\lambda }_1({\mathbf{T}}^{\inf })}{{\gamma }^{\inf }},$$

$$\frac{d{ϵ}_{\mathbf{i}}}{dt}$$

$${ϵ}_i^c\approx u_c{ϵ}_i,$$

$$\frac{d{ϵ}_i}{dt}$$

Full text

††journal: Physica A: Statistical Mechanics and its Applications

1 Introduction

Epidemiological models, developed as tools for analyzing the spread and control of infectious diseases, have also been adapted to model the dynamics of contagious entities as diverse as communicable diseases, cultural characteristics (such as religious beliefs, fads or innovations), addictions, or information spread (through rumors, e-mail messages, web blogs, peer-to-peer computer networks, etc). The susceptible-infected-susceptible (SIS) model is one of the simplest and well-studied model for emerging disease outbreaks (like influenza, chlamydia, gonhorrea, etc.) that does not give immunity upon recovery [1, 2, 3, 4, 5]. The SIS compartmental model divides the population into two compartments (classes) – susceptible to the infection of the pathogen (often denoted by S) and infected by the pathogen (given the symbol I). The disease is transmitted with a rate $\beta$ while infected individuals become susceptible again with a rate $\gamma$. The dynamics of the infectious class depends on a quantity called basic reproduction number, $R_{0}$, defined as the number of secondary infections caused by a single infective introduced into a population made up entirely of susceptible individuals over the course of the infection of this single infective. For the SIS model in a well-mixed population, the basic reproduction number is equal to $R_{0}={\beta}/{\gamma}$ [6]. By adopting so called quenched mean-field theory [7, 8], the basic reproduction number for the SIS model on networks has been estimated to be equal to

[TABLE]

where $A$ is the adjacency matrix of the network and $\lambda_{\rm max}(A)$ is the largest eigenvalue of the matrix $A$. The network structure, described with (expressed through) the largest eigenvalue of the network adjacency matrix, is encoded in the effective transmission rate $\beta_{\rm eff}=\lambda_{\text{max}}(A)\beta$. For the SIS model in homogeneous networks, such as the Erdos-Renyi random networks and random regular graphs, $\lambda_{\rm max}(A)=\langle k\rangle$. According to Ref. [9], for scale-free networks, $1/{\lambda_{\rm max}}=1/\sqrt{k_{\rm max}}$ if $\hat{\gamma}>5/2$, and $1/{\lambda_{\rm max}}=\langle k\rangle/\langle k^{2}\rangle$ if $2<\hat{\gamma}<5/2$, where the degree distribution follows $P(k)\sim k^{-\hat{\gamma}}$

In the thermodynamic limit when the number of nodes approaches infinity, a phase-transition occurs at $R_{0}=1$. If the basic reproduction number falls below the critical value ($R_{0}<1$), the infection dies out. For $R_{0}>1$ there is an epidemic in the population. At the phase transition, threshold for the transmission rate $\beta$ equals

[TABLE]

Note that in Eq. (2) for a fixed value of $\gamma$, large value of $\lambda_{max}(A)$ implies vanishing threshold ($\beta_{cr}\to 0$). In other words, the network structure causes an absence of an epidemic threshold and its associated critical behavior.

Models for multiple diseases that co-evolve in a network has recently been gaining attention [10, 11, 12, 13]. However, these diseases are often assumed to be mutually exclusive. While such models are usually discussed in the context of epidemics, they are more aptly used in studying belief propagation or product adoption, for example, in modeling competition in politics or competition in a marketplace. The generalization of the SIS model to arbitrary number of multiple contagions, however, has not yet been developed. In this paper we model the spread of multiple contagions on networks. We further assume that each contagion spreads over different spreading channels resulting in a multiple-contagion SIS model on a multilayer network.

1.1 Literature overview

1.1.1 Multiple-contagion compartmental models

Developing models for interacting strains of the same pathogen, such as influenza [14] or dengue [15], or interacting diseases such as HIV/AIDS and malaria [16], is one of the most theoretically challenging problems in infectious disease epidemiology. The central problem comes from the explosive growth in the number of state variables of the system with the linear increase in the number of strains or pathogens [17].

The spread of multiple diseases can be characterized by various factors such as contact with infected individuals, the interplay between diseases coupled by mutation [18] or cross-immunity [19, 20, 21, 18], and coexistence [20] or by the principle of competitive exclusion [22, 23] which states that the strain with the largest reproduction number drives other strains into extinction. Epidemiological models that study multiple strains fall into super-infection models [24, 25, 26] where strains cannot coexist in a host because the most virulent strain takes over, and co-infection models [27, 25] where coexistence is possible. For a discussion on how multiple infections have been modeled in evolutionary epidemiology, as well as co-and super-infection models in the standard setting, the reader is referred to [24]. An epi-evolutionary model on a $n$-regular network with $N$ pathogens which are allowed to mutate into each other is proposed in [28]. Besides spreading dynamics, the model takes into account the evolutionary properties of a parasite by keeping track of the change of the mean value. Authors demonstrate that the coupled dynamics lead to transient phenomena that cannot be described by standard invasion analyses.

Cross-immunity can be defined in at least two ways: one is reduced susceptibility and the other is reduced transmission (less likely to transmit) [29]. The biological basis of this is that multiple influenza infections occur during the life of individuals and one’s increased/decreased susceptibility depends on this history. The dimensionality of a history-based system is $O(2^{N})$ (the power set of the pathogens) and multiple methods have been proposed for its reduction, focusing on a symmetrized system [30, 21] or via age-structure by grouping the population into compartments which have seen a strain $i$ [31, 32, 33]. The work of [34] focuses on the stability analysis of a model of $n-$compartments, each having its own reproductive ratio and endemic threshold. Coexistence has been inspected in [35] when the strains are coupled by mutation.

1.1.2 Spreading processes on multilayer networks

Multilayer networks are fundamental for the understanding of dynamical processes on networked systems, including, for example, spreading processes, such as flows (and congestion) in transportation networks [36, 37], and information and disease spreading in social networks [38, 39, 40, 41, 42].

As reviewed in [43], there are two different categories of dynamical processes on multilayer networks: a single dynamical process and mixed (or coupled) dynamics, in which two or more dynamical processes, defined on each layer separately, are coupled together with inter-layer connections between nodes. One of the simplest types of dynamics is a diffusion process. In a random walk, a discrete diffusion process, a walker jumps between nodes. In a multilayer network, the walker switches between layers via an inter-layer edge, resulting in enriched random-walk dynamics [44, 45, 46]. The continuous diffusion process has also been analyzed in multiplex networks [47, 48] and novel phenomenon has been observed: diffusion can be faster in a multiplex network than in any of the layers considered independently. Congestion in multilayer networks has been studied recently [49, 37] for modeling multimodal transportation systems.

Coupled spreading processes on multilayer networks have recently been analyzed, including spreading dynamics of two concurrent diseases in two-layer multiplex networks [50, 51, 41, 52, 53, 54] and spread of disease coupled with the spread of information or behavior [55, 40, 56, 38, 39]. It was observed that two spreading processes can enhance each other (for example, one disease facilitates infection by the other), or one process can inhibit the spread of the other (for example, a disease can inhibit infection by another disease or the spreading of awareness about a disease can inhibit the spread of the disease).

1.2 Motivation and our contribution

This paper addresses multiple-contagion SIS model on multilayer networks. So far only multiple-contagion compartmental models have been studied as reviewed in the previous subsection. However, a simple example illustrates that the SIS model on a complex network produces different dynamics from the compartmental SIS model with the same parameter values. Neglecting the network structure of the model and working only with the compartment model can lead to completely different results concerning how the infected populations evolve in time. Figure 1 depicts the infected sub-populations over time for two models: two-contagion compartmental SIS model and a two-contagion SIS model on a two-layer Erdos-Renyi random network. The number of nodes in the network model is equal to the total population in the compartmental model and both models have same parameter values. The equations governing the evolution of both the compartmental and network models are given in the next sections – this is just an illustration that both models can exhibit different behavior and understanding multiple-contagion models on multilayer networks is the next step in modeling various spreading phenomena. The main contribution of this work is twofold: (1) a multiple-contagion SIS model on multilayer networks is developed; and (2) several properties of the model are analytically discussed including upper/lower bounds of basic reproduction number, rapid mutation processes and competition. Furthermore, we discuss how the epidemic threshold depends on the network structure. In particular, two classes of the model are studied in detail. In the first class, contagions do not compete and a node can be infected by an arbitrary number of strains, while the second class describes competing epidemics on networks implying that each node (agent) can only be infected by a single contagion. Coinfection with multiple pathogens is a common occurrence. Some known combinations are HIV and tuberculosis, HIV and hepatitis, HIV and malaria. To the best of our knowledge, an epidemic threshold for two-contagions models and an arbitrary number of network layers has not yet been derived. In a special case, we analytically analyze an example of a model with a mutation driven strain persistence characterized by the absence of an epidemic threshold.

2 Model description

We model multiple contagions on multilayer networks with three graphs: a multiplex graph, a bipartite graph, and a directed graph. The multiplex graph represents a social network consisting of a set of social actors and sets of dyadic ties. A multiplex graph/network is defined as a collection of a set of nodes $V$, a set $R=\{1,\ldots,r\}$ of relation types, and for each $l\in{R}$, a set $E_{l}$ of edges describing the presence or absence of edges of type $l$ between pairs of nodes. The graph $(V,E_{l})$ is also called a layer; we write $A_{l}=(A_{ij}^{l})$ for the adjacency matrix of this graph. We assume that each layer represents a different spreading channel through which contagions are spread. For example, in etiology, scientists have recognized five major modes of disease transmission: airborne, waterborne, bloodborne, by direct contact, and through vector (insects or other creatures that carry germs from one species to another). These modes of disease transmission not necessarily spread over the same ties between individuals. Let the set $H=\{1,2\ldots,m\}$ label the $m$ strains present in the system. For the bipartite graph $G=(R\cup H,E)$, the vertex set $R$ is the set of layers/channels, the vertex set $H$ is the set of contagions, and the set $E$ is the set of edges such that every edge connects a vertex in $R$ to one in $H$. In a classical SIS model, the parameter $\beta$ represents contact or infection rate of a disease. Here $\beta^{c}_{l}$ denotes the contact or infection rate of the contagion/disease $c$ through the layer/channel $l$. The bipartite graph can be represented with a weighted adjacency matrix defined by $[\beta^{c}_{l}]$. We assume that contagions can be changed/mutated and be transformed between each other. For contagious diseases (also called communicable diseases), scientists have documented that diseases could mutate to become more contagious. This is modeled as a directed graph with an $m\times m$ adjacency matrix $[\mu_{cb}]$; $\mu_{cb}$ is a rate at which $c$ mutates into $b$. In general, we assume that $\mu_{cb}\neq\mu_{bc}$.

In the classical SIS model, when a single contagion is presented in the system, the node’s state is defined with two variables: the probability that the node $i$ is susceptible $p^{\emptyset}_{i}$ and the probability that the node $i$ is infected $p^{I}_{i}$. Clearly, these two variables are not independent since $p^{\emptyset}_{i}+p^{I}_{i}=1$. Here, for $m$ co-circulating contagions, the number of variables each node has is $2^{m}$, out of which $2^{m}-1$ are independent. The power set of ${H}$, ${\cal P}({H})$, is the set of all subsets of ${H}$, including the empty set and ${H}$ itself. Let $p^{L}_{i}$ be the probability that the node $i$ is currently infected (contaminated) by all contagious entities in the set ${L}$, ${L}\in{\cal P}({H})$. Clearly, $\sum_{{L}\in{\cal P}({H})}p^{L}_{i}=1$. The current contagious set of a node can be changed by one of three mechanisms. The first mechanism is mutation: a contagion mutates into another contagion. We assume that this mechanism is described with a rate: $\mu_{cb}$ denotes the rate at which $c$ mutates into $b$. The second mechanism is analogous to the node healing from contagion $c$ and $\gamma^{c}$ represents the recovery rate. The third mechanism is transmission and is analogous to a node being infected with a disease from its neighbors. This mechanism is the only network-induced mechanism resulting in the change of the set ${L}$ due to contact with the neighbors and we also refer to it as the contact-induced change mechanism of the set ${L}$. This mechanism is described with a quantity $\beta^{c}_{l}$. Table 1 summarizes all quantities used in the paper.

For $i=1,\ldots,n$, $L\in{\cal P}({H})$, and $|H|>|L|>1$) the model consists of three parts. Equation (3) captures the dynamics of node $i$ when it is infected by only one contagion (denoted by c). The first term describes the rate at which node $i$ recovers from contagion $l\neq c$ if it is infected with both $c$ and $l$. The second term describes the rate at which contagion $l$ mutates into contagion $c$ if node $i$ is infected by contagion $l$. The third term describes the rate at which node $i$ will be infected by contagion $c$ if it is susceptible. The fourth and fifth term describe the rate at which node $i$ will either recover from $c$ or $c$ will mutate into $l\neq c$, if it is infected with $c$. The fifth term describes the rate at which node $i$ will be infected by another contagion $l$, if it is already infected with $c$.

[TABLE]

where $f^{c}_{i}$ and $f^{l}_{i}$ are quantities defined with (6). Equation (4) has a more complex mechanism. The first term describes the rate at which node $i$ recovers from $c$ if it is infected by all contagions in $L$ and also $c$. The second term describes the rate at which contagion $c\in H_{L}$ will mutate into contagion $b\in L$, if node $i$ is infected with contagions $L\cup c\setminus b$. The third term describes the rate at which node $i$ will be infected by contagion $c$ if it is already infected with contagions $L\setminus c$. The fourth term is the rate at which node $i$ will recover from $c$ if it is infected by all contagions in $L$. The fifth term is the rate at which node $i$ will be infected by contagion $c\in H_{L}$ if it is infected will all contagions in $L$. The sixth term describes the rate at which contagion $c\in L$ will mutate into contagion $b\in L$.

[TABLE]

Equation (5) describes the case when the node is infected by all contagions. The first term is the rate at which node $i$ will become infected by $c$ if it is infected by all contagions in $H_{c}$. The second term describes the rate at which node $i$ will recover from one contagions, if it is infected by all contagions.

[TABLE]

The contact mechanism defined in equation (6) goes over all neighbors of node $i$ in all layers. Contagion $c$ will be transmitted with corresponding rate $\beta^{c}_{l}$ for a layer/channel $l$ for those neighbors of $i$ which are infected by $c$, that is $c\in A$ where $A\in P(H)$ is the contagion set.

[TABLE]

2.1 Simulation on Networks

In this subsection we describe the Monte Carlo simulation procedure of the multiple contagion SIS model on multiple networks. The simulation is initialized by setting a certain percentage of nodes as infected. Then, for $t=1,\ldots,T$, every infected (assume infected with all contagions in the set $L\subseteq H$) node tries to infect their neighbors with every contagion $c\in L$, via every channel. This is done by checking whether the neighbor is not infected with contagion $c$. If they aren’t, the contagion is passed to the neighbor with probability $\Delta t\beta^{c}_{l}$. After this contact mechanism is finished, every infected node can either recover from contagion $c\in L$ with probability $\Delta t\gamma_{c}$, contagion $c\in L$ can mutate into contagion $d\in H\setminus L$ with probability $\Delta t\mu_{cd}$, or nothing occurs, and the node stays infected with the contagions in $L$. The corresponding continuous-time equations can be derived from this model leading to equations (3) – (5). This derivation is, however, out of scope of the paper; a more detailed comparison of deterministic and stochastic models will be provided in a forthcoming paper. Figure 2 shows a comparison between the aforementioned Monte Carlo simulation and the fourth order Runge-Kutta method applied to equations (3) – (5) with $\Delta t=0.01$. The comparison is carried out on a two-contagion SIS model on two Erdos-Renyi random networks. A node in the network can be in either a susceptible state (S) or in one of the three infected states denoted as 1, 2, or 12 (12 means the node is infected both with contagions 1 and 2).

3 Properties of the model

The model is a generalization of the classical SIS model on a single layer network. The equations (3)-(5) assume that the population remains constant, i.e., the concentrations satisfy the closure condition $\sum_{L}p^{L}_{i}(t)=1$ for every $t$ and every $i$. The equilibrium is obtained by setting $\dot{p}^{L}_{i}=0$. By acting so, the straightforward conclusion is that the point $p^{L}_{i}=0$, $L\in P(H)\backslash{\emptyset}$, $p^{L}_{i}=1$ $i=1,..,n$, i.e the point of epidemic origin, is an equilibrium point for the dynamical system (3), (4), (5), and (6).

It is a standard approach in the analysis of spreading processes occurring on complex networks, to derive the epidemic threshold relation, from the stability analysis of the point of epidemic origin. Locally stable epidemic origin implies that all initial perturbations of the system around the equilibrium point (epidemic outbreak) tend towards the origin (epidemic extinction) as time elapses. In the contrary, an unstable equilibrium point suggests that the outbreak will fully develop into persistent epidemic.

Due to the closure condition in the vicinity of the equilibrium, we restrict the stability analysis on the system of equations (4). By performing a classical perturbation analysis around the point of epidemic origin, one obtains:

[TABLE]

where

[TABLE]

The variables $\epsilon^{c}_{i}$ represent the probabilities that the node $i$ is infected by the strain $c$ at given instance in time, near the epidemic origin. One should note that the behavior of the system (3), (4), (5), and (6) near the epidemic origin restricts to the behaviour of a variant system in which only a single strain infection at each node $i$ may be present at a given time $t$.

This simplifies the stability analysis of the system of equations (4) from a system of dimension $n(2^{m}-1)$ to dimension $nm$.

Consider the vector $\mathbf{\epsilon}=[\epsilon^{i}_{j}]^{T}$. Let $\mathbf{D}=[d_{q,r}]$ be a $nm\times nm$ square matrix, which elements $d_{q,r}$ satisfy the following conditions:

$d_{q,r}=t_{i,j}^{c}$, if $q=cn+i$ and $r=cn+j$, $i\neq j$

• 2.

$d_{q,r}=\mu_{lc}$, if $q=cn+i$ and $r=ln+i$

• 3.

$d_{q,r}=-\gamma^{c}-\left(\sum_{l\in H}\mu_{cl}\right)$, if $q=r=cn+i$

• 4.

0, otherwise.

Under these conditions, the system of equations (7), may be written in the following vector form:

[TABLE]

Before we proceed, notice that the matrix $\mathbf{D}$ is characterized by negative diagonal elements, while all off-diagonal entries are non-negative. For completeness, we will state and proof the following statement: Matrix $\mathbf{D}$ has a distinct eigenvalue $\lambda_{1}(\mathbf{D})$, such that:

$\lambda_{1}(\mathbf{D})$ is real;

• 2.

$\Re\{\lambda_{j}(\mathbf{D})\}<\lambda_{1}(\mathbf{D})$, $j=\overline{2,nm}$

We will refer to $\lambda_{1}(\mathbf{D})$ as the largest eigenvalue of $\mathbf{D}$. Since all diagonal entries are finite, one may find $w>0,w\in\mathbb{R}$, such that the matrix $\mathbf{A}=w\mathbf{I}+\mathbf{D}\geq 0$. Since the analyzed network is well connected, matrix $\mathbf{A}$ is irreducible, and accordingly, by the Perron-Frobenius theorem for non-negative irreducible matrices, has one distinctive eigenvalue $\lambda_{1}(\mathbf{A})$ such that:

$\lambda_{1}(\mathbf{A})$ is real and positive;

• 2.

$||\lambda_{j}(\mathbf{A})||\leq||\lambda_{1}(\mathbf{A})||$, $j=\overline{2,nm}$

• 3.

$\Re\{\lambda_{j}(\mathbf{A})\}<\lambda_{1}(\mathbf{A})$, $j=\overline{2,nm}$

Since all eigenvectors of $\mathbf{D}$ are simultaneously eigenvectors of $\mathbf{A}$ as well, corresponding eigenvalues are related trough the equation $\lambda_{i}(\mathbf{A})=w+\lambda_{i}(\mathbf{D})$. From this equation it is straightforward clear that $\lambda_{1}(\mathbf{D})$ is real and $\Re\{\lambda_{j}(\mathbf{D})\}<\lambda_{1}(\mathbf{D})$, $j=\overline{2,nm}$. The proof is completed.

From linear systems theory it is a well known that the dynamical system described with the equation (9) converges towards the origin, $\mathbf{\epsilon}\rightarrow\mathbf{0}$, as time elapses, provided:

[TABLE]

or, since $\lambda_{1}(\mathbf{D})$ is real

[TABLE]

This implies that all initial perturbations of the dynamical system (3), (4), (5), and (6) around the point of epidemic origin (point $\mathbf{0}$) will tend towards the epidemic origin (subsequently the epidemic will vanish from the network), if relation (10) holds. Contrary, for $\lambda_{1}(\mathbf{D})>0$, all initial perturbations will move away from the origin and a state of persistent epidemic presence in the network exists.

We now show that the following statement holds: Providing relation (10) holds, epidemic origin is globally (on $[0,1]^{n(2^{m}-1)}$) stable point of equilibrium of the dynamical system (3), (4), (5), and (6).

The proof will be conducted on the system of relations (3), a variant that upper-bounds the system (3), (4), (5), and (6). In (3) stochastic independence of joint events is introduced in the term that represent mutation from arbitrary strain $l$ to strain $c$, provided strain $l$ is present and strain $c$ is absent at node $i$ at time $t$. Since, for stochastic events $A$ and $B$, $P(A,B)=P(B\backslash A)P(A)\leq P(A)P(B)$, implementation of stochastic independence of joint events in this case, results in an upper-bound system variant.

A similar approach is used through-out the literature that studies epidemic spreading on complex networks, primarily in relation to the term defining infection of susceptible node by infected neighbours. The latest, has already been taken into account in the definition of the system (3), (4), (5), and (6) as well, for example in the third and fifth term of equation (4). In what follows, by introducing this notion to one more term, we upper-bound relation (4) and even further upper-bound the true behaviour of the system.

Let $p^{c}_{i}\in[0,1]$ represent the probability that node i is infected with an agent from strain $c$ at an instance of time. Assuming statistical independence of joint events holds, the following system of relations may be considered:

[TABLE]

The system of relations (3) may be re-written in the following form:

[TABLE]

where $\mathbf{p}=[p_{k}]$, $k=cn+i$ is a column vector, and $\mathbf{D}$ is an $nm\times nm$ square matrix, as defined above. Consider the system:

[TABLE]

where $\mathbf{p^{\prime}}$ is an $mn\times 1$ vector. If (10) holds, then the point $\mathbf{p^{\prime}}=\mathbf{0}$ is globally stable point of equilibrium of the system (13). From the Gronwall-Bellman inequality, there exists a time instance $t_{0}$, such that for $t>t_{0}$, $\mathbf{p}\leq\mathbf{p^{\prime}}$ and, since $\mathbf{p}\geq 0$, $\mathbf{p}\rightarrow 0$, when $t\rightarrow\infty$.

Let $L=\{c_{a1},..c_{as}\}$ represent a particular infectious event comprising of $s$ different strain infections. Then:

[TABLE]

The proof is completed.

Equation (10) defines the epidemic threshold for the dynamical system (3), (4), (5), and (6). Though the obtained result is accurate, it gives no straightforward relation between infection parameters of the model and the structure of the network itself. Assessment of the epidemic state of the network requires numerical computation of $\lambda_{1}(\mathbf{D})$, that for spreading processes involving a large number of strains, might be computationally intensive.

In what follows, we present an analysis, that under some approximations, leads to expressions that roughly estimate the epidemic state of the system as an explicit function of epidemic parameters and the structure of the network. We start by grouping the system of equations (7) into $n$ groups of $m$ equations and then by summing the equations in each group:

[TABLE]

Notice that the third term in (15), $\sum_{c}\epsilon^{c}_{i}\sum_{l\in H^{c}}\left(\mu_{lc}-\mu_{cl}\right)=0$. Let $\epsilon_{i}=\sum_{c}\epsilon^{c}_{i}$ be the perturbation of the cumulative (all strains considered) infection probability of node $i$ in the vicinity of the epidemic origin. Subsequently, the relation (15) may be rewritten in the following form:

[TABLE]

One should note that the interpretation of $\epsilon_{i}$ as a cumulative infection probability is valid only in the vicinity of the origin. Only then one may assume that co-existence of two strains within a single node is highly improbable, and therefor treat infections with different species as mutually exclusive events. Starting from (16), from this point onward we consider several scenarios.

3.1 Bounding

From (16), the following relation may be formulated:

[TABLE]

where $\beta_{l}^{\rm sup}=\max_{c}\beta_{l}^{c}$, $l=1..R$ and $\gamma^{\rm sup}=\min_{c}\gamma_{c}$, are the values of $\beta_{l}^{c}$, $l=1..R$ and $\gamma_{c}$ that maximize the right-hand side of the equation (16) for given $i$. Notice that the choice of these parameters, under the given criteria is independent of $i$. We use this fact to rewrite the relation (15) in the following vector form:

[TABLE]

where $\mathbf{T^{\rm sup}}=[t_{ij}]$, is a $n\times n$ matrix, which elements $t_{ij}=\sum_{l\in R}a^{l}_{ij}\beta^{\rm sup}_{l}$ and $\mathbf{\epsilon}=[\epsilon_{1}\;\epsilon_{2}\;..\epsilon_{N}]^{T}$ is a cumulative (all strains considered) epidemic perturbation vector. Now, consider the system:

[TABLE]

The dynamical system (19), is a cumulative variant of the system (7), and in accordance with the relation (18), upper-bounds the behaviour of the dynamical system (3), (4), (5), and (6) near the region of the epidemic origin (i.e. vector $\mathbf{0}$). From (19), it follows that $\mathbf{\epsilon^{\rm sup}}\rightarrow\mathbf{0}$, when $t\rightarrow\infty$, if the larges eigenvalue of the matrix $\mathbf{Q}^{\rm sup}=\mathbf{T^{\rm sup}}-\gamma^{\rm sup}\mathbf{I}$, $\lambda_{1}(\mathbf{Q}^{\rm sup})$, satisfies the relation $\lambda_{1}(\mathbf{Q}^{\rm sup})<0$. Since each eigenvector of the matrix $\mathbf{T^{\rm sup}}$ is also an eigenvector of the matrix $\mathbf{Q}^{\rm sup}$, $\lambda_{1}(\mathbf{Q}^{\rm sup})$, can be expressed as:

[TABLE]

In accordance with the derived relations, one may claim that the point $\mathbf{p}=\mathbf{0}$ is a (locally) stable fixed point of the dynamical system (3), (4), (5), and (6), provided:

[TABLE]

We define:

[TABLE]

as an upper-bound basic reproductive number. From the discussion above, it is clear that for $R^{+}<1$ there are no infected nodes. One should notice that, since derivation of $R^{+}$ was conducted by the analysis of dynamical behaviour of a variant systems that upper-bounds the dynamical behaviour of the system (3), (4), (5), and (6) near the origin, the inverse does not hold, i.e $R^{+}>1$ does not necessarily imply that the epidemics exists.

Identically as before, starting from the equation (16), and by introducing $\epsilon_{i}=\sum_{c}\epsilon^{c}_{i}$, $\gamma^{\rm inf}=\max_{c}\gamma_{c}$, $\beta^{\rm inf}_{l}=\min_{c}\beta^{c}_{l}$, $\mathbf{T^{\rm inf}}=[t_{ij}]$, such that $t_{ij}=\sum_{l\in R}A^{l}_{ij}\beta^{\rm inf}_{l}$ and $\mathbf{Q}^{\rm inf}=\mathbf{T^{\rm inf}}-\gamma^{\rm inf}\mathbf{I_{N\times N}}$, one obtains:

[TABLE]

The dynamical system:

[TABLE]

lower-bounds the behaviour of an arbitrary perturbation of the system (3), (4), (5), and (6) near the epidemic origin, i.e. vector $\mathbf{0}$. Again, the well known result from linear system analysis suggests that $\mathbf{\epsilon^{\rm inf}}\rightarrow\mathbf{0}$, when $t\rightarrow\infty$, if the larges eigenvalue of the matrix $\mathbf{Q}^{\rm inf}$, satisfies the relation $\lambda_{1}(\mathbf{Q}^{\rm inf})<0$. By following similar reasoning as before and considering that in this case the analysis of the system (3), (4), (5), and (6) near the epidemic origin is conducted by analysis of a lower-bounding variant-system, one may claim that infection in the network will surely persist, if the following condition holds:

[TABLE]

By defining the lower-bound basic reproduction number as:

[TABLE]

one may claim that, for $R^{-}>1$ infection surely exists in the observed model described with (3), (4), (5), and (6). Contrary, $R^{-}<1$ does not guarantee that an ineffective process described with (3), (4), (5), and (6) will be eradicated from the network as time elapses. To summarize:

If $R^{+}<1$ an ineffective process described with the system of equation (3), (4), (5), and (6) will surely vanish as time elapses.

• 2.

If $R^{-}>1$ an ineffective process described with the system of equation (3), (4), (5), and (6) will persist in the network indefinitely.

• 3.

If $R^{+}>1$ and $R^{-}<1$ there is no way to ether confirm nor deny the existence of the epidemic in the network. Only numerical computation of $\lambda_{1}(\mathbf{D})$ may give a definite answer in this case.

In the analysis above, extreme cases scenarios are considered: strains either spread from node to node at maximal spreading rate and die-out with the rate of the most sensitive strain, or they spread at minimal spreading rate an die out with the rate of the most endurable strain.

3.2 Rapid mutation processes

As a special case, one may consider the system where mutation occurs on a much faster time scale then the infection/die-out process. There are many reasons why such an assumption should hold in practice. First, it is normal to assume that an entity that has capacity to mutate, is designed (biologically or otherwise) to do so several times before it dies out. Next, processes related to small entities (viruses) occur on a much faster time scale then processes related to the host (the node). Concretely, a virus or a bacteria will self replicate and the new entities will undergo a number of biological processes much faster then hosts immune system responds to the epidemic presence or the host is involved in interaction with other entities (nodes) that introduces spreading (these are, finally, the very basic premises on which epidemic spreading models are built).

For these cases one may assume that once a node is infected, the spreading entity will undergo several mutations before actually spreading to neighboring nodes or dying out. Consider an arbitrary node $i$. Let $\mathbf{\epsilon_{i}}=[\epsilon_{i}^{1}...\epsilon_{i}^{m}]$. Assuming only mutation occurs, within a reasonable time frame $[0,t_{mut}]$, the dynamical behavior of the infection within node $i$, within the time frame $[0,t_{mut}]$, near the epidemic origin, may be described with the equation:

[TABLE]

where $\mathbf{M}=[m_{i,j}]$ is an $R\times R$ matrix, which elements satisfy the relations:

$m_{i,j}=-\sum_{c}\mu_{ic}$ if $i=j$

• 2.

$m_{i,j}=\mu_{ji}$ otherwise.

It is straightforward obvious that all the columns in $\mathbf{M}$ sum up to zero, therefore $\lambda_{1}(\mathbf{M})=0$ and $\mathbf{\epsilon}\rightarrow\mathbf{u}=[u_{0}...u_{m-1}]$, when $t\rightarrow\infty$, where $\mathbf{u}$ is the eigenvector of $\mathbf{M}$ corresponding to $\lambda_{1}(\mathbf{M})$, and ${\sum_{c}u_{c}}=1$.

If now one considers equation (16), and assume that the rate of mutation process by far exceeds the rate of spreading/dying (probability of either spreading or dying of a particular specie, within time frame $[0,t_{mut}]$ significantly low), one may claim that the probability of node $i$ being infected with the strain $c$ near the epidemic origin, may roughly be estimated as:

[TABLE]

which transforms equation (16) into:

[TABLE]

By setting $\beta_{l}^{\rm eff}=\sum_{c}u_{c}\beta^{c}_{l}$, $\gamma^{\rm eff}=\sum_{c}\gamma^{c}u_{c}$, $\mathbf{T^{\rm eff}}=[t^{\rm eff}_{i,j}]$, such that $t^{\rm eff}_{i,j}=\sum_{l\in R}a^{l}_{ij}\beta_{l}^{\rm eff}$, we can re-write equation (24) in the vector form as:

[TABLE]

Using similar arguments as before, we conclude that, provided the rate of mutation exceed the rate of infection/entity extinction, the epidemic threshold of the system (3), (4), (5), and (6) is defined by:

[TABLE]

or in terms of the reproductive number $R^{\rm eff}$:

[TABLE]

4 Linear stability analysis of two-contingent models

Previous section discusses several general properties of our model (3), (4), (5), and (6). Here we derive close formulae for basic reproduction numbers and epidemic thresholds for various two-contingent models. The disease-free equilibrium is stable if the spectral radius of the Jacobian matrix $\lambda_{\rm max}(J)=\max\{\Re(\lambda):\lambda\in\text{spectrum of }J\}$ is less than [math]. Hence, the critical values of the transmission rates $(\beta^{c}_{l})$ can be found by examining when $\lambda_{\rm max}(J)$ traverses 0 as $(\beta^{c}_{l})$ are tuned.

4.1 Models on multilayer networks

We consider the case when the number of contagions is two, while the number of layers is arbitrary. The Jacobian matrix is equal to

[TABLE]

It is useful to partition the matrix $J$ so that

[TABLE]

with

[TABLE]

As we are interested in the case where the maximum eigenvalue of $J$ equals zero, we need to solve $\det J=0$. The determinant of $J$ takes the form

[TABLE]

Let $B=T_{2}-I(\gamma_{2}+\mu_{21})$ and assuming $B$ is invertible, we have

[TABLE]

Because invertible matrices have non-zero determinants ($\det B\neq 0$), the basic reproduction number is derived from $\det(T_{1}-I(\gamma_{1}+\mu_{12})-\mu_{21}IB^{-1}\mu_{12}I)=0$ to be

[TABLE]

where $H=T_{1}-\mu_{12}\mu_{21}B^{-1}$. Equation (29) is a generalization of Eq. (1). Assuming that the transmission rates of contagions 1 and 2 are $\beta^{1}$ and $\beta^{2}$, respectively, on all layers, $R_{0}=1$ in Eq. (29) defines on the $(\beta^{1},\beta^{2})$ plane a critical curve, which is, however, not analytically tractable. Expressions of the basic reproduction number in terms of eigenvalues of the adjacency matrices, and therefore, critical values of the transmission rates can be derived only for two points on the $(\beta^{1},\beta^{2})$ plane: $(\beta^{1}_{\rm cr},0)$ and $(0,\beta^{2}_{\rm cr})$. Indeed, for contagion 1, setting $\beta^{2}=0$ in Equation (28) we have

[TABLE]

Similarly for contagion 2, the following reproduction numbers are acquired

[TABLE]

The critical transmission rates can be derived by assuming that $\lambda_{\rm max}(T_{c})=\beta^{c}\lambda_{\rm max}(\sum_{l\in R}A_{l})$ and $R^{1}_{0}=R^{2}_{0}=1$. Therefore,

[TABLE]

This is a generalization of (2). For the single layer case and $\mu_{12}=0$ it can be seen that (33) reduces to the well known critical value $\beta_{\rm cr}=\gamma/\lambda_{\rm max}(A)$.

Now, we provide numerical simulations on random networks to investigate the correctness of the derived basic reproduction numbers and the critical transmission rates. The numerical simulations are conducted on a two-layer network with $n=50$ nodes where the first layer is an Erdos-Renyi network with parameter $p=0.05$, and the second layer is a Barabasi-Albert network with parameter $k=2$. Figure 3 depicts the sign of the spectral radius of the Jacobian matrix (27) as the transmission rates $\beta^{1}$ and $\beta^{2}$ are ranged from [math] to $0.1$. For every pair $(\beta^{1},\beta^{2})$ the basic reproduction number $R_{0}$ is computed using Equation (29) and the curve where $R_{0}=1$ is drawn. It can be seen that for $R_{0}<1$ the disease-free equilibrium stable and for $R_{0}>1$ it is not. The critical values, Eq. (33) and (34), are computed to be $\beta^{1}_{\rm cr}=0.041$ and $\beta^{2}_{\rm cr}=0.069$ for $\gamma_{1}=0.1,\gamma_{2}=0.3,\mu_{12}=0.9,\mu_{21}=0.3$, and $\beta^{1}_{\rm cr}=0.04$ and $\beta^{2}_{\rm cr}=0.061$ for $\gamma_{1}=0.1,\gamma_{2}=0.3,\mu_{12}=0.5,\mu_{21}=0.3$.

4.2 Competition models on multilayer networks

We now consider the case of our coinfection model when each node can either be susceptible or infected by only a single disease, yielding a competition model. We further assume that each disease can be transmitted through one layer (disease c transmits on layer c), so that $\beta^{c}\equiv\beta^{c}_{l}\neq 0$ only when $c=l$. Since in this case $|L|\leq 1$, the model equations (4) become

[TABLE]

for $c=1,\ldots,m$ and $i=1,\ldots,n$. Similarly to the previous section we are interested in the stability of the disease-free equilibrium of the competition model. The Jacobian matrix takes the form

[TABLE]

where $\bigoplus$ denotes the direct sum of $L$ matrices and $\Gamma$ is the mutation rate matrix

[TABLE]

For the two-contingent case the linear stability analysis of the disease-free equilibrium is similar to that of the coinfection model with two-contingents. The Jacobian takes the form

[TABLE]

the basic reproduction number becomes

[TABLE]

and the critical transmission rates have the following simplified form

[TABLE]

We conduct numerical simulations on two-layer Erdos-Renyi random networks with $n=50$ nodes and edge probability $p=0.1$. Figure 4 depicts the fraction of infected nodes at the steady state, as the transmission rates $\beta^{1}$ and $\beta^{2}$ are ranged from 0 to 0.2. For every pair $(\beta^{1},\beta^{2})$ the basic reproduction number $R_{0}$ is computed using Equation (39) and the curve where $R_{0}=1$ is drawn. It can be seen that for $R_{0}<1$ the fraction of infected nodes is 0, and for $R_{0}>1$, $\sum_{i=1}^{n}p^{1}_{i}+p^{2}_{i}>0$. The critical values, Eq. (40), are computed to be $\beta^{1}_{\rm cr}=0.06$ and $\beta^{2}_{\rm cr}=0.098$ for $\gamma_{1}=0.1,\gamma_{2}=0.3,\mu_{12}=0.2,\mu_{21}=0.2$, and $\beta^{1}_{\rm cr}=0.104$ and $\beta^{2}_{\rm cr}=0.114$ for $\gamma_{1}=0.5,\gamma_{2}=0.5,\mu_{12}=0.2,\mu_{21}=0.2$.

4.3 Compartmental models

Multiple-contagion compartmental model consists of $2^{m}$ compartments. Dynamical equations describing the time evolution of the compartmental model are acquired by removing variables that are network dependent. As an example, we consider 2 contingents (resulting in 4 compartments), $d$ layers, and write $p^{\emptyset}=p_{0}$, $p^{\{1\}}=p_{1}$, $p^{\{2\}}=p_{2}$, and $p^{\{1,2\}}=p_{3}$. Then the model equations can be written as:

[TABLE]

Since $p_{0}(t)+p_{1}(t)+p_{2}(t)+p_{3}(t)=1$ for all $t$, we keep only $p_{1},p_{2},p_{3}$ variables. The Jacobian matrix of the model evaluated at the disease-free equilibrium is a $3\times 3$ matrix:

[TABLE]

One eigenvalue is equal to $-(\gamma^{1}+\gamma^{2})$ and hence is always negative. The other two eigenvalues of the model are the eigenvalues of the matrix $A$ defined as:

[TABLE]

The largest eigenvalue of the Jacobian $J$ is, therefore, given by:

[TABLE]

where $a=\sum_{l=1}^{d}\beta^{1}_{l}-(\gamma^{1}+\mu^{12})$ and $b=\sum_{l=1}^{d}\beta^{2}_{l}-(\gamma^{2}+\mu^{21})$. The last equation can be rewritten as

[TABLE]

where $\epsilon$ is a function of $\mu^{12}$ and $\mu^{21}$. Setting $\lambda_{\rm max}=0$, the basic reproduction number becomes:

[TABLE]

We consider two examples. In the first example we assume $d=1$ and let $\lambda_{\rm max}=0$ (equivalently $R_{0}=1$). Then, the critical values of $\beta_{1}^{1}$ and $\beta_{2}^{2}$, $\beta^{1}_{cr}$ and $\beta^{2}_{cr}$, are related to each other through equation $ab=\mu^{12}\mu^{21}$, which can be rewritten as

[TABLE]

for $0\leq\beta^{1}_{cr}\leq c_{1}$ where $c_{1}>0$ is a constant that depends on the parameters $\mu^{12}$, $\mu^{21}$, $\gamma^{1}$, and $\gamma^{2}$. Assume now that $\gamma^{1}=0$, then for $\beta^{1}_{cr}=0$ and $\beta^{2}_{cr}=\gamma^{2}$, $\lambda_{\rm max}=0$, while for $\beta^{1}=0$ and $\beta^{2}>\gamma^{2}$, $\lambda_{\rm max}>0$. Therefore, even when $\beta^{1}=0$, the epidemic state is present in the model (when $\beta^{2}>\gamma^{2}$). This is illustrated in Figure 6. We fix the parameters to $\gamma^{1}=0.8$, $\gamma^{2}=0.2$, $\mu^{12}=0.3$, $\mu^{21}=0.2$, and plot $\beta_{cr}^{2}$ vs $\beta^{1}_{cr}$ on Figure 5. For $\beta^{1}_{cr}=0$, we have $\beta_{cr}^{2}=0.345$. Figure 6 shows how the size of the infected population evolves over time for $\beta^{1}=0$ and $\beta^{2}=0.352$: from a population of 400 approximately 7.142 of the population are infected with infection 2, 1.293 with infection 1, 0.005 with both infections 1 and 2, and 391.56 are susceptible. For any $\beta^{2}>\beta^{2}_{cr}$, the model approaches a fixed point solution different than the origin, although the values of $p^{1}$, $p^{2}$, and $p^{12}$ are small (of order of $\beta^{2}-\beta^{2}_{cr}$).

For the second example, we assume that $d=2$, and write $\beta=\beta^{1}_{1}=\beta^{2}_{2}$, and $\hat{\beta}=\beta^{1}_{2}=\beta^{2}_{1}$. Setting $\lambda_{\rm max}=0$ (equivalently $R_{0}=1$) leads to the following two relations between the critical values of $\beta$ and $\hat{\beta}$, $\beta_{cr}$ and $\hat{\beta}_{cr}$,

[TABLE]

which, after some algebra, can be written as $\hat{\beta}_{cr}=\beta_{cr}+c$, for $0\leq\beta^{1}_{cr}\leq c_{2}$ where $c_{2}>0$ and $c>0$ are constants that depend on the parameters $\mu^{12}$, $\mu^{21}$, $\gamma^{1}$, and $\gamma^{2}$. Again the absence of epidemic threshold ($\beta=0$) is observed while the epidemic state is presented (when $\hat{\beta}>c$).

Finally, we remark that in a special case when $\beta_{l}^{1}=\beta_{l}^{2}\equiv\beta_{l}$ for all $l$ and $\gamma_{1}=\gamma_{2}\equiv\gamma$, the basic reproduction number reduces to:

[TABLE]

for arbitrary values of $\mu^{12}$ and $\mu^{21}$. Therefore, in this example, mutation does not influence the basic reproduction number and the corresponding threshold as long as the transmission and recovery rates of two contingents are equal to each other. Moreover, if, in addition, $\beta_{l}=\beta$ for all $l$, the basic reproduction number reduces to

[TABLE]

5 Conclusion

In summary, we have suggested a multiple contagion SIS model on a multilayer network that incorporates different spreading channels and disease mutations. The model is analytically tractable; in particular, the absence of epidemic threshold and critical behavior even for compartmental models could further improve our understanding of epidemic spreading. The results obtained here could have implications in various disciplines including epidemiology, biology, and sociology.

Acknowledgment

This research was supported in part by DFG (grant Random search processes, Lévy flights, and random walks on complex networks) and ONR/ONR Global (grant No. N62909-16-1-2222).