Endemic state equivalence between non-Markovian SEIS and Markovian SIS model in complex networks.

In the light of several major epidemic events that emerged in the past two decades, and emphasized by the COVID-19 pandemics, the non-Markovian spreading models occurring on complex networks gained significant attention from the scientific community. Following this interest, in this article, we explore the relations that exist between the mean-field approximated non-Markovian SEIS (Susceptible-Exposed-Infectious-Susceptible) and the classical Markovian SIS, as basic reoccurring virus spreading models in complex networks. We investigate the similarities and seek for equivalences both for the discrete-time and the continuous-time forms. First, we formally introduce the continuous-time non-Markovian SEIS model, and derive the epidemic threshold in a strict mathematical procedure. Then we present the main result of the paper that, providing certain relations between process parameters hold, the stationary-state solutions of the status probabilities in the non-Markovian SEIS may be found from the stationary state probabilities of the Markovian SIS model. This result has a two-fold significance. First, it simplifies the computational complexity of the non-Markovian model in practical applications, where only the stationary distributions of the state probabilities are required. Next, it defines the epidemic threshold of the non-Markovian SEIS model, without the necessity of a thrall mathematical analysis. We present this result both in analytical form, and confirm the result through numerical simulations. Furthermore, as of secondary importance, in an analytical procedure we show that each Markovian SIS may be represented as non-Markovian SEIS model.

Introduction

Following the outbreak of several recent epidemics, the SARS, the MERS, Bird flu etc., and emphasized by the Covid-19 pandemics, the non-Markovian models captured the attention of the complex networks research community that studies stochastic spreading processes [1], [2], [3], [4], [5], [6], [7], [8], [9], [10]. The shift of interest from the Markovian to the non-Markovian realm was caused by the realization that no status transition, that an individual (node) undergoes following the contraction of the spread agent, may neither occur simultaneously, nor the probability of status transition on a daily bases is constant, as may be seen from the collected medical data [11], [12]. For example, when an individual contracts a virus, the individual has no capacity to instantly spread the agent to the neighboring nodes. A certain viral quantity should first be produced (the node acts as an incubator for the agent), in order for the infected individual to become infectious. The time required for the virus to replicate over the necessary viral load threshold, varies from individual to individual and follows a time dependent probability distribution [13]. On the other hand, the recovery process is a product of complex biochemical interactions within the hosts immune system, that take several days in order for proper defence response to be prepared and for the virus to be eradicated (as an example of the immune system modeling one may refer to [14]). Again, individual differences lead to specific probability distribution of time from Exposure to Recovery.

Markovian models do not consider these features. Markovian dynamics of spreading assumes that, the node that contracts the virus may immediately, or within the next time instance, spread the agent further to its neighbors and/or recover. Probabilities of transitions to subsequent statuses, following the agent contractions, are constant within equal time intervals, and independent of the time when contraction occurred. No memory exists regarding the duration of the period from exposure to the moment of subsequent transition.

In our recent paper [15], we introduced the mean-field discrete time non-Markovian SEIS (Susceptible–Exposed–Infectious–Susceptible) model as a basic mathematical non-Markovian form that describes re-occurring spreading processes, taking place on complex networks. In the model formulation, We assumed that status transitions from Exposed (non-Infectious) to Infectious status and Exposed (both non-Infectious and Infectious) back to Susceptible status, follow temporal distribution described with Discrete Time Probability Functions (DTPFs):

daily manifesting function : probability that an Exposed and previously non-Infectious node, becomes Infectious exactly at day ;

manifesting function : probability that an Exposed node, is Infectious at day ;

daily recovering function : probability that an Exposed node, recovers exactly at day ;

recovering function : probability that an Exposed node, is recovered by day ;

In this paper, We first extend the discrete-time concept to continuous-time non-Markovian model form. Adequately, the functions , , and in this scenario are continuous, and are further referred to as Continuous-Time Probability Functions (CTPFs), with and referred to as instance recovering probability and instance manifesting probability, correspondingly. For the continuous-time form, we derive the epidemic threshold in a strict mathematical procedure. Then, the main result of this paper is presented: that for each mean-field non-Markovian SEIS model (discrete-time or continuous-time), exists a Markovian SIS model, such that the stationary state probabilities of each node being exposed in the SEIS model, equals the stationary probability that the node is Infected in the SIS model, providing certain relations between process parameters hold. We consider this result to be of at-most importance for the following reason: non-Markovian models, although highly accurate in analyzing natural phenomena, are computationally sufficiently more demanding. Investigating these models with utilization of Markovian analogs (as shown here as possible), significantly reduces the computational complexity in acquiring significant data related to the endemic state of the diseases. The presented analysis directly leads to relations that define the epidemic threshold for the non-Markovian (SEIS) models occurring on complex networks, without the necessity of a thrall mathematical procedure. Similar type of equivalences, between the non-Markovian and the Markovian SIS model, using different settings and approaches, have been established by the authors in [1], [10]

As a result of secondary importance, it is shown that an arbitrary Markovian SIS model occurring on complex networks, may be represented an non-Markovian SEIS model. This equivalence is only vaguely mentioned and numerically illustrated in the Conclusions of [15]; here we show this feature through a rigorous mathematical procedure. One should note that similar analysis was conducted in respect to the non-Markovian and Markovian SIR model in [16].

The models

In this Section, we present a short introduction to the processes and models used in the analysis. First, a very brief description of the well known mean-field Markovian SIS model is presented. Then, following [15], a more detailed description of the discrete-time non-Markovian SEIS process and the accompanied mean-field model is given. Finally, a mean-field approximated continuous-time SEIS model is introduced, and for this model the epidemic threshold is derived.

