Bacteriophages—a new hope or a huge problem in the food industry

1.   Introduction

Appropriate mathematical models of infectious diseases are helpful to reveal the mechanism of disease transmission and to predict the final scale of the epidemic ^[1]. The social contact network of a population plays a significant role in controlling the propagation of infectious diseases ^[2,3]. However, directly analyze the propagation of disease on stochastic contact networks is very difficult, one often relies on deterministic mean-field models that are aimed at approximating some average quantities derived from the stochastic models ^[4,5,6,7].

Based on the probability generating function for the degree distribution on a contact network, Volz introduced an edge-based SIR epidemic model, which has the great advantage of having a system with only three differential equations ^[8]. This method describes the SIR disease spread process with heterogeneous connectivity by the contact network structure. Its core is to track the change of the probability that the neighboring nodes around the susceptible nodes are infected. Furthermore, Miller derived a single differential equation with only a single higher order term that governs the dynamics of the disease ^[9]. Recently, using the EBCM approach for the spread of susceptible-infected-recovered (SIR) diseases in contact networks, Miller et al. ^[10] derived some simple ordinary differential equation models capturing social heterogeneity and heterogeneous contact rates while explicitly considering the impact of partnership duration. However, the theoretical analysis of such models is usually difficult, and their validity is usually demonstrated through numerical simulations. For an edge-based epidemic model with large initial value, Miller ^[11] obtained the threshold of disease transmission and the final infection size. Using the same method, in a recent paper^[12], we established a sexually transmitted disease model on a bidirectional contact network with large initial value and analyzed in detail the global dynamics of the model. We refer the reader to references ^{[13,15,14,16]} and references therein for related works in this direction.

Vaccination in disease control strategies generally reduces the susceptibility of a node of network, but could either increase or decrease the infectiousness of a node by reducing the severity of symptoms. Modelling a disease with heterogeneous infectiousness and susceptibility in stochastic contact network can be seen in ^[17,18]. Using the EBCM technique, in a recent paper ^[19], Miller and Volz proposed an SIR model considering the heterogeneous susceptibility and infectiousness of individuals. Assume $ i $ is a parameter measuring a node's ability to become infected and cause infection, but that is does not influence the contact structure of the population. Without loss of generality, we assume that the population is divided into $ M $ types. The model proposed in ^[19] takes the following form:

with initial value $ \theta_{i, j}(0) = 1 $ and $ I(0) = R(0) = 0 $, where $ \tilde{\theta}_j = \sum_{l = 1}^M\theta_{j, l} \mathcal {Q}(l) $. The variables and parameters in model (1) are summarized in Table 1. For this model, the authors provided a numerical example to confirm its validity ^[19]. Here, we investigate the dynamics of the model (1) through mathematical analysis.

The organization of this paper is as follows. In the next section, we recall the derivation process of model (1), then we present our main results, including the derivation of the final size of an epidemic. In section 3, we provide the detailed proof of our main theorem. Finally, in section 4, we provide two simple scenarios for the special cases when $ M = 1 $ and $ M = 2 $ to make our results easier to understand.

2.   Derivation of model (1) and its global stability results

In this section, we recall the derivation process of model (1), then we present our main results.

2.1.   Recalling the derivation process of model (1)

Assume the network structure is unchanged, that is the degree distribution $ P(k) $ of the network is fixed. The nodes of the network are classified into $ M $ types according to their ability to become infected and cause infection. Now randomly choose a test node $ u $ from the network, and assume that the probability it is of type $ i $ is $ \mathcal{Q}(i) $, then $ \sum_{i = 1}^M\mathcal{Q}(i) = 1 $. It is easy to deduce the probability that the test node $ u $ is susceptible, infected or recovered. In fact,

where $ S_i $, $ I_i $ and $ R_i $ are respectively the probabilities that a type-$ i $ node is susceptible, infected and recovered at time $ t $, and therefore

Let $ \theta_{i, j} $ denote the probability that an edge from a type-$ j $ partner $ v $ to the test node $ u $ of type-$ i $ has not transmitted infection from $ v $ to $ u $, and let $ \phi_{S:i, j} $ be the probability $ v $ is still susceptible, $ \phi_{I:i, j} $ the probability $ v $ is infected but the edge has not transmitted, and $ \phi_{R:i, j} $ the probability that $ v $ has recovered without transmitting. Then $ \theta_{i, j} = \phi_{S:i, j}+\phi_{I:i, j}+\phi_{R:i, j} $, and therefore,

Define $ \tilde{\theta}_i $ as the probability that a random edge to a type-$ i $ node $ u $ has not transmitted infection to $ u $, then $ \tilde{\theta}_i = \sum_{j = 1}^M\theta_{i, j} \mathcal {Q}(j) $, and therefore,

