Building epidemic models for living populations and computer networks

Abstract

Accurate modeling of viral outbreaks in living populations and computer networks is a prominent research field. Many researchers are in search for simple and realistic models to manage preventive resources and implement effective measures against hazardous circumstances. The ongoing Covid-19 pandemic has revealed the fact about deficiencies in health resource planning of some countries having relatively high case count and death toll. A unique epidemic model incorporating stochastic processes and queuing theory is presented, which was evaluated by computer simulation using pre-processed data obtained from an urban clinic providing family health services. Covid-19 data from a local corona-center was used as the initial model parameters (e.g. $R_0$ , infection rate, local population size, number of contacts with infected individuals, and recovery rate). A long–run trend analysis for 1 year was simulated. The results fit well to the current case data of the sample corona center. Effective preventive and reactive resource planning basically depends on accurately designed models, tools, and techniques needed for the prediction of feature threats, risks, and mitigation costs. In order to sufficiently analyze the transmission and recovery dynamics of epidemics it is important to choose concise mathematical models. Hence, a unique stochastic modeling approach tied to queueing theory and computer simulation has been chosen. The methods used here can also serve as a guidance for accurate modeling and classification of stages (or compartments) of epidemics in general.

Introduction

Main objective of this work is to provide a simplified guide for building models that cover all the stages of an epidemic outbreak. Due to the discrete nature of its elements, epidemiology contains several randomness. This fact encourages us to model events of epidemiology as stochastic processes in order to increase the accuracy of the model at hand. The model presented here is not only eligible for producing reliable results, it also provides analysis of future tendencies. Another unique contribution of this paper is the formulation of the state transitions of the stages of an epidemic: (a) from healthy to susceptible, (b) from susceptible to infected, (c) from infected to quarantined, (d) from quarantined to recovered, and (e) from recovered to healthy. Furthermore, efficiency analysis of the recovery process dealing with different service/qualities is also discussed.

We do not specifically discuss the characteristics of some known approaches, for example, deterministic and stochastic approaches. We intended to design a more generic model, which can be applied to building models for human epidemics and computer epidemics with minor modifications. The approach presented here is based on stochastic processes and queueing theory. Stochasticity possessed inherently by epidemiology, and flexibility of applying stochastic models to a wide area of sciences encourages us to choose this approach for modeling epidemics on human and computer networks.

We computer engineers and medical doctor have limited knowledge of stochastic processes and queueing theory, which provide valuable methods and techniques for modeling stochastic systems. We usually describe event dynamics using models that are formalized mathematically. Although, inevitable, rigid mathematical formalization can sometimes be cumbersome and/or even disruptive when applying to experimentation and system implementation. There are, naturally, cases where complex theoretical proofs are necessary. However, there also exist areas of science that have already been discussed widely using a multitude of theories and their proofs. To obtain maximum benefit from the material presented here, we intend to simplify mathematical presentations and avoid theoretical proofs so that the models and approaches given can be easily adapted to various purposes, for example, experimentation, simulation, and resource management.

As detailed in Kondakci and Dincer, ¹ lifecycle of an epidemic system generally consists of five stages, healthy, susceptible, infected, quarantined, and recovered. These are modeled as distinct states taking place in a chain of interrelated stochastic processes. The states are in most cases recurrent, omitting the case of death (an absorbing state). In any case, an infection case goes dynamically through certain transitive states, that is, moving from one state to another. Each state is numerically described either by a probabilistic function or by an integer count. The probabilistic function represents the probability of estimating the value of a given parameter, for example, probability of being infected among 100 susceptible individuals. Or, it can also describe the probability of being in a state, for example, probability of being recovered within a specific time duration. A state can also be represented by the number of increments/decrements in sizes of a parameter, for example, infected population size, the state of infectious patients undergoing quarantine operations, and the number of recovered patients.

To estimate related parameters, we generally build mathematical models fitted to some data set. Unfortunately, this is more difficult than it sounds. Because errors in statistical fitting and uncertainty in the parameter estimations are inevitable. ² Epidemiology is often limited to the study of disease spreading in populations of humans or animals. Several commonalities between the human and computer epidemiology have recently emerged. In order to gain a general perspective, we consider that computer worms/viruses (malware) and epidemic threats to living populations (human, animal) are interchangeable. Therefore, assuming readers are from multidisciplinary environments, we avoid complex mathematical manipulations and in-depth proofs of theories, and instead focus on the applicability and correctness of the presented models. The aim is thus to provide a set of general formulae and methods that can be easily adapted to special cases in order to perform realistic trend analysis in epidemics. Simulations and experimental analysis of various cases will then be facilitated by simply modifying the general model provided.

We use stochastic approaches in conjunction with queueing theory to analyze the stages of epidemic outbreaks. Models discussed are generally substantiated by simulations and numerical examples, some of which are given here. If needed, one can refer to Taha ³ and Medhi ⁴ for a comprehensive coverage of Operations Research and Stochastic Models in Queueing Theory, and Cui et al. ⁵ for Markov models used in a multi-state repairable systems.

Motivations and related work

In–depth study of epidemiology is an extremely complex area of research. Contemporary infectious diseases evolve sporadically across a wide spectrum, for example, Covid-19 pandemics,^6–8 and fundamental concepts requiring detailed study are infection spread, surveillance, treatment, control, immunity in populations, and resource management. With the Covid-19 pandemics, the world has realized the importance of the efficient resource management for health institutions, where most countries have shown significant deficiencies in their health systems. For instance, even though it has relatively efficient health system, Turkey has locked down all the polyclinics of the major hospitals throughout the country and mobilized the resources to thwart the pandemics.

So called Recurrent Epidemic Model (REM) first introduced, by Kondakci and Dincer, ¹ applies a comprehensive stochastic analysis to computer epidemics, processing through a detailed analysis of a five stage stochastic model, where each stage is mathematically explored to obtain a satisfactory formalism in all dominating stages in the epidemiology. REM presents an extended structure of state transitions in order to highlight the shortcomings of the classical and deterministic models. It aims to obtain higher accuracy by use of a comprehensive model containing a sequence of stochastically dependent compartments. Other interesting approaches have been discussed in Amador and Artalejo,^9,10 Yang et al., ¹¹ and Sellke et al. ¹² An epidemic model, ¹³ applying the finite-Markovian process gives an overview of the steady-state analysis of the evolution of an epidemic case ¹³ has proposed concrete solutions for determining transition dynamics by formalizing the parameters for the infection and recovery rates.

Additionally, ¹⁴ discusses the detection of worm spread, with suffeciently detailed deterministic and stochastic models. Several stochastic approaches have also been proposed to analyze dynamics of epidemics, for example, Ortega et al. ¹⁵ To gain an insight into the stochastic models, the reader can examie the survey on stochastic epidemic models presented in Britton. ¹⁶ Mathematics and simulation are effective tools assisting us in almost every field of research, in particular, for building models and substantiating the model results with regard to the theoretical background. A review of mathematical models dealing with malware propagation in computer networks is presented in del Rey. ¹⁷

Classical epidemic models (e.g. SIR, SIS, SEIR) maintain their relevance within the field of epidemic research, providing a simpler and more structural overview of epidemic outbreaks. Related to this, Artalejo and Lopez-Herrero ¹⁸ discusses the dynamics of susceptible-infective-removed (SIR) model with the aid of continuous time Markov chains. A cellular automata based SEIR model for computer virus spreading is proposed in Batista et al. ¹⁹ A stochastic approach of analysis covering SIR models is considered in Kenah and Robins. ²⁰

Several important characteristics are common to almost all fields of research, such as mathematical biology, computer science, physics, economics, and the socail science. ²¹ Therefore, it is hard to confine epidemiology research to a particular field, without restricting the methods used to analyze cases that can lead to inaccurate results. For example, Mukhopadhyay and Bhattacharyya ²² discusses a comprehensive mathematical model considering the propagation of infections in human populations, which applies the theory of (traveling) waves to build the epidemic model. A comprehensive study of the mathematical epidemiology containing several models is presented in Hethcote. ²³ Other interesting models can be found in Paul and Mishra, ²⁴ Ren et al., ²⁵ Kudo et al., ²⁶ and Wanget al. ²⁷

A computer worm is an intelligent virus that can automatically spread through vulnerable (or unprotected) computer applications. Some protection systems are incapable of detecting certain viruses, which are classified as polymorphic viruses, such as Stuxnet. ²⁸ Worms are often deployed en–mass by mass–mailing in order to cause system unavailability. ²⁹ A survey on taxonomy and characterization of computer worms is presented in Fachkha and Debbabi. ³⁰ Numerous computer worm models exist in the literature, most of which are rooted in the classical biological models.^31–34

It is also important to study the efficiency of quarantine strategies under various infection profiles. To do so, we need to have an additional stage (compartment) dealing with modeling of the quarantine process. This paper deals with a stochastic model with five compartments progressing from the healthy state to the recovered state in order to cover all the major stages in general. Effect of the incubation time on the transmission of the p. vivax malaria is discussed in Nah et al. ³⁵