In this paper, it is considered that all processes occur on a network represented with the adjacency matrix . In the general case, the network is directed, weighted, and strongly connected; consequently the matrix is asymmetric, with , and irreducible (Perron–Frobenius theorem for non-negative irreducible matrices holds).

Markovian SIS model

Although well known and widely elaborated in the scientific literature, in this sub-section we give a short introduction to the Markovian SIS model occurring on complex networks.

The Markovian SIS model is a status model, in which each node, at given moment, may be in one of two distinctive statuses: Susceptible (S) and Infected (I). A node is infected at time , if the node was Susceptible at time , and acquired the spread agent from its Infected neighbors in the time interval , or was Infected at time and did not recover in the interval . Probability of the recovery, within a unity time interval is defined by the curing rate . Probability of acquiring the spread agent by a Susceptible node , from one of its infected neighbors is defined with the product , where is the averaged epidemiological significance of contact between nodes, acting in the direction, within a unity time interval, and is the infection rate, i.e. the probability of infection following a certain epidemiological contact with infected neighbor.

Following the formal definition stated above, and assuming statistical independence of joint events, the mean-field Markovian SIS model in discrete-time is described with [17], [18]: where is the probability that node is infected at time .

This equation form of the Markovian SIS model is often referred to as microscopic Markov chain approach in the literature [19], [20], [21]. It is the most common form of representation of the SIS process in the study of epidemic diseases.

The Markovian SIS model in continuous form is also well known [22], and may be written as:

Principles on which the mean-field continuous-time Markovian SIS model is derived are similar to those applied in formulation of the microscopic Markov chain approach, i.e. the assumption of statistical independence of joint events. However, one major difference exist between both formats: due to infinitesimal duration of the model updating period, in the continuous-time model form, it is fairly assumed that a susceptible node may acquire the spreading material from only one of its infected neighbors. This notion transforms the product-like term, that approximates the probability of transition of spread material from infected neighbors to a Susceptible node, into a sum-like term. The presented approach of mean-field modeling the Markovian SIS process is referred to as Quenched Mean-Field (QMF) or N-Intertwined Mean-Field Approximation (NIMFA), in the scientific literature [19].

Other methods of modeling the SIS process, have been studied in the past 20 years, as well. Early works on the subject, revolved around the degree-based mean-field (DBMF) approximation, introduced by Pastor-Satorras and Vespignani in [23], [24]. Recently, models that take into account the second-order statistical dependence, and therefore significantly improve on the accuracy, but increase the complexity of the mathematical models, were introduced in [25], [26].

The non-Markovian SEIS model

The discrete-time form of the SEIS model analyzed in this paper is originally introduced in [15]. For completeness, in what follows, we re-state the formal definition of the observed process, and for more details we refer the readers to the cited paper.

The SEIS model, as described in [15] and in this work, is a status model, which, in respect to the spreading agent, is characterized by two disjunctive statuses, Susceptible (S) and Exposed(E), and by status Infectious (I), that is a sub-state of the E-status. A node is in status Exposed, at time , if at the given instance it contains the spread agent; otherwise the node is considered to be Susceptible. Exposed node may be Infectious (manifesting infectiousness) or non-Infectious. Node is Infectious if it contains the agent (is Exposed) and is capable to spread the agent further to its neighbors.

Susceptible node contracts the spread agent from one of its Infectious neighbors and becomes Exposed at time , depending on two factors. First is the epidemiological importance (on average) of the contact at the given instance, denoted with the probability . Second, is a cluster of unknown variables that determine whether the transfer of spread agent will occur, following a certain epidemiological contact, with potential to evolve into an infectious process within the contracting node, and are described by the well known parameter , i.e. the infection rate. In this paper, as well as in [15], it is considered that factors described by the parameter are purely stochastic and by no means may be referenced to events related to the spreading dynamics; in that sense is considered constant throughout the process. As an example, one may consider to represent the fraction of the individuals that have no immunity, or are prone to developing an infection towards certain viral disease; average capacity of an individual to produce the viral load required for infecting others e.t.c. In a number of recent papers, that explore forms of the continuous-time non-Markovian SIS model [1], [4], [5], [6], [10], parameter is time dependent (p.d.f) and denotes the probability of the Infected node to transfer the spread agent to its neighbors at given moment. This role in our model is played by DTPFs/CTPFs and .

Similar discussion is relevant in respect to parameters : knowledge about potential infectious contact, followed by occurrence of symptoms, alters the contact dynamics of the individual (the spreading node). However, implementing these behavioral changes in the model, by far exceed the scope and will increase the complexity of this paper; for those reasons, in what follows, we consider parameters to be constant, as well. In this sense, the model evolves around the notion that the transition is purely Markovian in character. It is the authors impression that, disregarding behavioral changes of the Exposed/Infectious nodes, this setting adequately mimics the reality.

The process of agent contraction by the Susceptible node plays a role of a trigger event (): all consequent processes within the node are time-referenced to this transition. Exposed (but non-Infectious node) may become Infectious exactly at time after the trigger event with probability . Exposed node is Infectious at time following the trigger event with probability . To stress the difference between and , as explained in [15], the model allows for two different types of Infectiousness manifestation:

Cumulative manifestation — in this case has a character of a mass probability function in the discrete-time case scenario and density probability function in the continuous-time scenario. Adequately, (discrete-time), with , or (continuous-time), with , has a cumulative character, with the sign “” indicating that the Exposed node may not necessarily become Infectious prior to recovery. This type of behavior is typical for epidemic diseases;

Random manifestation — in this case has a random character, with being the only restriction.