and

where $ \psi(x) = \sum_k P(k)x^k $ is the probability generating function of degree distribution $ P(k) $. The transmission rate $ \beta_{i, j} $ from a type-$ j $ node to a type-$ i $ node and the recovery rate $ \gamma_i $ of a type-$ i $ node are type dependent, and they are assumed to be positive. From the probability fluxes between $ S_i $, $ I_i $ and $ R_i $ compartments (upper) and those between $ \phi_{S:i, j} $, $ \phi_{I:i, j} $ and $ \phi_{R:i.j} $ (lower) as shown in Figure 1, we have that

We can derive from the last two equations in Eq (6) with initial value

that $ \phi_{R_{i, j}} = \frac{\gamma_j}{\beta_{i, j}}(1-\theta_{i, j}) $. Combining Eqs (3) and (5) and the second equation of $ \phi_{R:i, j} $ in Eq (6), we then derive system (7) from the last equation of $ \theta_{i, j} $ in Eq (6) (see ^[19] for more details). Thus, the whole dynamics of the disease transmission are summarized by model (1).

2.2.   Main results

Notice that the dynamics of model (1) is determined by the dynamics of the following subsystem:

where $ i, j = 1, 2, \ldots, M $. Define

Then $ \Omega $ is a positive invariant set of system (7). We only need to consider the dynamics of system (7) constrained in $ \Omega $. Notice that system (7) always contains the disease free equilibrium $ E_0 = (1, 1, \dots, 1, 1, \ldots, 1) $.

Let $ D_{\theta}F(\theta) = (D_{\theta}F_{i, j})(\theta) $ be the jacobian matrix of system (7), where

and denote

where $ D_{\theta}F(E_0) $ is the jacobian matrix $ D_{\theta}F(\theta) $ at $ E_0 $. Our main result can be summarized in the following theorem.

Theorem 2.1. System (7) has at most two equilibria in the region $ \Omega $. Moreover,

(a) When $ \rho_0 < 0 $, then the system has only the disease free equilibrium

which is globally asymptotically stable in $ \Omega $.

(b) When $ \rho_0 > 0 $, then the disease free equilibrium $ E_0 $ becomes unstable, and another outbreak equilibrium

appears which is globally asymptotically stable in the region $ \Omega\setminus\{E_0\} $.

Remark 1. Based on Theorem Theorem 2.1, we know from model (1) that if $ \rho_0 > 0 $, then after the epidemic is over, the fraction of individuals who have never been infected is

where $ \tilde{\theta}_i^{(*)} = \sum_{j = 1}^M\theta_{i, j}^{(*)}\mathcal {Q}(j) $. Thus, the final size of an epidemic in the population is

3.   Proof of Theorem Theorem 2.1

We first introduce the following lemma which is used in the proof of our main theorem.

Lemma 3.1 (^[20]). For system $ \frac{d}{dt} x = F(x) $, let F(x) be a $ C^1 $ cooperative vector field in $ R^n $, whose flow $ \phi $ preserves $ R_+^n $ for $ t\geq0 $ and is irreducible in $ R_+^n $. Assume that the origin is an equilibrium and that all trajectories in $ R_+^n $ are bounded. Suppose the matrix-valued map $ DF $: $ R^n\rightarrow R^{n\times n} $ is strictly antimonotone, in the sense that if $ x > y $, then $ DF(x) < $DF(y). Then either all trajectories in $ R_+^n $ tend to the origin, or else there is a unique equilibrium $ p\in \text{Int}\ R_+^n $ and all trajectories in $ R_+^n\setminus \{0\} $ tend to $ p $.

Proof. The jacobian matrix $ D_{\theta}F(\theta) $ of system (7) is

where $ i, \ j, \ l, \ m = 1, 2\ldots, M $, and when $ i = l $,

when $ i\neq l $,

The off-diagonal elements of $ D_{\theta}F(\theta) $ are non-negative, and it is easy to check that $ D_{\theta}F(\theta) $ is strongly connected^* when $ \beta_{i, j} > 0 $. Thus, system (7) is a irreducible cooperative system in $ \Omega $ ^[20]. Consequently, when $ \rho_0 < 0 $, the equilibrium $ E_0 $ is locally asymptotically stable, and when $\rho_0>$0, the equilibrium $E_0$ is unstable.

^*The definition of strongly connected and irreducible cooperative systems, refer to ^[20].

To prove the global results of system (7) in $ \Omega $ applying Lemma 3.1, we make the coordinate transformation $ \vartheta_{i, j} = 1-\theta_{i, j} $, under which system (7) becomes