As also discussed in Kondakci, ¹³ many solutions assume homogeneous systems with constant transition rates, where states change with invariant probabilities. However, it has been frequently observed that both the threatened objects and their responses are heterogeneous, where states change with varying probabilities.

Communication patterns in a network is important for the propagation of viral infections. For example, the total number of bites made by mosquitoes has an important impact on the incubation time and transmission probability of malaria. ³⁶ Immunization is naturally an important part during the course of an epidemic spreading. Accordingly, Pastor-Satorras and Vespignani ³⁷ considers epidemics and immunization in scale–free networks applying a rigid mathematical analysis. Scale–free networks related to epidemic dynamics is considered in Fu et al. ³⁸ Network traffic characteristics for some computer viruses may not be as important as we assume. Regarding this, Sharif et al. ³⁹ discusses some approaches used in computer epidemics to generate propagation methods that result in successful worm outbreaks, ignoring the traffic characteristics of networks.

Models that do not incorporate the incubation period, and assume only deterministic distributions of the compartmental parameters and state changes will generally fail when fitting data to the initial outbreak data.

Susceptible nodes do not have latent periods. Since the infection cases behave differently regarding the population and infection types, the related model must accommodate the evolution of each specific case. The model presented in Kondakci and Dincer ¹ discusses methods to overcome the limitations by introducing additional parameters to model the incubation periods.

For example, incubation period of computer viruses generally depends on the deployment of the individual virus and the immunity (or protection level) of the victim system. Assuming a malcode (virus, worm, trojan) embedded in an e-mail attachment, and the victims (users) of the receiving host delay opening the malcoded attachments, then the latent time might eventually be very long. Some aggressive worms like the well-known CodeRedII ⁴⁰ worm, W32.Goner.A@mm, and SQL–Slammer ⁴¹ may immediately infect systems (zero latent time) without sufficient protection.

Population structures

Epidemic analysis deals with interrelated classes of populations from different networks, where the elements of each network are subject to different treatments in different stages during the course of an epidemic case. Two different types of networks are considered here, centralized and scale–free networks. Nodes in a scale–free network are connected through short paths, where the number of connections emanating from a given object obeys the power-law distribution. In epidemiology, “short path” means that the connectivity weight of two neighboring nodes is reasonably short, or the frequency of contact is sufficient to cause an infection, and therefore proper control of the epidemic spread requires more effort and resources. In contrast to scale–free networks, centralized networks have low degree of connectivity among the nodes, that is, fewer number of individuals in a population. Hence, regarding the epidemiology, nodes in a centralized networks have low rate of infection compared to the scale–free networks.

Centralized networks

Centralized networks are those that encompass many elements with low degree of connectivity. Typically, human epidemiology from small-sized communities exhibits the characteristics of centralized networks, because connectivity factors among the individuals are small and the epidemic spread is more under the control of individuals. Humans’ inherent intelligence under many circumstances can limit their access to the pathogenic environment. There are also exceptional populations who are unable to contact, for example, military wards, jails, mass production plants, complex residential areas, farms, and animal flocks. Thus, a parameter denoting the degree of contact for such networks has a significant role in the epidemiology of dense networks (e.g. intensely crowded populations).

Scale-free networks

Scale-free networks (SFN) are those that follow the properties of the power low distribution, where there is an exponential relation between the degree of connectivity of a node and the frequency of connections it makes. More on discussions about power low can be found in Newman. ⁴² Scale-free networks are regularly present in all forms of networks such as computer networks, social networks, and biological systems. SFNs have the property of power law, where a quantity changes as a power of another, that is, $P(k)=k^{-\gamma }$ , which denotes a tight relationship between two variables, where a change in one causes a change in the other. The changes are mutual, relatively proportional, and independent of the initial conditions. For example, the probability that an infectious node contacts $k$ susceptible hosts obeys the power-law distribution as $P(k)=k^{-2.06}$ . Unfortunately, analysis and characterization of power-law data is very complex; least-square fitting method is often used to analyze the power low data, which usually leads to uncertainty and inaccuracy in the results. Therefore, it is important to apply a through statistical analyses on the empirical data in order to recognize its actual characteristics, and whether it accurately fits the power-law distribution. A statistical framework applying the method of maximum-likelihood fitting with the Kolmogorov–Smirnov statistics is used in Clauset et al. ⁴³ for discriminating power-law patterns in some empirical data.

The power low property is also associated with preferential attachment, which refers to the distribution of properties among the nodes in a network. Scale-free networks generally grow progressively by randomly joined nodes, where each node links to current hosts with preferential attachment. That is, the probability of connecting to host $i$ is proportional to the number of currently connected links $k_i$ that host $i$ connected to Barabási, ⁴⁴ that is,

$$P(\mathrm{joining}\mathrm{node}i)=\frac{k_i}{{\sum }_j{k}_j}.$$

For example, distribution of resources in a society will often centralize at a node with a higher degree of connections. A typical example of this is business people’s preference for commercials on highly attractive environments (e.g. TV channels, newspapers, or Web sites). Another example is the prevalence of malware attacks on highly connected computer networks in order to gain access to greater numbers of computers. Lower immunity facilitates the “invitation” of infectious diseases.

The rate of an infection a node undergoes depends on the number of contacts with infected nodes that it may have, as well as with the duration (or frequency) of the contacts. Assuming that each contact has the same infection rate denoted by $\beta$ , then the rate at which a node from the susceptible group $i$ gets infected is given by

$$\lambda =\beta \sum _j^N{C}_{i,j}{I}_j,$$

(1)

where $I_j$ is the proportion of infectious nodes in population $j$ , $C_{i,j}$ denotes the number of contacts per time slot each node in group $i$ contacts nodes in population $j$ having the size of $N$ .

Growth and extinction dynamics

Let us consider the classical SIR (Kermack–McKendrick model ³² ) model and expand it to contain the dynamics of healthy $h(t)$ and quarantined $q(t)$ compartments. With this expanded SIR model (ESIR), we can explore significant details about the long-run dynamics of the state transitions of the system. Given the total number of nodes, $N(t)$ , at time $t$ , the number of healthy nodes is determined by summing the respective fractions from all compartments as

$$h(t)=N(t)-[s(t)+i(t)+q(t)+r(t)],$$

(2)

where $s(t)$ denotes the number of susceptible, $i(t)$ infected, $q(t)$ quarantined, and $r(t)$ recovered (removed) nodes. The fallowing differential equation set, equation (3), represents the dynamic behavior of the ESIR model.

$$\begin{array}{c}\mathrm{Healthy}:\frac{d}{\mathrm{dt}}h(t)=\phi N(t)-\delta h(t)+\mu r(t)+\epsilon s(t)\hfill \\ \mathrm{Susceptible}:\frac{d}{\mathrm{dt}}s(t)=\delta h(t)+{ϵ}_{s}r(t)-\lambda s(t)-\epsilon s(t)\hfill \\ \mathrm{Infected}:\frac{d}{\mathrm{dt}}i(t)=\lambda s(t)-\alpha i(t),\hfill \\ \mathrm{Quarantined}:\frac{d}{\mathrm{dt}}q(t)=\alpha i(t)-\gamma q(t),\hfill \\ \mathrm{Recovered}:\frac{d}{\mathrm{dt}}r(t)=\gamma q(t)-\mu r(t).\hfill \end{array}$$

(3)

The parameters are self explanatory, for example, population change in a compartment is interrelated with the parameters of its predecessor and ancestor compartments. For example at time $t$ , the growth rate of the infected population i(t) depends on the number of the susceptible objects by a factor of $\lambda$ and the current size of the susceptible population $s(t)$ , minus the number of nodes removed from the infected population to the quarantine population $q(t)$ having the rate $\gamma$ , that is, $\gamma q(t)$ .

A distinct property of the ESIR model is that a very small proportion, ${ϵ}_{s}r(t)$ , of the repaired hosts that will be moved to the healthy population may again become susceptible. This can be considered as a loose user vigilance or degraded immunity due to lack of vaccination.

Basic reproduction number

The concept of basic regeneration number was first used in conjunction with demography in 1880s to predict the growth rate in specific sections of populations. ⁴⁵ Regarding the epidemiology, it expresses the number of potential transmissions from a single infected node during its transmitter period, and is denoted by $R_0$ , meaning an offspring that starts at its first birth ( $0$ ) and reproduces itself with rate $R_0$ in time. There is no fixed formula for $R_0$ , which varies considerably from case to case, and even from population to population, taking on different parameters. Therefore, it is generally difficult to match its theoretical form to its practical use. In its simplest form, it is defined as a function of three parameters as

$$R_0=〈\beta ,\lambda ,\mu 〉.$$