Exposed node recovers and becomes Susceptible again exactly at time following the exposure, with probability , and is recovered at time with probability , with , in the discrete-time model, and , with in the continuous case. The probability function (DTPF/CTPF) is a m.p.f in the discrete-time scenario, and p.d.f. in the continuous-time case, with being a cumulative probability function. In the modeling of the SEIS process, we widely use the complement .

In the formal sense, the dynamical behavior of the model is defined as follows: node is Exposed at time if it contracted the agent at time and did not recover in the time interval ; node is Infectious at time if it contracted the agent at time , did not recover in the time interval and is capable to spread the agent to its neighbors at time . We consider that the process lasts for maximum time units, with , in the discrete-time scenario, and , in the continuous-time version.

Let denote the probability that node is Exposed, and denote the probability that node is Infectious, at time . Considering the definitions stated above, and following [15], the mean-field discrete-time SEIS model is mathematically defined in the following form: with representing the product-like term: that denotes the probability that a Susceptible node will contract the spread agent from its neighbors, at time . The system of Eqs. (3) is derived relying on the well known assumption of statistical independence of joint events, used as basis in developing microscopic Markov chain models [17], [18], [27], [28].

For details related to derivation of Eq. (3), we refer the readers to [15].

One should note the difference in the formulation of the non-Markovian SEIS model, as introduced in [15] and presented in this paper, compared to the classical Markovian SEIS model. In the classical version the E-stage refers only to nodes that contracted the spread agent, but are non-Infectious [29], [30], [31]. This distinguishes between E and I statuses as disjunctive, leading to conservation equation . In the non-Markovian SEIS model, the I status is a sub-state of the E status, consequently the conservation equation is of form . The reason for this lays in the character of the processes. In the Markovian scenario, the probability of each of the transitions , , and is constant at every time instance, regardless of the moment of initiation of the cycle , i.e. the time when transition, occurred. In the non-Markovian scenario presented here, the probability of occurrence of each subsequent transition, following the transition depends on the moment in time when the cycle begun. Treating status E as a compound of sub-statuses (in the SEIS model only sub-status I), creates a unique time-frame in which all other processes that an individual node undergoes, following the contraction of the spread agent, occur. In that sense, one should note that separation of status Exposed (but non-Infectious) and Infectious, as well as other statuses that may occur in extended non-Markovian scenarios (models), into disjunctive statuses is possible (with adequate adjustments to the DTPFs/CTPFs); however that requires introduction of multiple time-frames (each subsequent status should be referenced to the previous). Even further, each of these time-frames should be referenced to the origin of the initial time-frame (i.e. the moment of the transition). This will make the analytical model more complicated.

Continuous-time SEIS model

In this section, we introduce the continuous-time non-Markovian SEIS model, as an extension to the discrete-time model (3). Prior to the formal introduction of the model, we stress the major difference that exist between discrete-time and continuous-time modeling approach. It is a standard practice in modeling spreading phenomena in continuous time to assume that the transfer of the spread agent, within an infinitesimal time interval , may occur from a single sources (neighbor), only (no-multiple infectious events assumption). This notion transforms the product-like term in the following manner: For sufficiently small , the term is neglected, and the appropriate sum-like term [22], obtained.

Bearing in mind the differences, We may now re-write the system of Eqs. (3) as follows:

When , providing no discontinuities of first kind exist in or , the terms multiplied by may be neglected. In what follows, we show that this term may be neglected even in the presence of finite number of first order discontinuities. The analyses is focused around the second set of equation in the system (4), related to the , variables; by analogy, the same analysis is valid for the set of equations related to the variables.

Consider a point , such that a discontinuity of first kind or exists at . Let be an integration constant, such that the series provides a proper sampling of and . Let be an index such that . Under these assumptions, the following inequality may be considered:

Let R be a total number of first kind discontinuities of either and . Then, in accordance with the relation above:

The preceding analyses indicates that the terms multiplied by in the set of equations, related to the in (4), may be neglected, since for finite , may be maid arbitrary small, with the right choice of . Similar analysis, leading to the same conclusion, may be conducted for the set of equations related to in (4). Consequently, from (4) and considering , one obtains the integral form of equations for the non-Markovian SEIS model occurring on complex networks in continuous-time: with being the Heaviside function. In what follows we consider both and to be smooth around the point , and the Heaviside function may be neglected in the system of Eqs. (5).

Differential form

In this segment, we show that the non-Markovian SEIS model, represented with (5), may be written in a differential form, as well. The purpose of this model-form is to relate the non-Markovian SEIS model and the Markovian SIS model, in order to investigate the circumstances under which an arbitrary Markovian SIS may be presented as non-Markovian SEIS.

Starting from the system of Eqs. (4), one obtains:

Considering that , by dividing both sides of the equation with and by letting , the following relation may be written :

Similarly for one obtains:

Epidemic threshold for the continuous-time model

One of the main results that are derived in the theoretical analysis of the stochastic spreading processes occurring on complex networks, is the determination of the epidemic threshold. Epidemic threshold defines the critical relation between the process parameters and network topology, that separate the parametric region in which the network is disease free, from the region in which a permanent epidemic exists.

To find the epidemic threshold for the continuous-time non-Markovian SEIS model, we resort to the investigation of the stability criteria of the dynamical system (5), around the point of epidemic origin, i.e. , , for all . In that sense, we consider the following:

Consider a directed, weighted and strongly connected graph, represented with the adjacency matrix that is, consequently, non-negative and irreducible. The vector , i.e. the epidemic origin, is a globally asymptotically stable point of equilibrium of the dynamical system (5) , providing the following relation holds: with being the leading eigenvalue of the matrix .

Consider the system (5). Since , for all , and , the argument under the integral is strictly positive, so the following relation holds: In other words, the dynamical behavior of the system (5) is bounded from below by the epidemic origin, and from above, by the dynamical system:

Consequently, if , , then , , as well. In that sense, the proof of the global stability of the (epidemic origin of the) system (5), reduces to proof of the global stability of the system (8).

One may notice that the second set of equations in (8) is self-sufficient. In that sense the dynamical stability of the system (8) reduces to the dynamical stability of this set of equations only: from (8) follows that, if , then .

By conducting a Laplace transform on both sides of each of the equations from the second set of equations in (8), following the methodology in [16], one obtains: with being the Laplace transform of the product , and the term that takes into account the initial conditions. The initial conditions and the nature of this term are discussed separately in Remark 1, at the end of this segment.

Consider the vector and . The system of Eq. (9) may be re-written in a vector form as: leading to a solution in Laplace domain, in the following form: where is a matrix, which elements are the minors of the matrix .

It is a well known result in the dynamical system theory, that the stability of the (origin of the) dynamical system is determined by the position of the poles of the system in the complex plane. If all poles of the dynamical system lie within the left-half of the complex plane, i.e. , the dynamical system is globally asymptotically stable. From Eq. (12), one obtains that the poles of the system (8) may be determined from the zeros of the equation:

On the other hand: with , , being the eigenvalues of the adjacency matrix . From the last relation, the position of the poles of the dynamical system (8) are determined from the set of equations:

Let , be a pole of the dynamical system (8), associated with the th eigenvalue in the following manner: Index allows for multiple poles associated with a single eigenvalue .

From the definition of the Laplace transform of , i.e. Eq. (10), the following conclusions hold:

if , i.e. if a pole of the dynamical system (8) lies on the right-half of the complex plane, the value of the term lies within the circle in the complex plain;

if , i.e. if a pole of the dynamical system (8) lies on the left-half of the complex plane, the value of the term lies outside the circle ;

Bearing in mind the preceding discussion and the relation (13), the poles of the dynamical system (8) will lie within the left-half of the complex plane, i.e. , providing: for all and .

From the Perron–Frobenius theorem for non-negative and irreducible matrices, the leading eigenvalue, , of the matrix is distinct, real and largest by module, compared to all other eigenvalues; therefore it minimizes the term . For this reasons, providing: holds, the poles of the second set of equations of the dynamical system (8) lie within the left-half of the complex plane, resulting in . This yields , , , and the point of epidemic origin of the dynamical system (5) is globally asymptotically stable.

The Proof is completed. □

In accordance with the Theorem, the relation: defines the boundary between the parametric region related to the state of permanent epidemic presence in the network and the region of epidemic absence. In that sense, Eq. (15) represents the epidemic threshold for the non-Markovian SEIS model occurring on complex networks.

The immense importance of proper inclusion of the initial conditions in the non-Markovian SEIS model, is discussed, for the discrete-time case, in details in [15]. The focus of this article is set around the model analysis in the endemic state, i.e. for circumstances in which initial conditions play a minor role. Therefore, it is fairly assumed that the initial epidemic outbreak occurred at moment preceding the beginning of analysis (time ), for both systems (3), (5) and that properly collected set of initial condition, for the time period , exists.

Equivalence of the stationary states of the non-Markovian SEIS and Markovian SIS model

In this section, we present the main result of the paper — that the stationary state solutions of the mean field non-Markovian SEIS model, may be found from the stationary state solutions of the mean-field Markovian SIS model, providing certain relations between process parameters – infection rates for both Markovian SIS and non-Markovian SEIS model , curing rate for the Markovian SIS model and the DTPFs (CTPFs) and for the non-Markovian SEIS model – hold. The purpose of the following analysis is to show that the steady state probabilities that the node is Exposed, , for the non-Markovian SEIS model, may be directly related to the steady state probability of the node being Infected for the Markovian SIS model. Knowing steady-state value of , one may easily calculate the steady state probability of the node being Infectious, , for the non-Markovian SEIS model.

In the following text, we would frequently re-direct the attention of the reader between the Markovian SIS and the non-Markovian SEIS model. In order to avoid any confusion in respect to the form we are referring to, in this and in the next Section, labels for the Markovian SIS and for the non-Markovian SEIS form would be used for the variables and the parameters (except for and ), in the form of superscripts.

Discrete-time model

Starting from the system of Eqs. (3), the stationary state solutions of the discrete-time non-Markovian SEIS model, , , may be found from the relations:

By dividing equations in (16): One should note that relation (17) is not a consequence of the modeling, rather a general feature of the SEIS process, as shown in the Appendix, Eq. (A.4).

From Eqs. (16), (17), the stationary state solution of the non-Markovian SEIS model, in respect to the variable , may be written in the following form: with defined with:

The stationary state solution of the discrete-time Markovian SIS model, may be found from the system of Eqs. (1), when , and is of form:

We formulate the stationary state equivalence between the Markovian SIS, and non-Markovian SEIS model in a sense that . By taking: and then dividing Eqs. (18), (20), assuming equality holds, one obtains: The analysis conducted above, leads to the following conclusion: providing relations (21), (22) hold, the stationary probabilities of each node being Exposed in the discrete non-Markovian mean-field SEIS model equals the stationary probabilities of the nodes being in status Infected, in the discrete Markovian mean-field SIS model.

In the stationary state, , , the relation (5) takes the following form:

The last pair of equations represent nonlinear system from which one can determine the stationary probabilities. By dividing equations in (23), We show that in the stationary state, the variables and are related in the following fashion: We would again draw the attention of the reader to the fact that, identically as in the case of relation (17), the Eq. (24) is a general feature of the SEIS process, independent of the modeling procedure, as shown in Appendix.

If one substitutes the relation (24) into the first equation of the system (23), one obtains:

From (2), the stationary state solution, , of the Markovian SIS model may be written in the form:

We seek equivalence between the models by equalizing the stationary probability of the arbitrary node being Infected in the Markovian SIS model with the stationary probability of the corresponding node being Exposed in the non-Markovian SEIS model, i.e. . From relations (25), (26) (dividing the two equations under the equality assumption), one obtains that the stationary state of the non-Markovian SEIS spreading process may be obtained form a Markovian SIS process providing following relation holds:

One should note the interesting difference between the relations (21), (22), (27). While (21), (22), in the discrete-time case, fully define both the and for the Markovian SIS equivalent, relation (27), in the continuous-time case, leaves certain degree of freedom in the choice of one of these parameters.

As an interesting consequence, one should note that from the second relationship in (23), when the epidemic is weak , one has

In the matrix form, if is the vector of probabilities of infectious state, and is the network connectivity matrix, one has

The last expression indicates that the Infectiousness probability vector is eigenvector of scaled connectivity matrix in weak epidemic. This is similar to the previous result that the principal eigenvector determines the probabilities of infection in SEAIR model [32].

Determination of the epidemic threshold from the endemic state equivalence

Standard approach in determining the epidemic threshold for stochastic spreading processes, requires thrall mathematical procedure that is based either on establishing the stability criteria for the system around the point of epidemic origin, as done here in Section 2.3.2, or complex statistical analysis. In this segment, we show that the endemic model equivalence, enables us to derive the epidemic threshold for non-Markovian SEIS model, directly from the well-known epidemic threshold for the Markovian SIS model.

The epidemic threshold for the Markovian SIS model is defined as: with being the largest (leading) eigenvalue of the adjacency matrix . This relation holds for both discrete-time SIS model [17], [18], as well as continuous-time SIS model [22], [33].

Stationary-state equivalence between the Markovian SIS and non-Markovian SEIS implies that, providing Eqs. (21), (22) in the discrete-time, and (27) in continuous-time SEIS model hold, then: Consequently, if , then , as well. Since, in the general case, holds for all if the system is parameterized “under” the epidemic threshold, in accordance with relations (21),(22), (27), the epidemic threshold for the non-Markovian SEIS model is defined with: in the discrete-time case, and: in the continuous case. For the discrete-time case the relation (29) is derived in precise analytical procedure in [15]. The relation (30) is identical to Eq. (15) derived in Section 2.3.2.

As presented in this subsection, the stationary-state model equivalence, leads directly to the result for the epidemic threshold of the non-Markovian model. No thrall statistical or system stability analysis is required to obtain this result. This feature further emphasizes the importance of the main result of the paper.

Numerical analysis

In this section we present the results of the numerical analysis, in order to validate the theoretical results obtained in the previous sections. The analysis is focused around the main result of the paper, i.e. the stationary state equivalence between the mean-field non-Markovian SEIS and the Markovian SIS model. In addition, following the formal introduction of the continuous-time non-Markovian SEIS model, and the derivation of the theoretical epidemic threshold (14), we present a short analysis of how the theoretical value compares to reality, based on an extended Monte Carlo simulation.

The analyses in the paper are conducted on the following networks:

Barabási–Albert [34] directed and weighted graph with nodes, total of uni-directional links, and the largest eigenvalue of the graph’s adjacency matrix . The reference BA1000, will be used for this network thought the paper. The network is derived from a symmetrical BA(1000,1996) graph, with nodes, generated with parameters , ;

Watts–Strogatz [35] directed and weighted graph with nodes, total of uni-directional links, and the largest eigenvalue of the graphs adjacency matrix . This network would be further referenced as WS1000. The network is derived from a symmetrical WS(1000,3000) graph, with nodes, generated with parameters , ;

In order to confirm the results for the discrete-time model, following [15], we consider three different sets of DTPFs. The first two sets, labeled CM1 and CM2, are related to the case of cumulative-like manifestation of Infectiousness. In the third set, labeled as RM, the manifestation has random character. The sets are presented in Table 1.

Table 1

DTPFs of the three discrete-time analyzed cases: Cumulative manifestations 1 and 2 and Random manifestation.

	τ	0	1	2	3	4	5	6	7	8	9	10
Cumulative manifestation 1 (CM1)	B(τ)	0.0	0.0	0.0	0.1	0.2	0.3	0.5	0.8	1.0	1.0	1.0
	Γ(τ)	0.0	0.0	0.0	0.0	0.1	0.2	0.3	0.5	0.8	0.9	1.0
Cumulative manifestation 2 (CM2)	B(τ)	0.1	0.2	0.3	0.5	0.8	1.0	1.0	1.0	1.0	1.0	1.0
	Γ(τ)	0.0	0.0	0.0	0.0	0.1	0.2	0.3	0.5	0.8	0.9	1.0
Random manifestation (RM)	B(τ)	0.0	1.0	1.0	0.0	0.0	1.0	1.0	0.0	0.0	1.0	1.0
	Γ(τ)	0.0	0.0	0.0	0.0	0.1	0.2	0.3	0.5	0.8	0.9	1.0

The numerical analysis consists of two sub-parts: parametric dependence analysis and time series analysis. For both sub-parts, in order to compare the behavior of the mean-field models, as well as to compare the models to reality, two types of simulations were conducted: model simulations and stochastic simulations.

In the parametric dependence analysis, the parametric region was swiped with a step of 0.01. Using relations (21), (22), parameter , varied in the parametric region , in the CM1 case, , in the CM2 case and in the RM case, with in all cases. For each value of and the corresponding values of and , the number of the Exposed/Infected/Infectious nodes in the stationary state for both models and processes (stochastic simulations) were compared to each other as a function of . In the time series analysis, for given a value of the parameter and the corresponding and , outcome of both the discrete-time non-Markovian SEIS and Markovian SIS, as well as the stochastic simulations of the processes described by the models, were compared as function of time, including the transient period from the moment of initial infection, until full stationarity is achieved. In addition, both in parametric dependence and in the time series analysis, a synthetic variable is included, as a number of Infectious individuals of the non-Markovian SEIS process, calculated from the number of Infected individuals for the SIS process, assuming process equivalence and using relation (A.4).