where $ i, j = 1, 2, \cdots, M $, and $ \tilde{\vartheta}_j = \sum_{l = 1}^M(1-\vartheta_{j, l})\mathcal {Q}(l) = 1-\sum_{l = 1}^M\vartheta_{j, l}\mathcal {Q}(l) $. Notice that $ \Omega $ is still the positive invariant set of system (10), and the disease free equilibrium $ E_0 $ is now transformed into the zero equilibrium $ \mathscr{E}_0 $ of system (10), where

Denote $ \mathcal {F}(\vartheta) = (\mathcal {F}_{1, 1}, \mathcal {F}_{1, 2}, \cdots, \mathcal {F}_{1, M}, \mathcal {F}_{2, 1}, \cdots, \mathcal {F}_{M, M})(\vartheta) $. Noticing that $ \mathcal{F}(\vartheta) = -F(1_{M^2}-\vartheta) $, we can deduce that

where $ D_{\vartheta}\mathcal{F}(\mathscr{E}_0) $ is the jacobian matrix $ D_{\vartheta}\mathcal{F}(\vartheta) $ of system (10) at $ \mathscr{E}_0 $. Moreover, from Lemma 3.1 we know that $ \mathscr{E}_0 $ is locally asymptotically stable when $ \rho(D_{\theta}F(E_0))\leq 0 $, and unstable when $ \rho(D_{\theta}(E_0)) > 0 $.

We first consider the global stability results of system (10). Notice the following:

(a) System (10) is still an irreducible cooperative system in the bounded region $ \Omega $ as

(b) $ D_{\vartheta}\mathcal {F}(\vartheta) $ is strictly antimonotone in $ \Omega $. Since $ \psi(\tilde{\theta}_i) $ and $ \psi''(\tilde{\theta}_i) $ are both monotone increasing functions with respect to $ \theta_{i, j} $ in the bounded region $ \Omega $, and therefore both $ \psi(\tilde{\vartheta}_i) $ and $ \psi''(\tilde{\vartheta}_i) $ are monotone decreasing functions of $ \vartheta_{i, j} $ in $ \Omega $ due to the relationship $ \tilde{\vartheta}_i = 1-\sum_{j = 1}^M\vartheta_{i, j}\mathcal {Q}(j) $. By comparison, we can check that when $ \vartheta_1 < \vartheta_2 $, we have $ D_{\vartheta}\mathcal {F}(\vartheta_1) > D_{\vartheta}\mathcal {F}(\vartheta_2) $.

According to Lemma 3.1, either all trajectories in $ \Omega $ tend to the zero equilibrium $ \mathscr{E}_0 $ of system (10), or else there is a unique equilibrium

and all trajectories in $ \Omega\setminus\mathscr{E}_0 $ tend to $ \mathscr{E}_* $. These, together with the local stability results of $ \mathscr{E}_0 $, indicates that when $ \rho_0 < 0 $, $ \mathscr{E}_0 $ is globally asymptotically stable in $ \Omega $, and when $ \rho_0 > 0 $, it becomes unstable and the equilibrium $ \mathscr{E}_* $ is globally asymptotically stable in $ \Omega\setminus\mathscr{E}_0 $.

Now let $ E_* = 1_{M^2}-\mathscr{E}_* $, then $ E_*\in \text{Int}\ \Omega $ is the unique outbreak equilibrium of system (7). Based on the arguments above, we deduce that when $ \rho_0 < 0 $, $ E_0 $ is globally asymptotically stable in $ \Omega $, and when $ \rho_0 > 0 $, it becomes unstable and the outbreak equilibrium $ E_* $ is globally asymptotically stable in $ \Omega\setminus\mathscr{E}_0 $. The proof of Theorem Theorem 2.1 is thus completed.

Remark 2. Theorem 2.1 provides a theoretical basis for the prevention and control of heterogeneous infectious and heterogeneous infectious diseases. For example, vaccination generally reduces the susceptibility of a node in community, then the population can be divided into those who have or have not received vaccination. The results of our analysis can be used to analyze the effectiveness of vaccination strategies.

4.   Two simple scenarios

Scenario 4.1. In the special case when $ M = 1 $, i.e., there is only one type of node in contact network, system (7) is reduced to the model established in Miller ^[9]:

Here $ E_0 = (1) $ and $ \rho_0 = -(\beta_{1, 1}+\gamma_1)+\beta_{1, 1}\frac{\psi''(1)}{\psi'(1)} $. Moreover, $ \rho_0 = 0 $ can be equivalently written as $ \frac{\beta_{1, 1}}{(\beta_{1, 1}+\gamma_1)}\frac{\psi''(1)}{\psi'(1)} = 1 $. Define

