Stochastic SIS Modelling: Coinfection of Two Pathogens in Two-Host Communities.

A pathogen can infect multiple hosts. For example, zoonotic diseases like rabies often colonize both humans and animals. Meanwhile, a single host can sometimes be infected with many pathogens, such as malaria and meningitis. Therefore, we studied two susceptible classes S 1 ( t ) and S 2 ( t ) , each of which can be infected when interacting with two different infectious groups I 1 ( t ) and I 2 ( t ) . The stochastic models were formulated through the continuous time Markov chain (CTMC) along with their deterministic analogues. The statistics for the developed model were studied using the multi-type branching process. Since each epidemic class was assumed to transmit only its own type of pathogen, two reproduction numbers were obtained, in addition to the probability-generating functions of offspring. Thus, these, together with the mean number of infections, were used to estimate the probability of extinction. The initial population of infectious classes can influence their probability of extinction. Understanding the disease extinctions and outbreaks could result in rapid intervention by the management for effective control measures.

1. Introduction

Animal diseases such as rabies and hantavirus can be transmitted to humans. The pathogens involved therein cause zoonotic epidemics [1]. Since two or more hosts are included in the process, the spread of diseases is sometimes called a multi-host epidemic [2]. For example, rabies is commonly transmitted by domestic dogs [3], and wild rodents play a key role in transmitting hantavirus to humans. Meanwhile, the co-infection of two different pathogens in a single host is common [4]. This can be seen when a host is simultaneously infected by HIV and malaria [5] or by cholera and typhoid [6], for example. Even though many mathematical models describing the dynamics of this process were previously studied [7,8,9], models involving co-transmission and co-infection in multiple hosts are lacking. The probabilities of disease extinction associated with such models are not well-studied. Therefore, this study is concerned with extending the idea of the co-infection of two pathogens to two-host communities such as the zoonotic diseases described above.

In epidemiological studies, in order to capture the infection dynamics of a pathogen spreading into a given host community, two approaches—stochastic and deterministic techniques—are mostly used. As reported in [10,11,12], both of these have their advantages and disadvantages. Though in most scenarios the two techniques are used to analyze a similar problem, their formulations differ. The stochastic models are often derived through continuous time Markov chains (CTMCs); for example, some well-known results obtained by this technique have been used to investigate a variety of problems, such as host–vector malaria in susceptible–infective–susceptible (SIS) and susceptible–infective–recovered (SIR) epidemic models [13]. Meanwhile, the deterministic models formulated through ordinary differential equations (ODEs) can be good at approximating the overall epidemic dynamics when the population size is large [14]; however, when the population is small, stochastic models can be more appropriate to capture the randomness associated with the system [15]. The epidemic extinction can be estimated even if the population is growing exponentially [16]. The differences between stochastic and deterministic models are sometimes determined numerically [17]. When these two models were compared, the stochastic model was able to show that the time until epidemic extinction was longer [18]. Stochastic models have been reported to address some important epidemic questions through parameter estimation on the available data [19].

Another important class of stochastic model is the branching process. For example, a multi-type branching process can be applied to a wide range of epidemic models since it can approximate the CTMC models when answering epidemic questions. Exact solutions to problems of interest are not always available; therefore, asymptotic results derived through the branching process can serve as an alternative means [20]. This can be seen when the asymptotic formulae for the epidemic extinction [21,22], duration of epidemic outbreak [23], and the mean number of infectious individuals [24] were obtained by approximating the CTMC with the multi-type branching process.

A result, the epidemic outbreak is linked with the basic reproduction number [13]. This is an important ratio in epidemiological studies. The management’s intervention, as well as control measures, are concerned with how the epidemic can be minimized within the shortest possible time from the outbreak. This can be achieved when the basic reproduction number is successfully reduced to a threshold value less than one [12]. An epidemic outbreak model incorporating this ratio was used to investigate the stochastic patches, through which the movement of infectious individuals from higher-risk areas to those with lower-risk can serve as another effective control strategy [17]. Similarly, using the basic reproduction number, a stochastic model of two spreading pathogens was studied, through which the effects of key parameters responsible for spreading resistant bacteria in hospitals were determined [25]. The threshold for predicting epidemic extinctions in a single infectious class is , where is the basic reproduction number, and is the initial population count of the class [23]. However, this asymptotic result cannot hold for multiple infectious classes [12].