Model simulations were conducted using directly Eqs. (3) for the non-Markovian SEIS and (1) for the Markovian SIS model.

The stochastic simulations of the discrete-time SEIS model were conducted as follows: a Susceptible node becomes Exposed, when in contact with infected neighbor , as a result of an outcome of a Bernoulli random generator with probability . The moment the node becomes Exposed, timer is activated, recovery time is set by a random generator with m.p.f. , and in the cumulative case a (first) moment of manifestation is set by a random generator with m.p.f. . The node remains Exposed until the chosen recovery time is reached, while being Infectious from the moment of (first) manifestation. In the random case, at each time step Exposed node may become Infectious as a result of a Bernoulli random generator with probability . In the parametric dependence analysis, the process is iterated for 50k cycles, and the obtained values for the number of Exposed and Infectious nodes, at each time step, averaged in the last 5k cycles. In the time series analysis a classical Monte Carlo simulation was conducted (average over 100 individual trials).

Stochastic simulations on the SIS model, were conducted, following classical methodology: a Susceptible node becomes Exposed, when in contact with infected neighbor , as a result of an outcome of a Bernoulli random generator with probability . Infected node recovers, within a unity time interval, as a result of a decision made by a Bernoulli random generator, with probability . In the parametric dependence analysis, the process is iterated for 20k cycles, with the last 5k values for the number of Infected nodes, at each time step, being averaged. In the time series analysis a classical Monte Carlo simulation was conducted (average over 100 individual trials).

Results of the parametric dependence analysis for the discrete-time case are presented in Fig. 1, with corresponding time series analysis presented in Fig. 2. Analysis of the epidemic threshold for the different DTPFs and graphs is summarized in Table 2.

Fig. 1

Comparison between the steady-state values of the number of Exposed/Infectious/Infected individuals in the discrete-time non-Markovian SEIS model and the Markovian SIS model. Meaning of the symbols in the legend: nM — non-Markovian, M — Markovian; E — Exposed (number of), I — Infected (SIS)/Infectious(SEIS); calculated number of Infectious individuals in the SEIS process, from number of Infected individuals in the SIS process, assuming process equivalence; SS — stochastic simulations; the critical value of that defines the Epidemic threshold (refer to Table 2 for details): (a) BA1000 graph, CM1; (b) BA1000 graph, CM2; (c) WS1000 graph, CM1; (d) WS1000 graph, CM2; (e) BA1000 graph, RM; (f) WS1000 graph, RM.

Fig. 2

Comparison between the typical time series of the number of Exposed/Infectious/Infected individuals in the discrete-time non-Markovian SEIS model and the Markovian SIS model. Meaning of the symbols in the legend: nM — non-Markovian, M — Markovian; E — Exposed (number of), I — Infected (SIS)/Infectious(SEIS); calculated number of Infectious individuals in the SEIS process, from number of Infected individuals in the SIS process, assuming process equivalence; SS — stochastic simulations: (a) BA1000 graph, CM1, (); (b) BA1000 graph, CM2, (); (c) WS1000 graph, CM1, (); (d) WS1000 graph, CM2, (); (e) BA1000 graph, RM, (); (f) WS1000 graph, RM, (). In all cases .

Table 2

Analysis of the epidemic threshold for the three analyzed discrete-time cases and two different graphs: — critical value of the parameter , calculated directly from Eq. (29), and following [15]; — critical value of the parameter , calculated from the value of and Eq. (21); — calculated from Eq. (22).

DTPF set/Graph	β^NM,c	Corresponding β^M,c, γ^M	γ^M/β^M,c Ratio
CM1/BA1000	0.120355	0.026244, 0.138889	5.2922
CM2/BA1000	0.045863	0.026244, 0.138889	5.2922
RM/BA1000	0.052488	0.026244, 0.138889	5.2922
CM1/WS1000	0.194785	0.042474, 0.138889	3.26998
CM2/WS1000	0.074227	0.042474, 0.138889	3.26998
RM/WS1000	0.084948	0.042474, 0.138889	3.26998

The results from the analysis presented in Fig. 1 indicate that there is a perfect overlap of stationary state solutions of both the discrete-time non-Markovian SEIS model and the Markovian SIS model, providing relations (21),(22) hold, in all simulated cases. Also, there is a good match between the (mean-field) models result and the reality (stochastic simulations), as well as between the real (stochastic) SIS and SEIS processes, in all, but for the case of the CM1 DTPFs, Fig. 1, Fig. 1, for the Exposed nodes. These findings are further supported by the results presented in Fig. 2. The results imply that the (mean-field) model equivalence indicates a process equivalence as well, except in regions where there is a certain model error. The nature of this error is investigated, in details, in the Section 4.3.

The threshold analysis presented in Table 2 indicate that there is a perfect match between the epidemic threshold for the non-Markovian SEIS model and the derived Markovian SIS equivalent. Notice, from the table, that the value of is identical in all six cases, since the DTPF is identical for CM1, CM2 and RM, and in accordance with relation (22) is specific. The ratio (fourth column in Table 2), when is calculated from Eq. (21) for (third column in Table 2), equals (with small truncation error) the largest eigenvalue of the analyzed graph — i.e satisfies the epidemic threshold relation for the SIS model, Eq. (28). Equality of the both thresholds may be observed from Fig. 1, as well (the critical value of is represented with the dashed red line).