Then $ \mathcal{R}_0 $ is the basic reproduction number of the disease. Notice that $ \mathcal{R}_0\leq 1 $ if and only if $ \rho_0\leq 0 $. We have the following result.

Corollary 1. System (11) has at most two equilibria in the region

Moreover,

(a) When $ \mathcal{R}_0 \leq 1 $, then the system has only the disease free equilibrium $ E_0 $ which is globally asymptotically stable in $ \Omega $.

(b) When $ \mathcal{R}_0 > 1 $, then the disease free equilibrium $ E_0 $ is unstable, and the outbreak equilibrium $ E_* = (\theta_{1, 1}^{(*)})\in \text{Int} \Omega $ appears and it is globally asymptotically stable in the region $ \Omega\setminus\{E_0\} = \{\theta_{1, 1}: 0\leq \theta_{1, 1} < 1\} $.

Proof. We only need to prove that when $ \mathcal{R}_0 = 1 $ (i.e., $ \rho_0 = 0 $), $ E_0 = (1)_{1\times 1} $ is globally asymptotically stable. In fact, denote

We compute that

for $ \theta_{1, 1}\in[0, 1) $ due to $ \rho_0 = 0 $. Thus, $ F_{1, 1}(\theta_{1, 1}) $ is strictly monotonically decreasing on the interval $ [0, 1] $ with $ F_{1, 1}(1) = 0 $ and therefore $ E_0 = (1) $ is the unique equilibrium of system (11). It then follows from $ \frac{ {\bf{d}}\theta_{1, 1}}{ {\bf{d}}t} > 0 $ that for any solution with initial value $ \theta_{1, 1}(0)\in [0, 1) $, $ E_0 $ is globally asymptotically stable in $ \Omega = [0, 1] $. The proof is thus completed.

Scenario 4.2. In the special case when $ M = 2 $, i.e., there have two difference type of node in contact network, system (7) takes the form:

Here $ E_0 = (1, \; 1, \; 1, \; 1) $. The jacobian matrix $ D_{\theta}F(\theta) $ of system (13) at $ E_0 $ is

where $ \triangle_{i, i}(1) = (\beta_{i, i}+\gamma_i)-\beta_{i, i}\frac{\psi''(1)}{\psi'(1)}\mathcal {Q}(i), \ i = 1, 2 $. Define $ \rho_0 $ as in Eq (8). We have the following result.

Corollary 2. System (7) has at most two equilibria in the region

Moreover,

(a) When $ \rho_0 < 0 $, then the system has only the disease free equilibrium $ E_0 $ which is globally asymptotically stable in $ \Omega $.

(b) When $ \rho_0 > 0 $, then the disease free equilibrium $ E_0 $ is unstable, and the outbreak equilibrium

appears and it is globally asymptotically stable in the region $ \Omega\setminus\{E_0\} $.

Remark 3. From Scenario 1 we know that for the special case when $ M = 1 $, the disease free equilibrium $ E_0 = (1) $ is globally asymptotically stable (GAS) when $ \rho_0 = 0 $ (i.e., $ \mathcal{R}_0 = 1 $), but we can not prove that a similar result holds for the special case when $ M = 2 $. We guess the disease free equilibrium $ E_0 $ of model (1) is also GAS for all $ M\neq 1 $. We leave it as an open problem.

5.   Discussions and conclusions

Considering the fact that the individuals may have different susceptibility and infectiousness to the disease spreading in the population, in a recent paper ^[19], Miller and Volz proposed an SIR disease network model (i.e., model (1)) with heterogeneous infectiousness and susceptibility. The authors have provided a numerical example to demonstrate its validity but its mathematical analysis remain unsolved. In this paper, with the aid of the nature of irreducible cooperative system in the theory of monotonic dynamical system, we prove that the dynamics of the model are completely determined by a critical value $ \rho_0 $: When $ \rho_0 > 0 $, the disease persists in a globally stable outbreak equilibrium; while when $ \rho_0 < 0 $, the disease dies out in the population and the disease free equilibrium is globally stable.

Notice that it is assumed that in model (1) the random network considered is statically fixed and does not take into account the influence of new nodes, deleted nodes and other factors on the network. This is an approximation of the reality. More proper and reasonable model should consider these facts and be based on dynamic random networks. Moreover, how to characterize the epidemic process of heterogeneous infectious networks through numerical simulation is also a challenging problem. We leave all these for our future consideration.

Acknowledgments

Research is supported by the National Natural Science Foundation of China (No. 11671260).

Conflict of interest

The authors declare there is no conflict of interest.