Furthermore, discrete models such as difference equations, which are often used in modeling species with non-overlapping generations (e.g., insects) [26], are not used in this work since the pathogens herein are assumed to be infecting both humans and animals. Therefore, the stochastic models are used to investigate the co-infection of two infectious classes due to their robustness in computing disease extinction. This can be achieved by approximating the multi-type branching process through the transition probabilities. The stochastic model used in this study is of SIS type. The environmental fluctuations of this model were reported to suppress the disease outbreak [27,28,29], determine the coexistence of two infectious diseases [30], and estimate the length of the epidemic disease outbreak [31]. In an attempt to estimate the probability of extinction of two infectious classes co-infecting two-host communities, we derived the deterministic system of ODEs using the transition probabilities of the CTMC models, determined the epidemic extinction and the mean number of infections using asymptotic formulae derived through the multi-type branching process technique, and compared the stochastic models and their deterministic counterparts through numerical simulations.

2. Materials and Methods

2.1. Continuous Time Markov Chain Model

Suppose that and are population sizes of two infectious classes infecting two susceptible host communities and . This can be written as the continuous time Markov chain [32], {, where t is the continuous parameter of the process. Assuming the probability of the population count of the four classes, taking is whereby Writing the schematic reactions of the given system allows us to derive the deterministic models through the transition probabilities of the CTMC.

The rate at which the two infectious classes, and , can be recovered are and respectively,

and .

The two susceptible classes, and , can move to infectious classes, and , at the disease transmission rates , and , such that

,

,

,

.

Taking the infinitesimal time , we list the events happening in the interval , together with their corresponding transition rates [33] (see Table 1).

Table 1

Assumptions of the continuous time Markov chain (CTMC) model.

Event	Transition Between t and t+δt	Probability
(i) Mortality of I_1(t)	(n_1,n_2+1,n_3,n_4)→(n_1,n_2,n_3,n_4)	Φ_1(n_2+1)δt+0(δt)^2
(ii) Mortality of I_2(t)	(n_1,n_2,n_3,n_4+1)→(n_1,n_2,n_3,n_4)	Φ_2(n_4+1)δt+0(δt)^2
(iii)I_1(t) infects S_1(t)	(n_1+1,n_2−1,n_3,n_4)→(n_1,n_2,n_3,n_4)	β_11(n_1+1)(n_2−1)δt+0(δt)^2
(iv)I_2(t) infects S_1(t)	(n_1+1,n_2,n_3,n_4−1)→(n_1,n_2,n_3,n_4)	β_12(n_1+1)(n_4−1)δt+0(δt)^2
(v)I_1(t) infects S_2(t)	(n_1,n_2−1,n_3+1,n_4)→(n_1,n_2,n_3,n_4)	β_21(n_2−1)(n_3+1)δt+0(δt)^2
(vi)I_2(t) infects S_2(t)	(n_1,n_2,n_3+1,n_4−1)→(n_1,n_2,n_3,n_4)	β_22(n_3+1)(n_4−1)δt+0(δt)^2

Substituting the probabilities given in Table 1, we get the following master equation known as the Kolmogorov forward differential equation: The master Equation (2) can be simplified with the aid of the following probability-generating function: Taking the derivative of Equation (3) with respect to t and substituting Equation (2) into the result obtained, we get the following partial differential equation, Using the moment-generating function [34], Equation (4) can be written as follows: Given that , we obtain from (5), for , and . This gives the following equations:

The heuristic approach allows us to write these equations as follows: The compartmental differential Equations (6) describe the dynamics of two infectious classes and mixing with two susceptible groups and . The birth and immigration of individuals are not captured by the model.

If the recovery rates, and , are not considered in the process, we obtain the following master equation: Applying similar techniques to Equation (7), we get the following system of ordinary differential equations: Thus, Equations (8) describe the dynamics of two infectious classes excluding the recovery terms in the disease compartments.

2.2. Basic Reproduction Number

As introduced earlier, the basic reproduction number is an important threshold for measuring how diseases are spreading in communities. When determining the basic reproduction number of the derived model (6), we used the matrix , which is referred to as the next-generation matrix. While F represents the matrix of infection rates, V is that of the transfer rates in the disease compartments. These can be analytically obtained by linearizing the system of ordinary differential equations near the disease free-equilibrium [21,35] as follows: and

The spectral radius of the next-generation matrix gives the following results [4,36]: and

Therefore, the basic reproduction number associated with the system can be written as

2.3. Multi-Type Branching Process

To examine the probability of epidemic extinction for the system of ODE Equations (6), we employed the multi-type branching process. An important theorem for determining such probabilities along with other statistics is stated as follows.

Suppose that the function whereby whereby

This theorem guarantees that we can determine the fixed points of equations of offspring as well as their mean numbers through the multi-step branching process. This process could be used to approximate our CTMC model. We assumed that represents the four classes of discrete random variables . Two classes of interest, through which the epidemic extinction and outbreak can be estimated, are and .

The probability of the epidemic extinction for the two classes and can be computed by taking the probability-generating function of the two random variables and , Considering the transition rates for the CTMC model in Table 1, we notice that the individual of type 1 can only give birth to its own type. This is similar to the individual of type 2.

The mean number of infections can be determined using the formula [24]: as

The probability of epidemic extinction for the system of ODEs (6) can be obtained since matrix M is consistent with the basic reproduction number (11). For the sub-critical case, the probability can be investigated by applying the jury condition [37], which states that the spectral radius of matrix M, if and only if It is clear from matrix M, and det. Thus, the jury condition holds for matrix M, and this also satisfies the first part of Theorem 1.

Since matrix M is reducible, the epidemic extinction in the super-critical case can be determined by considering the two basic reproduction numbers and . Solving for and in (12),(13), we get , , and . While the epidemic transmission and its recovery rates for are and , respectively, those of can be written as and , respectively. Thus, the fixed points for the probability-generating functions of the offspring, and the two non-zero entries and of matrix M can be used as follows:

The epidemic extinction in the super-critical case can be determined if , , the fixed point is Similarly, if , , we apply the fixed point

These satisfy the second part of Theorem 1. Since each infectious class transmits only its own type of pathogen, the determination of inverses for the two basic reproduction numbers (Equations (14), (15)) allows us to write the probabilities of extinction for the two infectious classes separately. Thus, for , we obtained the following [16,23]:

Meanwhile, the epidemic extinction of the second class for can be written as

A minor outbreak or epidemic extinction occurs if , , whereas a major outbreak can exist only if , [13].

3. Results

Numerical Examples

Since it can be difficult to analytically solve the transition probabilities from the forward Kolmogorov Equations (2) and (7), the Gillespie realizations conveniently give their numerical approximations. This can be achieved by drawing two pseudo-random numbers between zero and one from the uniform distribution. While the first random number allows jumps to the next states, the second is the inter-event time representing the time elapsing between events in the process [38]. Each stochastic realization is obtained with the successive number of these steps.

The stochastic model was found to agree with its deterministic counterpart. Using arbitrary parameters, the disease transmission rates, , , the infectious recovery rates , the initial number of susceptible individuals , , and the initial number of infectious classes, , , one Gillespie realization from each class approximately converges with their deterministic analogues (see Figure 1a,c). Though any other suitable values can be used to test the dynamics of the two infectious classes, the infection rates were held constant to allow us to examine the changes in and when varying the initial conditions of the two host communities. Thus, reducing the population of susceptible individuals affects the number of each infectious class (Figure 1). Removing the recovery terms in disease compartments (Equations (8)), we observed the disease persistence in the system (see Figure 1b,c).

Figure 1

One stochastic realisation approximates the deterministic models for , , , , given (a) the sample paths of Equations (6) with , , , , and ; (b) the sample paths of Equations (8) with the initial conditions , , , and ; (c) the sample paths of Equations (6) for , , , , and ; and (d) the sample paths of Equations (8) with the initial conditions , , , and .

To examine the behavior of the stochastic realizations after a long period of time, the distribution of the process can serve as an alternative solution. Using the parameters described above, , , , , , , one can see that at , the population of infectious individuals is nearly extinct (Figure 2). This indicates that both Figure 1 and Figure 2 capture the dynamics of the co-infection models, Equations (6) and (8).

Figure 2

Solutions of the compartments, (a) , (b) , (c) , (d) , of the epidemic model (6) at based on 1000 stochastic realisations. The parameters used for determining the distributions are , , , , , , , , and .

To determine the probability of extinction of each infectious class, we employed both the asymptotic formulae P and P, as well as the CTMC simulations. Though not all the entries of Table 2 were the same for both P, P and CTMC simulations (Approx 1 and Approx 2), they reveal that the extinction of each infectious class depended on its initial population size, when their numbers were small [21]. This is because the probability of disease extinction decreases with the increase in the initial population size of each of the two infectious classes (see Table 2).

Table 2

The probabilities of extinction P and P of two infectious classes and , while Approx 1 and 2 are their respective numerical estimations using CTMC for , = = 0.3 based on sample paths.

I_1(0)	I_2(0)	P_01	Approx 1	P_02	Approx 2
1	1	0.7500	0.4997	0.7500	0.5003
2	1	0.5625	0.3333	0.7500	0.6667
1	2	0.7500	0.6664	0.5625	0.3336
2	2	0.5625	0.5001	0.5625	0.4999
3	2	0.4219	0.3995	0.5625	0.6004
2	3	0.5625	0.6000	0.4219	0.4000
3	3	0.5000	0.5000	0.4219	0.4219

4. Discussion

In this study, the probability of epidemic extinction of two infectious classes co-infecting two-host communities was estimated. The basic reproduction number, a threshold for determining epidemic extinction, can be applied to a single infectious class to study problems related to probabilities of extinctions and outbreaks [23]. However, when many infectious classes are to be studied, their probabilities of extinction can be estimated through the matrix of mean number of infections obtained from the multi-type branching process [16]. When determining this result in this work, we reported a reducible matrix. Therefore, no single explicit formula could be used to estimate the extinctions of the two infectious classes related to our formulation (Equation (6)).

Furthermore, since the transition probabilities for our CTMC model (Table 1) were used when formulating the probability-generating functions of the offspring for the multi-type branching process, our study reported the emergence of two basic reproduction numbers, and , for the two infectious classes and , respectively. This is consistent with what has been obtained using the next-generation matrix in this study and others [4,12,36]. Using the two basic reproduction numbers and CTMC simulations, the epidemic extinctions of two infectious classes which infect two-host communities was estimated. Our results show that the CTMC simulations could nearly approximate the asymptotic results involving reproduction numbers when many sample paths were used (i.e., ). This is similar to what was reported in [16]. Our results equally show that the extinction of each infectious class depended on the initial number of individuals therein (Table 2).

The Gillespie realizations were reported to approximate the deterministic model (see Figure 1). This is similar to what has been reported in many epidemiological studies comparing stochastic and deterministic models, such as [11,18]. When comparing the solutions of models (6) and (8), our results showed that the disease persistence could be suppressed by introducing recovery terms in the disease compartments.

5. Conclusions

The formulation and analyses of the stochastic models of two infectious classes co-infecting two different susceptible groups can give additional information on how pathogens spread in host communities. Estimation of the disease extinction through the multi-type branching process supplements our results obtained by the deterministic model, such as the basic reproduction number.

Specifically, the stochastic formulation of SIS models provides essential information on the extinctions of two pathogens spreading in two host communities. The transition rates of the CTMC not only produce the deterministic system of ODEs through the master equations, but also approximate those (transition rates) of the multi-type branching process. This provides information about the disease extinction using a small number of initial population sizes of the infectious classes. Both the stochastic SIS simulations and the solutions of ODEs suggest that removing the recovery terms in the disease compartments can lead to the persistence of diseases.