Continuous-time model

The CTPF’s used in the continuous-time case (cumulative manifestation only) were constructed as follows: it is assumed that the process lasts for total of time units. The instance manifestation probabilities were obtained from the Weibull p.d.f.: , with parameters and , following Eq. (2) in [36] (updated in [12]), normalized to 65 days: . The daily recovering probabilities were obtained from log-normal distribution , , normalized to 60 days, and then time-shifted for 4 days, obtaining: , for and , for . Parameters were chosen to match mean recovery time of days.

In the analysis, the continuous-time non-Markovian SEIS model was simulated using the integral form (5) (using Riemann sum (4)), while varying the parameter in the range. The Markovian SIS model (2) was numerically investigated, using forward Euler method, with parameter arbitrarily chosen, and parameter that varied in the region, as calculated from relation (27). Both model and stochastic simulations were conducted identically as in the discrete-time case, with the adequate adaptations. Integration constant was used in the analyzes.

Results of the parametric dependence analysis of the continuous-time case is presented in Fig. 3, accompanied by the results for the time series analysis in Fig. 4. In addition, the epidemic threshold analysis is presented in Table 3.

Fig. 3

Comparison between the steady-state values of the number of Exposed/Infectious/Infected individuals in the continuous-time non-Markovian SEIS model and the Markovian SIS model (parametric dependence analysis). Meaning of the symbols in the legend: nM — non-Markovian, M — Markovian; E — Exposed (number of), I — Infected (SIS)/Infectious(SEIS); calculated number of Infectious individuals in the SEIS process, from number of Infected individuals in the SIS process, assuming process equivalence; SS — stochastic simulations; the critical value of that defines the Epidemic threshold (refer to Table 3 for details): (a) BA1000 graph; (b) WS1000 graph.

Fig. 4

Comparison between the typical time series obtained for the number of Exposed/Infectious/Infected individuals in the continuous-time non-Markovian SEIS model and the Markovian SIS model. Meaning of the symbols in the legend: nM — non-Markovian, M — Markovian; E — Exposed (number of), I — Infected (SIS)/Infectious(SEIS); calculated number of Infectious individuals in the SEIS process, from number of Infected individuals in the SIS process, assuming process equivalence; SS — stochastic simulations: (a) BA1000 graph; (b) WS1000 graph. In all cases , ().

Table 3

Analysis of the epidemic threshold for the continuous-time cases and two different graphs: — critical value of the parameter , calculated directly from Eq. (15); — critical value of the parameter , calculated from the value of and Eq. (27), with .

DTPF set/Graph	β^NM,c	β^M,c	γ^M/β^M,c Ratio
BA1000	0.011508	0.0018895	5.2924
WS1000	0.018625	0.0030581	3.27000

As in the discrete-time case, the results from the analysis indicate that there is a perfect overlap of stationary state solutions, as well as the epidemic thresholds, of both the continuous-time non-Markovian SEIS model and the Markovian SIS model, providing relation (27) holds. However, one should exercise caution when considering the arbitrary choice of one of the parameters or , in the continuous time case: the application of relation (27), should result in or , providing or are arbitrary chosen, respectively. The minor model differences, in respect to the reality (stochastic simulations), indicate that model equivalence virtually implies a process equivalence, as well.

Epidemic threshold for the continuous-time non-Markovian SEIS model — extended analysis

In order to compare the theoretical result for the epidemic threshold of the continuous-time non-Markovian SEIS model, denoted with Eq. (14) and the reality, we have conducted a number of extended Monte Carlo simulation for different values of the parameter in the vicinity of the theoretical value. For each analyzed value of , one-hundred stochastic simulations (trials) were performed, with 500 (50%) of the nodes initially set as Exposed at . Two parameters were followed during the simulations: the number of non-zero outcomes at the end of each trial (at ), and the average number of Exposed nodes from the 100 trials, at each , . As critical, two values of were recorded: value at which at least one of the trials completed with non-zero outcome and the value at which at least 50% of the trials completed with non-zero outcome.

For the BA1000 graph (theoretical threshold calculated at ), for two trials completed with non-zero (2 and 14) number of Exposed nodes, with the average number of Exposed nodes, in 100 trials equaling 3.99 at and 0.16 at . For , 60 out of one-hundred trials completed with non-zero number of Exposed node, with the average number of Exposed nodes, in 100 trials equaling 30.78 at and 18.81 at .

For the WS1000 graph (theoretical threshold calculated at ), for two trials completed with non-zero (5 and 9) number of Exposed nodes, with average number of Exposed nodes, in 100 trials, equaling 5.42 at and 0.09 at . For , 51 out of one-hundred trials completed with non-zero number of Exposed nodes, with the average number of Exposed nodes, in 100 trials, equaling 29.45 at and 9.68 at .

Error analysis

Assumption of statistical independence of joint events, widely implemented in mathematical analysis of spreading phenomena on complex networks, is known to impose certain degree of error into analyzed models. Though efforts to assess the degree of error it generates, have been made in analysis of some models (for example see Appendix B in [18]), there is no definite answer to the size of this error depending on the state of the system. Authors in [25] concluded, that for the SIS process, the model obtained by utilizing this approximation upper-bounds the reality. In our previous experience, and supported by the analysis below, the size of this error is more noticeable in regions of low Infectiousness, while the approximated models provide a good assessment of reality in mid- to high- Infectiousness regions. One should note though, that in spite of its shortcomings, assumption of statistical independence of joint events is still considered to represent a golden standard in modeling spreading processes.