Where, $\beta$ denotes the effective contact rate (i.e. efficacy of infection), $\lambda$ infection, and $\mu$ recovery rate. We can then search for the changes in the parameters that lead endemic equilibrium, where the infection is constantly maintained with reproduction period of the disease.

There exist numerous methods to define $R_0$ , however which of these is the most accurate and applicable to a given case is frequently disputed by researchers. However, there is no doubt that $R_0$ has a remarkable role during the initial phase of epidemic spreading. As a quick reference, one can refer to Driessche and Watmough ⁴⁶ for an elaborated work on the basic reproduction number, used in compartmental infection models. A preliminary estimation of $R_0$ of the novel Coronavirus (sars-cov-2) is presented in Zhao et al. ⁴⁷ Considering varies categories of hosts, it is claimed by Diekmann et al. ⁴⁸ that the next-generation matrix is the basis for defining and determining $R_0$ , in which authors follow an instructive process to derive $R_0$ . A useful approach to dealing with the definition of $R_0$ for a variety of infectious diseases is also considered in Driessche. ⁴⁹ As another initiative, ⁵⁰ determines $R_0$ from the data obtained during the initial phase of the infection outbreak, where the discussion is particularly focused on computing $R_0$ for vector-borne diseases. Given the part of the SEIR model below, also detailed in Heffernan et al., ⁵¹ we can easily determine $R_0$ using the next-generation method.

$$\begin{array}{c}\frac{d}{\mathrm{dt}}e(t)=\beta s(t)i(t)-(\mu +k)e(t)\hfill \\ \frac{d}{\mathrm{dt}}i(t)=\mathrm{ke}(t)-(\gamma +\mu )i(t),\hfill \\ R_0=\frac{k\beta \zeta }{(\mu +k)(\mu +\gamma )\mu }.\hfill \end{array}$$

(4)

Where, $\mu =$  death rate, $\beta =$  efficacy of infection, $k=$  infection rate, $\gamma =$  recovery rate, and $\zeta =$  birth rate of susceptible hosts. ⁵¹

According to World Health Organization (WHO), Vector-borne infections make up 17% of all infectious diseases; these are transmitted by mosquitoes, triatomine bugs, blackflies, tsetse flies, sandflies, mites, ticks, lice, and snails. Hence, finding an appropriate formula of $R_0$ regarding these conglomerate populations would be NP–hard, if possible. Extremely complex mathematical analysis are currently popular, but may be cumbersome to apply in general. For example, Yang ⁵² uses jacobian and next generation matrices to define $R_0$ for the propagation of Dengue (a flavivirus) virus transmitted by the interaction between mosquito and human populations in urban areas. Defining $R_0$ for complex diseases in wildlife is presented in Hartemink et al., ⁵³ which focuses on the complexity of Tick-borne infections that have multiple transmission routes. The pros and cons of approaches used for estimating $R_0$ based on the initial rate of spread of influenza A (H1N1) 2009 in Japan are discussed in Nishiura et al. ⁵⁴ A model for the highly pathogenic H5N1 avian influenza virus is discussed in Colizza et al., ⁵⁵ where each compartment of the influenza case is described by a set of stochastic equations. The basic reproduction number was also specifically defined for the case occurring at the beginning of 2006 as

$$R_0=\frac{\beta }{\mu }(r_{\beta }{P}_a+1-P_a).$$

(5)

Where, $\beta$ denotes the transmission parameter, which depends on both time and geographical zone. Increase and decrease in population probabilities are denoted by $P_a$ and $1-P_a$ , respectively. The equation is based on the method given in Diekmann and Heesterbeek. ³⁴

Epidemic threshold

This is conditional on the fact that prior to reaching an epidemic outbreak, there should be a large number of susceptible hosts exposed to infection, regardless of the probability of an initial outbreak. An epidemic threshold can be reached when at least one infectious host is in contact with susceptible hosts for a specific duration. Branching process is a useful approach for modeling many physical events, as well as epidemic outbreaks. An initial epidemic outbreak often evolves as a branching process in a crowded population without preventive control. For example, infection has a slower changing stochastic behavior in a human population compared to poultry. On the other hand, unprotected computer networks have a rapid evolving of infection characterized as a typical branching process. In most cases, effective infectious periods for computer networks are brief, whereas these periods may last a lifetime for human populations. Predicting the epidemic threshold and the basic reproductive number in such cases can be extremely complicated, if a certain threshold of observation data is not available.

The probability of a single infection producing a new infected node with rate $\lambda$ during $\mathrm{\Delta }t$ is defined as, ⁵⁶

$$P_1(\mathrm{\Delta }t)=\lambda \mathrm{\Delta }t+o(\mathrm{\Delta }t),$$

whereas the probability of a node without infection is expressed by

$$P_0(\mathrm{\Delta }t)=1-\lambda \mathrm{\Delta }t+o(\mathrm{\Delta }t),$$

$o(\mathrm{\Delta }t)$ denotes an infinitesimal higher order than $\mathrm{\Delta }t$ , saying that $o(\mathrm{\Delta }t)$ has the property $\underset{\mathrm{\Delta }t\to 0}{lim}\frac{o(\mathrm{\Delta }t)}{\mathrm{\Delta }t}\to 0$ , that is, it converges to zero much faster than $\mathrm{\Delta }t$ . We can now express the infection growth as an exponential function by setting ${\lambda }_1=-\lambda$ and ${\lambda }_2=\lambda$ so that $\sum _j{\lambda }_j=0$ , hence, the probability $P_{\mathrm{kn}}(t)$ producing $n$ infections at time $t$ is obtained by

$$P_{\mathrm{kn}}(t)=e^{-\lambda t}(1-e^{-\lambda t}{)}^{n-1},n\ge 1.$$

(6)

Now, let

$$P(t)=\mu \mathrm{\Delta }t+o(\mathrm{\Delta }t)$$

be the probability of recovering from each disease with rate $\mu$ in a small time slot $\mathrm{\Delta }t$ . We can also consider the recovery process as a diagnose process with negative outcomes (i.e. no infection). Based on the balance property of a Markov process we have

$${\lambda }_0=\mu ,{\lambda }_1=-(\lambda +\mu ),\mathrm{and}{\lambda }_2=\lambda .$$

Thus, the probabilities for zero and $n$ infections occurring after the time slot $t$ are obtained as follows:

$$\begin{array}{c}\xi =\lambda -\mu ,\hfill \\ P_0(t)=\mu \gamma ,\hfill \\ P_n(t)=(1-\lambda \gamma )(1-\mu \gamma )(\lambda \gamma {)}^{n-1},\forall n\ge 0,\hfill \end{array}$$

(7)

where