In order to support the notion that the assumption of statistical independence of joint events introduces larger amount of error in regime of low ineffectively, we conducted the following analysis: the obtained results for the cumulative manifestation cases were rearranged in order to quantify the relative model error as a function of the number of Infected nodes in the SIS case, and the number of Infectious nodes in the SEIS case, obtained from the stochastic simulations. The role of the Infectious nodes in respect to generally Exposed nodes, in the SEIS model, is emphasized in the error analysis, since Infectious nodes are directly involved in the process of agent spreading, i.e. the part of the process that is approximated and generates the error in the modeling.

Let denote the fraction of Infected/Infectious nodes obtained from the stochastic simulation of the SIS/SEIS process, and the accompanying model results. We define the relative model error as: The results of the analysis are presented in Fig. 5.

Fig. 5

Analysis of the model error as a function of the fraction of Infected (SIS) and Infectious (SEIS) individuals. Data obtained from CM1 and CM2 cases simulations.

Results of the analysis indicate that, given identical fraction of Infected (SIS) and Infectious (SEIS) nodes, the mean-field SEIS model introduces smaller degree of error, compared to the SIS model. However, since the endemic equivalence, as defined here, understands equality between the E state (SEIS) and I state (SIS), the SEIS model error tends to become evident when there is huge discrepancy between number of Exposed and number of Infectious nodes in the SEIS process. In the case of CM1, this ratio is 4.586. For the other DTPFs/CTPFs, this ratio is much lower (1.75 for CM2, 2.0 for RM, and 1.52 in the continuous analysis). Notice that this ratio is DTPF/CTPF specific, and is fully defined with the relation (A.4). From Eq. (A.4), one obtains that, for SEIS process in which level of overlap between recovery and manifestation time distribution is high (larger amount of nodes tend to recover, prior to becoming Infectious), ratio between number of E and I nodes increases, and the error imposed on the modeling is higher. Contrary, as the degree of overlap reduces, i.e. nodes tend to become Infectious, prior to recovery, the model error decreases, adequately, for these cases, one can claim that the stationary-state model equivalence, implies process equivalence as well. This is well supported by the results for the CM2,RM DTPFs (Figs. 1(b), 1(d), 1(e), 1(f), 2(b), 2(d), 2(e), 2(f)), and the utilized CTPFs in the continuous case (Fig. 3, Fig. 4).

Representation of Markovian SIS models as non-Markovian SEIS models

In this Section, as a result of secondary importance, We show that every classical (Markov) SIS model occurring on complex networks, may be represented as non-Markovian SEIS model, with proper selection of DTPFs/CTPFs and . Illustrative numerical example of this feature for the discrete-time case, has already been presented in [15]. In this article, this feature is theoretically investigated and shown for both model forms in a concise mathematical procedure.

As a note to the readers, in what follows, we avoid the use of label superscripts, since the whole procedure is conducted on the non-Markovian SEIS model, that under the investigated circumstances reduces to the Markovian SIS form.

Discrete-time case

Consider, as suggested in [15], , for all , and , with , with being the Heaviside function. Consequently, , , , . Acting similar as in [16], we obtain:

When , the last term in the equation vanishes. Bearing in mind that , the following relation holds: The last equation is the equation of the Markovian SIS model (1). To summarize, for arbitrary values of parameters and , the discrete-time SIS model, may be represented as non-Markovian SEIS model, providing and satisfy the relations defined in the introduction of this Subsection.

Continuous-time case

By considering , consequently , and and (see [16] for details) one obtains: The last relation is identical with the classical formulation of the Markovian SIS model occurring on complex network, in continuous time (2). The relation confirms that, for a given Markovian SIS, defined by parameters and , exists a non-Markovian SEIS model defined with the identical parameter and CTPFs and , such that the non-Markovian SEIS mimics the behavior of the Markovian SIS model on complex networks.

Typical time series generated by both models and the co-responding stochastic processes (reality) are presented in Fig. 6. The figure illustrates a perfect match between the models, with minor discrepancies in respect to the stochastic simulations (reality).

Fig. 6

Comparison between the typical time series obtained for the number of Exposed/Infectious/Infected individuals in the continuous-time non-Markovian SEIS model with and and the Markovian SIS model, for and . Meaning of the symbols in the legend: nM — non-Markovian, M — Markovian; E — Exposed (number of), I — Infected (SIS)/Infectious(SEIS); SS — stochastic simulations: (a) BA1000 graph; (b) WS1000 graph.

Conclusions

The non-Markovian systems more accurately address the spreading processes in comparison with Markovian models. This characteristic of non-Markovian models originates in their basic definition — to consider the status transitions that accompany the spreading as non-Poisonous processes, as confirmed by every-day practices. On the adverse side, the numerical analysis of non-Markovian processes is computationally more demanding. The computational complexity increases with the process memory, i.e. with parameter . Even further, in continuous models the choice of an accurate integration step may require substantial amount of computer memory, in order to store and manipulate with an excessive number of preceding states. When these features are accompanied by a huge network, the analysis of a non-Markovian, in this case SEIS, model, may become an overwhelming task for all, but a fairly small number of computing devices.

In this article we have shown that for the basic non-Markovian re-occurring model, the SEIS model, the stationary state distributions of nodes being Exposed/Infectious, may be found from a Markovian SIS equivalent. This result is of at-most importance, since it allows the numerical analysis of the stationary state solutions of the spreading processes to be conducted on systems with standard computational capabilities. In that sense, the result contributes to bringing the computational analysis of non-Markovian models closer to more users, especially individual researchers, small research teams and organizations, that may not be able to acquire a sufficiently powerful computing equipment.

CRediT authorship contribution statement

Igor Tomovski: Conceptualization, Formal analysis, Methodology, Software, Writing – original draft, Writing – review & editing. Lasko Basnarkov: Formal analysis, Methodology, Supervision. Alajdin Abazi: Supervision, Project administration.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.