$$\gamma =\{\begin{array}{c}\frac{1-e^{\xi t}}{\mu -\lambda e^{\xi t}},\mathrm{if}\lambda \ne \mu ,\hfill \\ \multicolumn{1}{c}{\frac{t}{1+\lambda t},\mathrm{if}\lambda =\mu .}\end{array}$$

The probability of zero production, $P_0$ , is conditional on the relation between the infection and recovery rates.

$$P_0=\{\begin{array}{c}1,\mathrm{if}\mu >\lambda ,\hfill \\ \multicolumn{1}{c}{\frac{\lambda }{\mu },\mathrm{if}\mu <\lambda .}\end{array}$$

Rather than the notion of the basic reproductive number $R_0$ , an effective alternative indicator for an epidemic outbreak can be considered as epidemic threshold. Fitting $R_0$ to general cases cannot be as accurate as that of the point where an epidemic threshold is identified because the threshold method can be easily estimated by observation or collecting an appropriate amount of data of a given epidemic case. Figure 1 illustrates the threshold observations during eight time slots (e.g. days) of a pandemic for different population structures. As can be noted from the shaded (threshold) region, each of the populations shows a different growth rate within the first two time epochs. For example, the first curve (population (a)) has a sharper slope, Covid-19 initial spread rate in Turkey, which immediately jumps beyond the threshold level during the first two time slots. On the other hand, population (c) has a slowly increasing growth rate, where the threshold level is reached at the end of the second time slot. If we have sufficient data about different infection cases, we can easily determine the corresponding epidemic threshold using the method of differentiation from first principles (Delta method). It should be noted that threshold regions are specific to epidemic cases and often determined empirically by use of earlier data set, if available. That is, similar to $R_0$ , every case has its own epidemic threshold and one should avoid the generalization before the data set has been strictly validated whether it fits the case under consideration.

Figure 1.

Infection growth rates for different epidemic threshold values: (a) Covid-19 initial spread rate in Turkey, (b and c) different forms of typical influenza.

Stages in epidemic cases

Spreading and recovery operations i a typical epidemic case follow a set of contiguous and interdependent stages, where each stage is modeled as a dynamic subsystem with adequately chosen state parameters. Each subsystem has some input and output parameters which change dynamically throughout the lifecycle of the epidemic case. In most cases, the subsystems are transient, and each subsystem feeds the input of its ancestor subsystem with data, while receiving data from its predecessor subsystem. Data exchanges between the subsystems take place through the individual state parameters of the associated subsystem. Each subsystem has its own population of elements, where the number of elements in a subsystem is denoted by a numerical value called “state.” For example, state $n$ means that the related subsystem contains $n$ elements, and state transition denotes the change in the number of elements in a subsystem. Table 1 contains the state space for all populations of the generic model. A typical state transitions diagram shown in Figure 2, illustrates the interaction between the model subsystems. Following each state transition, the individual population sizes and some parameter values are subject to change. These parameters denote state transition probabilities affected by factors such as immunity, contact rate, and quarantine and recovery some efficiency of the hospitalization unit. Figure 2 shows the chain of state transitions. An associated probability distribution function is denoted together with the respective block (compartment). Each of these distributions/equations describes the state transition dynamics of the related compartment.

Table 1.

Compartmental states associated rate parameters.

Parameter space for state changes
Flow rate	Description	Population size changed
${\lambda }_{\mathrm{hs}}$	Healthy to susceptible	Susceptible: $s(t)$
${\lambda }_{\mathrm{sh}}$	Susceptible to healthy	Healthy: $h(t)$
${\lambda }_{\mathrm{si}}$	Susceptible to infected	Infected: $i(t)$
${\lambda }_{\mathrm{iq}}$	Infected to quarantine	Quarantined: $q(t)$
${\lambda }_{\mathrm{qr}}$	Quarantined to recovery	Recovered: $r(t)$
${\mu }_{\mathrm{rh}}$	Recovered to healthy	Healthy: $h(t)$

Figure 2.

A compartmental representation of the states of the recurrent model.

As can be noted, elements of the populations iterate from state to state depicting a recurrent characteristic, where the population sizes are precisely expressed by the following state equations, equation (8).

$$\begin{array}{c}h(0,T]=h(t)+\sum _{t=0}^T[{\lambda }_{\mathrm{sh}}\times s(t)+{\mu }_{\mathrm{rh}}\times r(t)-{\lambda }_{\mathrm{hs}}\times h(t)],\hfill \\ s(0,T]=s(t)+\sum _{t=0}^T[{\lambda }_{\mathrm{hs}}\times h(t)-{\lambda }_{\mathrm{sh}}\times s(t)-{\lambda }_{\mathrm{si}}\times s(t)],\hfill \\ i(0,T]=i(t)+\sum _{t=0}^T[{\lambda }_{\mathrm{si}}\times s(t)-{\lambda }_{\mathrm{iq}}\times i(t)],\hfill \\ q(0,T]=q(t)+\sum _{t=0}^T[{\lambda }_{\mathrm{iq}}\times i(t)-{\lambda }_{\mathrm{qr}}\times q(t)],\hfill \\ r(0,T]=r(t)+\sum _{t=0}^T[{\lambda }_{\mathrm{qr}}\times q(t)-{\mu }_{\mathrm{rh}}\times r(t)].\hfill \end{array}$$

(8)

As also illustrated in Figure 2, some fractions of elements are flowed from population to population at different intensities, for example, from healthy to susceptible at rate ${\mu }_{h_s}$ and from susceptible to infected at rate ${\lambda }_{\mathrm{si}}$ .

State variables contain the time-dependent number of objects in a given epoch $(t,T]$ for the related state. For example, as depicted by equation (8), during the time slot $(0,T]$ , a proportion of the recovered population is moved to the healthy population, while some healthy objects remain susceptible. Likewise, a proportion of the susceptible population will become infected, some of the infected objects are quarantined, and a proportion of the quarantined objects recover. Consequently, at time $t+T_x$ , the number of objects in the respective population will be changed in a recursive manner. ¹

Using notations from queuing theory, the stages (compartments) of the generic model will be described in the following sections. Arrivals of an object to a stage and operation completion (departure) from a stage are modeled as an M/M/1:(FCFS) queueing system. Based on this first-come first-served (FCFS) principle, each stage in the system encompasses a cascaded queueing system, where each stage serves a unique function. The overall behavior of the system is illustrated by the flow diagram shown in Figure 2, where each stage is a dynamic component of the recurrent operations. Initially, at least one host must be infected in the Susceptible population for further development of a disease. Often a stochastic growth of the susceptibility takes place as a branching process. Thereafter, an epidemic outbreak (exponentially distributed process) is initiated in the susceptible population. Arrivals to the infection stage and most other stages are modeled as a Poisson process denoted as ${\lambda }_{x_y}$ , meaning that the arrival is from stage $x$ to stage $y$ . The compartments of the recurrent model, shown in Figure 2, will be detailed in the following sub-sections.

A generic model

In order to facilitate the analysis in a broader scope, a generic model is necessary, which can elucidate the most encountered stages of epidemiology. In order to support this idea, the model presented here is built upon five types of populations: (i) healthy, (ii) susceptible, (iii) infected, (iv) quarantined, and (v) recovered. Furthermore, we have two major assumptions: (a) all individuals/nodes in a sample population are initially susceptible, (b) if a node gets infected then it becomes the source (transmitter) of a branching process for spreading the infection. At the next time slot, some of the hosts in the population are susceptible, and some are transmitter, which can cause infections. Initialization of the susceptible population with empirically obtained data is necessary in order to keep analysis within an acceptable dynamic range until the equilibrium state is achieved.

Unless otherwise specified, the process of infection through recovery is considered as an M/M/1 queueing system with a single server having Poisson arrival (infection case) and exponentially distributed recovery service time (recovery case). More specifically, the system is an FCFS Markovian queue with exponentially distributed interarrival and departure (or service) rates: arrivals and departures are independently and exponentially distributed with rates $\lambda$ $\mu$ , respectively. Hence, the service utilization factor is $\rho =$  “arrival rate/service rate,” that is, $\rho =\lambda /\mu$ .

Two distinct populations (biological and computer) may have different levels of susceptibility, where the susceptibility of an object depends on the protection level for computers and immunity level for living individuals. The immunity level for living individuals has several parameters such as age, locality (e.g. kindergarten, chicken poultry), mobility, vaccination, and so on.

Spread and extinction dynamics of a disease can be analyzed by two fundamental methods, time-dependent and limiting (or equilibrium) behavior analyses. Analysis of time-dependent behavior can help determine the number of patients (or hosts) under care by computing the probability that at time $t$ there are $n$ hosts in the system. System dynamics can be easily traced by the continuous Markov chain modeled as a birth–death process. A process is in state $n$ means that population (e.g. infected population) size is $n$ . Transition from state $n$ to state $n+1$ is denoted as a birth, whereas a transition from state $n$ to $n-1$ depicts a death. Both births and deaths occur only to the adjacent state (e.g. from $x$ to $x+1$ or from $x$ to $x-1$ ) with rates $\lambda$ and $\mu$ , respectively. Birth and death rates are independent of time and of each other, and only dependent on the current state. The following definitions are used to describe the basic functions for the state transitions taking place in a given phase.

Definitions

• $P_0(t)\equiv$  Probability of no state change during time $t$ ,

• $P_1(t)\equiv$  Probability of a single state change during time $t$ ,

• $\lambda \mathrm{\Delta }t\equiv$  Probability of an arrival (birth) with rate $\lambda$ in $\mathrm{\Delta }t$ ,

• $\mu \mathrm{\Delta }t\equiv$  Probability of a departure (death) with rate $\mu$ in $\mathrm{\Delta }t$ ,

• $1-\lambda \mathrm{\Delta }t\equiv$  Probability of no arrival in $\mathrm{\Delta }t$ ,

• $1-\mu \mathrm{\Delta }t\equiv$  Probability of no departure in $\mathrm{\Delta }t$ .

Healthy to susceptible

All starts with a susceptible object having contact with one or more transmitter objects. The susceptible objects can be infected and become an active transmitter after the latent period. Following this, an epidemic outbreak will occur. The outbreak is best modeled as a branching process, see Athreya and Ney ⁵⁷ and Kalinkin ⁵⁸ for a general review of branching processes. A branching process is a regenerative Markov process, where an individual in a population goes through a recursive regeneration process that proliferates in the next generation with a specific probability.

Let us assume a group of transmitter objects have contact with some susceptible objects, where the contacts and eventual infections are independent of each other. Here, the probability of a contact causing $k$ infections is denoted by $p_k(t)$ . Following the initial infection phase, $n-k$ nodes out of $n$ become susceptible, which will then initiate an epidemic outbreak as a branching process ( $\xi$ ). The branching process evolves as a pure birth (Poisson) process with the basic transition rates expressed as

$${\xi }_{\lambda }=\{\begin{array}{c}{p}_{i,i}=-\lambda ,\hfill \\ \multicolumn{1}{c}{p_{i,i+1}=\lambda ,}\\ \multicolumn{1}{c}{p_{i,j}=0,\mathrm{for}j\ne (i,i+1).}\end{array}$$

Besides these one-step transitions, transitions from $i$ to $j$ occurring in $n$ steps will approximate to a Bernoulli process performed $n$ times having the probability of success $\lambda$ given by

$$P\{\xi =j|n\}=\{\begin{array}{c}\left(\begin{array}{c}n\hfill \\ \multicolumn{1}{c}{j-1}\end{array}\right){\lambda }^{j-1}(-\lambda {)}^{n-(j-i)}\hfill \\ \multicolumn{1}{c}{0,\mathrm{otherwise},\mathrm{for}0\le j-1\le n}\end{array}$$

(9)

The Poisson approximation of the above distribution is

$$P_{\xi }(j)~\frac{{\lambda }^j}{j!}{e}^{-\lambda },\lambda =n{p}_s.$$

(10)

A small proportion of the susceptible population may be well-protected or may have sufficient level of immunity, and hence, will not become infected. Hence, the proportion of the susceptible hosts returning to the healthy population is depicted by $R_{s_h}$ , which depends on the current population sizes ( $R_{s_i}$ , $R_{\mathrm{ss}}$ ) and rate parameters, susceptible–to–infected ( ${\lambda }_{\mathrm{si}}$ ) and susceptible–to–susceptible ( ${\lambda }_{\mathrm{ss}}$ ) given by

$$\frac{\partial R_{s_h}}{\partial t}=R_{s_h}-R_{s_i}(t){\lambda }_{\mathrm{si}}-R_{\mathrm{ss}}(t){\lambda }_{\mathrm{ss}}.$$

(11)

Susceptible to infection

Assume that we have a population consisting of $N$ susceptible and $N-k$ homogeneously infected hosts each infecting some susceptible hosts independently of all the others. Each of the transmitters initially cause a small probability $p(t)$ of infection during time slice $t$ . Assuming a Poisson process with very large $N$ then the probability $p(t)$ of infection becomes very small (property of Poisson distribution contra Binomial), then the distribution of infection described by process $\xi (t)$ becomes accurately a Poisson distribution, especially when $N$ approaches to infinity. Thus, the probability of finding exactly $k$ infections will become

$$P[\xi (t)=k]=\frac{n^{k}t}{k!}{e}^{-\mathrm{nt}},k=0,1,2,\dots ,$$

(12)

However, if the population has a small number of individuals, in which each infection being considered as a “success” in a Bernoulli process. Then, $\xi (t)$ will hold the number of successful outputs in a series of $N$ Bernoulli trials performed in $t$ time units.

$$P[\xi (t)=k]=\left(\begin{array}{c}N\hfill \\ \multicolumn{1}{c}{k}\end{array}\right)p(t{)}^k(1-p{)}^{N-k},k=0,1,2,\dots ,$$

(13)

As can be noted, the infection phase consists of two sub–phases: first, some susceptible hosts, which are vulnerable (or unprotected), will be infected with a rate depending on the infection probability and density of the susceptible nodes. A second sub-phase consisting of a virus spread that starts through the currently infected nodes, as a branching process expressed by the following distribution, ¹

$$P_k(t)=e^{-\lambda t}(1-e^{-\lambda t}{)}^{(k-1)};k\ge 1,$$

(14)

where $\lambda$ denotes the branching rate and $k$ is the current number of offspring generated from a single infection at time $t$ . Equation (14) is based on a single infectious node that initiates the branching of new infections. It should be reminded that, each infectious node causing a new infection in a small epoch $\mathrm{\Delta }t$ is a Markovian birth, that is, $\lambda \mathrm{\Delta }t+o(\mathrm{\Delta }t)$ .

Example

Assume we have a record of 180 infection cases a month. We need to determine (a) the probability of eight incidents per day, (b) the probability of at most five cases in a day

• (a)There are $\lambda =180/30=6$ cases per day on average, which is the rate of infection per day. Then,

$$P\{x=8\}=\frac{{\lambda }^x}{x!}{e}^{-\lambda }=\frac{6^8}{8!}{e}^{-8}=0.103,$$

• (b)  $\begin{array}{c}P\{x\ge 4\}=1-(P\{x=0\}+\dots +P\{x=4\})\hfill \\ =1-(\frac{6^0}{0!}{e}^{-6}+\dots +\frac{6^4}{4!}{e}^{-6})=0.725.\hfill \end{array}$

The results of a simulation with the Poisson distributed infection threshold and branching convoluted with Gamma distribution are plotted in Figure 3.

Figure 3.

Illustration of rates for the spread of infection evolving as a branching process.

As can be noted, recovery starts soon after the infection starts. This enforces a small diminish in branching rates modeled as Gamma function. Gamma distribution is tied with the immunity growth that enforces annihilation of the spread, with the degree of immunity. The curve in Figure 3(b) shows the branching rate changes with an insignificant incubation time, whereas 3(c) shows the effect of latent incubation times.

Infection to quarantine

Assume that, following an epidemic outbreak, we have a group of arrivals of patients with a Poisson rate ${\lambda }_{\mathrm{iq}}$ . This process entails the discovery process, where a rapid diagnose is performed on each request in order to determine false positives. It is a fact that, the probability of no infection discovery during $t+\mathrm{\Delta }$ is

$$p_0(t+\mathrm{\Delta }t)=p_0(0,t)[1-{\lambda }_{\mathrm{iq}}\mathrm{\Delta }t],$$

a single discovery in the time interval $(t,t+\mathrm{\Delta }t)$ is

$$p_1(t,t+\mathrm{\Delta }t)={\lambda }_{\mathrm{iq}}\mathrm{\Delta }t,$$

the probability for more than one cases is expressed by

$$p_n(t+\mathrm{\Delta }t)=p_n(0,t)[1-{\lambda }_{\mathrm{iq}}\mathrm{\Delta }t]+p_{(n-1)}(0,t){\lambda }_{\mathrm{iq}}\mathrm{\Delta }t.$$

The diagnose process proceeds as follows. It will be shown that $p_0(t)=e^{-{\lambda }_{\mathrm{iq}}t}$ , and thus, the probability of $n$ arrivals during some time interval is proved by induction as:

$$P\{I_q=n|{\lambda }_{\mathrm{iq}}\}=\frac{{({\lambda }_{\mathrm{iq}}t)}^n}{n!}{e}^{-{\lambda }_{\mathrm{iq}}t}.$$

(16)

The infection–to–quarantine rate parameter, ${\lambda }_{\mathrm{iq}}$ , denotes the mean arrival rate to the quarantine service.

Let us now determine the distribution of the transition time from infection to the quarantine stage, which we model as a memoryless process, where let the state $s_0$ be the state at time $t=t_0$ and $\upsilon$ be the random time it takes the process to leave the state $s_0$ by going to some other state, then

$$\phi (t)=P\{\upsilon >t\}=e^{-{\lambda }_{\mathrm{iq}}t},t\ge 0,$$

(17)

where ${\lambda }_{\mathrm{iq}}$ denotes the density of the departure from state $s_0$ . If ${\lambda }_{\mathrm{iq}}=0$ , the process will remain forever in $s_0$ , however, if ${\lambda }_{\mathrm{iq}}>0$ , the probability will be ${\lambda }_{\mathrm{iq}}\mathrm{\Delta }t+o(\mathrm{\Delta }t)$ , denoting the probability of a change in the state during $\mathrm{\Delta }t$ . It follows from equation (17) that $P\{\upsilon >t\}$ is some function of $t$ . If, assume, $\phi (s)=P\{\upsilon >s\}$ , then the process will be in the same state at time $t_0+s$ as at time $t_0$ , and hence its subsequent behavior will be similar to that of $s=0$ . Especially,

$$\phi (t)=P\{\upsilon >s+t|\upsilon >s\}$$

(18)

will obtain the probability of event $\upsilon >s+t$ given that $\upsilon >s$ . Then,

$$\phi (t)\phi (s)=P\{\upsilon >s+t|\upsilon >s\}P\{\upsilon >s\}=P\{\upsilon >s+t\},$$

(19)

and hence,

$$\phi (s+t)=\phi (t)\phi (s).$$

(20)

Then, the random variable $\upsilon$ , called sojourn time in the current state, has the probability density given by

$$\begin{array}{c}{p}_{\upsilon }(t)=\frac{P[(\upsilon >s+t)\cap (\upsilon >s)]}{P(\upsilon >s)}=\frac{P(\upsilon >s+t)}{P(\upsilon >s)}\hfill \\ =\frac{e^{-{\lambda }_{\mathrm{iq}}(s+t)}}{e^{-{\lambda }_{\mathrm{iq}}s}}=e^{-{\lambda }_{\mathrm{iq}}t}=P\{\upsilon >t\},\hfill \end{array}$$

(21)

$$p_{\upsilon }(t)=\{\begin{array}{c}{\lambda }_{\mathrm{iq}}{e}^{-{\lambda }_{\mathrm{iq}}t},\mathrm{if}t\ge 1,\hfill \\ \multicolumn{1}{c}{0,\mathrm{if}t<0.}\end{array}$$

Thus, the infection to quarantine time is exponentially distributed.

Quarantine to recovery

As shown in the previous sub-section, nodes in the infected population (or network) are inspected, where false positives are discovered in significantly short times. Regarding the computer viruses, users generally discover malfunctioning of computers at random when they turn on the computer. A simple birth–death process is appropriate for modeling the discrimination of real infection events and false positives, where a successful diagnose of an infection and a false positive can be obtained by

$$\begin{array}{c}{\alpha }_k=\varphi (m-k)+\upsilon ,\hfill \\ \mathrm{birth}:p_{k,k+1}=[\varphi (m-k)]\frac{1}{{\alpha }_k},\hfill \\ \mathrm{death}:p_{k,k-j}=\left[\psi \left(\begin{array}{c}k\hfill \\ \multicolumn{1}{c}{j}\end{array}\right)q^j{(1-q)}^{k-j}\right]\frac{1}{{\alpha }_k}.\hfill \end{array}$$

(22)

Here, $q$ is the probability that a detection is indeed a false detection (death), which is considered as a diagnostic process consuming negligible (zero) time. The population size is $m$ , $\varphi$ denotes the rate of the infection discovery, that is, rate of births in $t$ time units, and $\psi$ denotes the inspection rate.

For the stationary state, the probability of exactly $j$ nodes being detected as infected is then obtained by

$$P_{\xi }\{\xi =j\}=p_j=\frac{m!}{(m-j)!}{\left(\frac{\varphi }{\psi }\right)}^j{p}_0,j=0,1,\dots ,m,$$

(23)

where, $p_0$ is determined from the condition

$$\begin{array}{c}\sum _{j=0}^m{p}_j=p_0\sum _{j=0}^m{\rho }^j=p_0\frac{1}{1-\rho )}=1,\hfill \\ p_0=1-\rho ,\rho =\varphi \psi .\hfill \end{array}$$

The average number of nodes discovered as infected and awaiting for recovery is given by

$$\mathbf{E}\left[I_k\right]=k-\frac{\varphi +\psi }{\varphi }(1-p_0)=k-\frac{\varphi +\psi }{\varphi }\left(\frac{\varphi }{\psi }\right).$$

(24)

Mortality rate

The probability $p_j$ from equation (23) incorporates the exponential infection rate $\varphi$ from the related Poisson distribution, where $\varphi =1-e^{-\varphi t}$ , then

$$n=\mathbf{E}[\xi (t)]=m{p}_j(t)=m(1-e^{-\varphi t})$$

as the average number of infected objects, whereas the number of lost objects is simply given by

$$\begin{array}{c}\overline{n}=m-\xi (t)=\mathbf{E}[m-\xi (t)]\hfill \\ =m-m{p}_j(t)=m{e}^{-\varphi t}.\hfill \end{array}$$

The reader is encouraged to determine the mortality in terms of $\mathbf{E}\left[I_k\right]$ from equation (24).

Recovery to healthy

Total recovery (or hosting) time for individuals is denoted by an exponentially distributed random variable, ${\mu }_{\mathrm{rh}}$ , where each individual waits for at least $T$ random time slots and eventually requires more time unit denoted by $\delta$ during the hosting period, then,

$$P\{T\ge \delta \}=e^{-{\mu }_{\mathrm{rh}}\delta },t\ge 0,$$

(25)

will give the probability for determining the total time required for the overall hosting of a patient. Let the times between various hostings be denoted by a Poisson process, $H$ , with parameter ${\mu }_{\mathrm{rh}}$ denoting the rate of recovery of an individual. Thus, $H$ has the probability density function

$$f_h(t)={\mu }_{\mathrm{rh}}{e}^{-{\mu }_{\mathrm{rh}}t}.$$

Hence, the mean time between recoveries is given by

$$E(H)=\frac{1}{{\mu }_{\mathrm{rh}}}.$$

The recovery process described by equation (25) can be simply modeled with two states, idle $s_0$ and busy $s_1$ , respectively. Assume a recovery facility is in an idle state at time $t_0$ , and due to the memoryless property of the Poisson arrivals, the probability $p_{0,1}(\mathrm{\Delta }t)$ of the system going from state $s_0$ to state $s_1$ during $\mathrm{\Delta }t$ is just the probability as described earlier by $\lambda \mathrm{\Delta }t+o(\mathrm{\Delta }t)$ . As can be noted, the density of the transition from $s_0$ to $s_1$ is denoted by $\lambda$ . However, assume the system is in the busy state $s_1$ (processing an input request) at time $t$ , then the probability of a system moving from the busy state ( $s_1$ ) to the idle state ( $s_0$ ) is $p_{1,0}(t)$ . Assuming the recovery service has been busy for $\upsilon$ time units at time $t$ , then the probability for the idle state of the recovery unit is

$$p_{1,0}(t)=1-P\{T>\upsilon +t|T>\upsilon \}=1-\frac{P\{T>\upsilon +t\}}{P\{T>\upsilon \}}.$$

(26)

Hence, plugging equation (25) into equation (26), we get the simplified equation for the transition probability from busy to idle state as

$$p_{1,0}(t)=1-\frac{e^{-{\mu }_{\mathrm{rh}}(\upsilon +t)}}{e^{-{\mu }_{\mathrm{rh}}\upsilon }}=1-e^{-{\mu }_{\mathrm{rh}}\upsilon }={\mu }_{\mathrm{rh}}\mathrm{\Delta }t+o(\mathrm{\Delta })t.$$

(27)

Obviously, we have a Markov process satisfying the following conditions

$$p_{0,1}(t)=1-p_{0,0}(t),p_{1,0}(t)=1-p_{1,1}(t).$$

(28)

We have the facts

$$\left(\begin{array}{cc}{\lambda }_{0,0}=-\lambda \hfill & {\lambda }_{0,1}=\lambda \hfill \\ \multicolumn{1}{c}{{\lambda }_{1,0}=\mu }& {\lambda }_{1,1}=-\mu \hfill \end{array}\right),$$

(29)

and

$$p_{i,j}(0)=\{\begin{array}{c}1,\mathrm{if}j=i,\hfill \\ \multicolumn{1}{c}{0,\mathrm{if}j\ne i}\end{array}$$

(30)

to our disposal. Transition density from state $s_i$ is denoted by ${\lambda }_i$ and from state $s_i$ to state $s_j$ is denoted by ${\lambda }_{i,j}$ , respectively. As can be noted service times for recovery are exponentially distributed with parameter ${\mu }_{\mathrm{rh}}$ and requests arrive as Poisson distributed expressed as

$$\begin{array}{c}{p}_{0,1}(\mathrm{\Delta }t)=1-e^{-\lambda \mathrm{\Delta }t}=\lambda \mathrm{\Delta }t+o(\mathrm{\Delta }t),\hfill \\ p_{1,0}(\mathrm{\Delta }t)=1-e^{-\mu \mathrm{\Delta }t}=\mu \mathrm{\Delta }t+o(\mathrm{\Delta }t).\hfill \end{array}$$

(31)

Recovery facility with single serviceman

We need to determine the limiting solutions so that we can easily determine the mean requests awaiting for the recovery at any time. Suppose we have a service facility with a single repairman (server) handling at most $m$ identical incoming requests. We assume also the incoming requests are Poisson distributed with mean arrival rate $\lambda$ , where each request has probability $\lambda \mathrm{\Delta }t+o(\mathrm{\Delta }t)$ and their arrivals are independent of each other. Thus, the recovery process is modeled as a Markov process, where recovering a request is exponentially distributed with parameter $\mu$ . To find the limiting probabilities for the recovery process we consider the following relations associated to transition densities both for the arrivals and departures of the requests.

$$\begin{array}{c}\underset{t\to 0}{lim}\frac{p_{j,j+1}(t)}{t}=(m-j)\lambda ,\mathrm{for}j=0,1,\dots ,m,\\ \underset{t\to 0}{lim}\frac{p_{j,j-1}(t)}{t}=\{\begin{array}{c}\mu ,\mathrm{for}j=1,2,\dots ,m,\hfill \\ \multicolumn{1}{c}{0,\mathrm{for}j=0,}\end{array}\\ \underset{t\to 0}{lim}\frac{1-p_{\mathrm{jj}}(t)}{t}=\{\begin{array}{c}(m-j)\lambda +\mu ,\mathrm{for}j=1,2,\dots ,m,\hfill \\ \multicolumn{1}{c}{m\lambda ,\mathrm{for}j=0.}\end{array}\end{array}$$

We can obtain the limiting probabilities by solving the following system of equations

$$\begin{array}{c}m\lambda p_0=\mu p_1,\hfill \\ \mu p_m=\lambda p_{m-1},\mathrm{for}j=m,\hfill \\ [(k)\lambda +\mu ]p_j=(k+1)\lambda p_{j-1}+\mu p_{j+1},\hfill \\ \mathrm{for}k=m-j,j=1,\dots ,m-1.\hfill \end{array}$$

(32)

Hence, solving this system of equations, the limiting probability of exactly $j$ requests present in the system is obtained by

$$p_j=\frac{m!}{(m-j)!}{\left(\frac{\lambda }{\mu }\right)}^j{p}_0,j=1,2\dots ,m.$$

(33)

Thus, based on the condition ${\sum }_{j=0}^m{p}_j=1$ ,

$$p_0=\frac{1}{{\sum }_{j=0}^m[\frac{m!}{(m-j)!}{(\frac{\lambda }{\mu })}^j]}.$$

We can now determine the limiting probability of exactly $j$ requests being recovered by solving the following system of equations.

$$\begin{array}{c}m\lambda p_0=\mu p_1,\hfill \\ m\mu p_m=\lambda p_{m-1},\mathrm{for}j=m,\hfill \\ [(k)\lambda +j\mu ]p_j=(k+1)\lambda p_{j-1}+(j+1)\mu p_{j+1},\hfill \\ \mathrm{for}k=m-j,j=1,2\dots ,m-1.\hfill \end{array}$$

(34)

Thus, the limiting probability $p_j$ leading to exactly $j$ recoveries out of $m$ requests is

$$p_j=\left(\begin{array}{c}m\hfill \\ \multicolumn{1}{c}{j}\end{array}\right){\left(\frac{\mu }{\lambda +\mu }\right)}^{m-j}{\left(\frac{\lambda }{\lambda +\mu }\right)}^j,\mathrm{for}j=0,1,\dots ,m.$$

(35)

It should be noted that these analyses assume $m$ identical requests each handled independently of the others.

Example: Resource utilization

Let us assume that a healthcare center randomly receives an average of 128 patients each 24 h. The center spends on average 8 min on each patient, where each individual is treated one at a time. Based on this information, we can determine a set of useful results regarding (a) proportion of the time the center is busy, (b) system utilization, and (c) mean service rate, (d) probability of the number of patients at the center, (f) idle time of the system.

Assume that the care center is modeled as an $M/M/1$ queue system with Poisson arrivals of patients and exponentially distributed service (diagnose + care) times. First, average arrival rate of patients per hour is

$$\lambda =128/24=5.333/h.$$

Average service rate per hour is

$$\mu =(8/60{)}^{-1}=(2/15{)}^{-1}=7.500/h.$$

System utilization (proportion of busy time) is

$$\rho =\lambda /\mu =5.333/7.500=0.7111/h.$$

Probability of having the system idle (i.e. exactly 0 patient under care) is

$$p_0=(1-\rho )=(1-0.7111)=0.289.$$

Probability of having exactly one patient under treatment is

$$p_1=(1-\rho ){\rho }^1=(1-0.7111)0.7111=0.205.$$

Probability of having more than one patient under treatment is

$$\begin{array}{c}{p}_{(n>1)}=1-p_{(n\le 1)}=(1-p_0)\rho =1-(p_0+p_1)\hfill \\ =1-(0.289+0.205)=0.506.\hfill \end{array}$$

(36)

Recovery facility with random service rates

Throughput of a service facility is one of the dominating factors for jobs with large and randomly fluctuating backlogs. It is also important for us to introduce this case here. Generally, higher workload may lead to a situation of “panicking” server resulting in diminished throughput, whereas lower workload may divert system to a “lazy” server, also resulting in diminished throughput. It is often suitable to assume constant mean service rates for systems involving automated or machinated services. However, in reality, queueing systems involving human servers cannot always have constant rates. For example, increased incoming requests (e.g. patient flow into an emergency unit in a hospital) may result in increased pressure on the doctors/nurses undertaking the service. This will, in turn, result in increased effort and higher service rate, or even cause a worsened situation due to the higher pressure. On the other hand, when the backlog (queue length) is small, then the service rate will decrease, due to the “lazy” server syndrome. More on the related issues can be found in Gross et al., ⁵⁹ Stewart, ⁶⁰ and Asmussen. ⁶¹ One can also refer to Schwarz ⁶² for more detailed analysis of time-dependent queueing systems.

Suppose we have a system with state space $\xi [0,1,2,\dots ,n]$ behaving as a continuous-time birth–death process having birth rates ${\lambda }_j$ and death rates ${\mu }_j$ , where $j$ is the current number of requests present in the system. The transition probabilities

$$p_{i,j}(t)=P\{\xi (t)=j|\xi (0)=i\},(i,j)\in \xi [0,1,2,\dots ]$$

need to satisfy the forward Kolmogorov equations

$$\begin{array}{c}\frac{\partial p_{i,0}(t)}{\mathrm{dt}}=-{\lambda }_0{p}_{i,0}(t)+{\mu }_1{p}_{i,1}(t)\hfill \\ \frac{\partial p_{i,j}(t)}{\mathrm{dt}}={\lambda }_{j-1}{p}_{i,j-1}(t)-({\lambda }_j+{\mu }_j)p_{i,j}(t)+{\mu }_{j+1}{p}_{i,j+1}(t).\hfill \end{array}$$

We then have for a system having stationary characteristics

$$\underset{t\to \infty }{lim}{p}_{i,n}(t)=p_n=P\{\xi =n\}.$$

Hence, we have the limiting distribution $\forall p_i,i=0,1,2\dots$ , and as detailed in Appendix, we have

$$p_n=p_0\underset{i=1}{\overset{n}{\Pi }}\frac{{\lambda }_{i-1}}{{\mu }_i},$$

where

$$p_0={[1+\sum _{k=1}^{\infty }\underset{i=1}{\overset{k}{\Pi }}\frac{{\lambda }_{i-1}}{{\mu }_i}]}^{-1}.$$

On the basis of the non-catastrophic constrain we also have the following relation.

$$\left(\frac{p_i}{p_{i-1}}\right)=\left(\frac{{\lambda }_{i-1}}{{\mu }_i}\right),p_i=\left(\frac{{\lambda }_{i-1}}{{\mu }_i}\right)p_{i-1},\mathrm{for}i=0,1,2,\dots$$

Again, assuming Poisson arrivals with $\lambda$ and exponential service times denoted by ${\mu }_n$ . A new parameter ( $R$ ) is now needed to represent the degree to which the server is affected due to the large backlog, where the queue length is $n$ . Therefore, the mean service rate for $n$ requests will be

$${\mu }_n=n^R{\mu }_1,R=0,\dots ,1,n=1,2,\dots ,$$

(37)

where, $R$ represents a random coefficient denoting the stress factor and ${\mu }_1$ stands for the expected service time when there is only one request in the queue. Hence, the overall throughput ratio of a queue with $n$ requests will be

$${\sum }_n^R=\frac{1}{{(n!)}^R}{\left(\frac{\lambda }{{\mu }_1}\right)}^n,n=1,2,\dots ,.$$

(38)

The negative effect of large backlogs leads to decrease in the service rate. This can be due to the annoyance of individuals waiting longer than expected, or diverting the requests to other facilities, if available. Such cases will typically lead to decreased arrivals, where a random coefficient, $D$ , affects the growth rate of the input queue given in

$${\lambda }_n^D={\lambda }_0(n+1{)}^{-D}n=1,2,\dots ,.$$

(39)

This random fluctuating with large backlogs is illustrated in Figure 4, which depicts the result of a simulation with random mean service rates satisfying the relation of $\rho \le \frac{\lambda }{\mu }$ .

Figure 4.

The effect of fluctuating arrivals causing increased service times: (a) input queue, (b) increased service time with random fluctuations affected by large backlog, and (c) decreased steady-state departure rime.

Experiments and simulations

One may need to experiment a chosen model with different cases and scenarios in order to validate it, that is, to determine whether it can produce expected results. Experiments should be better combined with simulations so that they can be numerically verified based on various scenarios. It is usually trivial to simulate computer epidemics due to more concrete definition and fewer number of the system parameters compared to that of human epidemics. In fact, some human viral infections behave similarly to that of computer infections, and can be easily adapted to computer simulations. For example, some human/animal viruses attack certain cells and proliferate within the attacked victims in the same way as a computer worm spreads itself. Therefore, we can modify some validated computer virus scenarios to match the biological viral infection cases and scenarios. The scenarios can then be experimented with in a network laboratory to observe the spread and control of infections, and the results are then justified with computer simulations. The simulation can then be expanded to evaluate a variety of biological epidemics, if desired. For example, CMV: cytomegalovirus, HPV: Human Papilloma, HTLV: Human T-Lymphotropic, rubella, influenza, Herpes 1 and Herpes 2 are the family of human viruses, where their propagation patterns can be justified by computer simulations, if appropriate models with sufficient amount of data are available.

In order to emphasize the diversity (or similarity) between the human and computer epidemics, a simulation of two related cases was conducted, where both processes have Poisson input, but their recovery processes evolved differently. It was shown that the human recovery process involves the Gamma distributed recovery process. This result is realistic, due to the automatically grown immunity found in human intelligence. ⁶³ Figure 5 illustrates results of the simulations with constant infection rate of $\lambda =10$ distributed by a Poisson Distribution. A large susceptible population size with minimal probability of infection can be modeled by a Poisson process, whereas higher probability and smaller population size will approximate to Binomial distribution.

Figure 5.

Simulation results depicting the infection and recovery characteristics of human and computer cases each having 1000 hosts.

The recovery is exponentially convolved with Gamma distribution. Gamma convolution increases the accuracy of recovery due to increased immunity of the infected objects while undergoing recovery. Based on a simulation, the changes in a population of 1000 nodes are shown in Figure 5.

Treatment evaluation

The efficiency of a healing method (or protection) can be evaluated to determine whether it is reliable during the course of its application. The model given here is straightforward, where the probability of a single failure and success of a treatment develops as a simple birth–death Markov chain. It is assumed that healing fails with rate $\lambda$ and succeeds with rate $\mu$ , where both rates are constant. Constant rate means that a given event occurs with a certain frequency during a specified time unit. For example, drug A cures 65% of patients of age 60–65 after 14 dosages. Thus, we can define two mutually exclusive states (success and fail) for the healing method at hand. Based on the definitions for the state transitions defined earlier, we can obtain the probability for the success or failure of a treatment. Let us first assume the following definitions for two time variables, $t$ and $\mathrm{\Delta }t$ :

• $P_s(t)$ : probability of success,

• $P_f(t)$ : probability of failure,

• $\lambda \mathrm{\Delta }t$ : probability of failure in $\mathrm{\Delta }t$ , that is, probability of transition from state success to failure in $\mathrm{\Delta }t$ .

• $\mu \mathrm{\Delta }t$ : probability of success in $\mathrm{\Delta }t$ , that is, probability of transition from state failure to success in $\mathrm{\Delta }t$ ,

• $1-\lambda \mathrm{\Delta }t$ : probability of not failing in $\mathrm{\Delta }t$ ,

• $1-\lambda \mathrm{\Delta }t$ : probability of not succeeding in $\mathrm{\Delta }t$ .

Hence, we can deduce the following set of equations to describe the probability of state transitions occurring during $t+\mathrm{\Delta }t$ .

$$\begin{array}{c}{P}_s(t+\mathrm{\Delta }t)=P_s(t)[1-\lambda \mathrm{\Delta }t]+P_f(t)\mu \mathrm{\Delta }t\hfill \\ P_f(t+\mathrm{\Delta }t)=P_f(t)[1-\mu \mathrm{\Delta }t]+P_s(t)\lambda \mathrm{\Delta }t\hfill \end{array}$$

(40)

These linear equations deal with a single item, which can be in any of two mutually exclusive Markovian states at a time, either failure (birth) or success (death). Hence, failure of a treatment is depicted by a birth process, whereas a success is depicted by a death process. Prior to solving these equations we rewrite them in the form of differential equations as.

$$\begin{array}{c}\frac{\partial P_s(t)}{\mathrm{dt}}=-\lambda P_s(t)+\mu P_f(t)\hfill \\ \frac{\partial P_f(t)}{\mathrm{dt}}=-\mu P_f(t)+\lambda P_s(t).\hfill \end{array}$$

Solutions of these equations can be obtained by applying various methods, for example, Laplace transform, recursive numerical method, Langranges, and transition matrix. Thus, by setting $P_s(0)=1$ at $t=0$ and applying Laplace transform, we can obtain a general solution to express the probability of failure $P_f(t)$ and success $P_s(t)$ , respectively. Setting $\psi =\lambda +\mu$ , we have

$$\mathrm{Success}:P_s(t)=\frac{\mu }{\psi }+\frac{\lambda }{\psi }{e}^{-\psi t}$$

$$\mathrm{Failure}:P_f(t)=1-P_s(t)=\frac{\lambda }{\psi }-\frac{\lambda }{\psi }{e}^{-\psi t}.$$

(41)

Predator–prey cases

In computer networks, we can have predator software that can fight against malware in real-time. Internet pop-up and Trojan blockers should work in this way. However, most anti-virus programs are usually configured to run periodically or manually, which then cannot be considered as dynamic predators. Because computer malware (worm, virus) are frequently and randomly deployed in different forms, discovering them in real-time becomes extremely difficult. The predator-prey model has a periodic behavior because the effect of hunting and surviving changes relative to the respective population growth rates. Lotka-Volterra predator-prey model has been frequently referred to for describing the relation between a predator and its prey, where the predator consumes the prey to survive. Increasing the level of consumption leads to higher starvation rates, especially if a certain prey is the only source of survival for many different predator types.

Let us simulate a scenario dealing with computer viruses in order to observe the power of predators behaving as a benign computer virus. The following differential equation, equation (42), models an anti-virus as the prey, where a computer virus plays the role of the predator, which is also programmed to destroy (stop) as many as anti-viruses as possible. Population parameters are $<a,A>$ and $<b,B>$ for the predators and preys, respectively. Each time the predator destroys an anti-virus (prey), the predator population grows following the Poisson distribution with the density parameter $\lambda$ . However, the prey has an immunity growth feature against the predator attacks, see Attack-obstacle model in Kondakci. ⁶³ Preys gain intelligence that proliferates them with a certain rate each time they survive an attack. Hence, the immunity growth of the prey follows the Gamma distribution with $\alpha$ and $\beta$ as the shape and scale parameters, respectively.

$$\begin{array}{c}\frac{\partial a}{\mathrm{dt}}=a[1-\frac{b}{B}]\left[\frac{{\lambda }^x}{x!}{e}^{-\lambda }\right]\hfill \\ \frac{\partial b}{\mathrm{dt}}=-b[1-\frac{a}{A}]\left[\frac{x^{\alpha -1}}{{\beta }^{\alpha }\Gamma (\alpha )}{e}^{-x/\beta }\right].\hfill \end{array}$$

(42)

Figure 6 illustrates the simulation results obtained from the model given in equation (42).

Figure 6.

Illustration of the predator-prey relation.

As shown in the figure, the rate changes have also random characteristics, because virus attacks are random, that is, by nature hunting and surviving occur in a random fashion. Ideally, the system can come to an equilibrium point, where the effect of infection rate will be governed by the prey population as shown in Figure 7.

Figure 7.

Illustration of the steady–state in the predator-prey extinction state.

Discussion

Designing accurate models can significantly mitigate the number of cases and mortalities. Models developed for computer epidemics cannot be directly applied to human epidemics, and vice versa. Computer epidemics differ from human epidemics in some aspects. For example, incubation (latent) time for computer viruses can be extremely short, and even in some cases unpredictably long. On the other hand, latent time depends generally on the type of virus for human epidemics, and depends minimally on the individual’s immunity.

We generally consider analysis of epidemics under two different numerical approaches, deterministic, and stochastic. Deterministic models can be simple to apply, in which the limiting distributions can be obtained more easily with fewer parameters. On the other hand, stochastic models have more sophisticated stages of analysis, which allow fine-graining of parameters for more precise results. This feature leads naturally to more time needed for analysis in order to achieve the equilibrium state, hence manual computation of results is infeasible for large data set, for example, more than 10 observations. Another limitation is the level of learning threshold, which may hinder its application in some environments.

Designing of an extended stochastic model is planned as a future work, in which the extension will deal with new births (connections coming from external networks), mortality, and cases of pandemics dealing with infectious agents exported to external networks.

Conclusions

It is important to build models reflecting realistic characteristics of a given case, including effective resource planning, estimation of future risks, and reliable measures for mitigation. This paper presented a generic model for analyzing the epidemic stages and their interaction with each other. In most cases, epidemiology is a stochastic process involving several stages of evolution. Therefore, applying queueing theory enhanced with stochastic analysis can give a better overview over the dynamics of epidemic cases. Although the correctness of other models (e.g. deterministic, diffusion, differential analysis) is undeniable, features of stochastic queue models can guide us to simpler and more structured solutions that can be realized using recursive computer algorithms and other numerical approaches. This paper has modeled all stages of an epidemic case, where each stage has been formulated thoroughly within a stochastic context, so that the approach can serve as a guide for building models applicable to different environments with stochastic behavior.

The approaches given here are applicable to both human and computer epidemics. Test and verification of models used for human epidemics can be validated by an associated computer model, where computes can be used as test objects representing the individuals of a human population. By this manner, we can acquire some key parameters that can be precisely used for the human model. For example, parameters such as immunity, incubation time, and infection rate can be adjusted to match the characteristics of a given virus attacking humans, and then run simulation with computers to obtain results reflecting the human epidemics caused by that virus. In order to increase the accuracy of a model we can set up various infection sceneries and configure a wide range of parameters and network sizes/types in a computer simulation. Hence, we can precisely obtain the characteristics of a given type of viral outbreak in a human population by fine-adjusting of the variables used for the simulation with computers as victims.