The network-level reproduction number and extinction threshold for vector-borne diseases

The reproduction number of deterministic models is an essential quantity to predict whether an epidemic will spread or die out. Thresholds for disease extinction contribute crucial knowledge on disease control, elimination, and mitigation of infectious diseases. Relationships between the basic reproduction numbers of two network-based ordinary differential equation vector-host models, and extinction thresholds of corresponding continuous-time Markov chain models are derived under some assumptions. Numerical simulation results for malaria and Rift Valley fever transmission on heterogeneous networks are in agreement with analytical results without any assumptions, reinforcing the relationships may always exist and proposing a mathematical problem of proving their existences in general. Moreover, numerical simulations show that the reproduction number is not monotonically increasing or decreasing with the extinction threshold. Key parameters in predicting uncertainty of extinction thresholds are identified using Latin Hypercube Sampling/Partial Rank Correlation Coefficient. Consistent trends of extinction probability observed through numerical simulations provide novel insights into mitigation strategies to increase the disease extinction probability. Research findings may improve understandings of thresholds for disease persistence in order to control vector-borne diseases.


Introduction
Vector-borne diseases greatly impact health of humans and animals and are among the leading causes of worldwide death every year [16]; almost half of the world's population is infected with at least one type of vector-borne diseases and millions of people die of vector-borne diseases every year [9]. These diseases also cause significant economic losses in areas of animal trade, agriculture, health care, tourism, as well as destroy ecosystems of throughout the world. Therefore, control and prevention of vector-borne diseases are both economical and humane. Efficient interventions require a good understanding of disease transmission and persistence, and dynamic modeling of vector-borne diseases may contribute greatly to this end [14]. A model may be used to learn many characteristics of an outbreak such as: whether or not an outbreak may occur, the size of the outbreak, the duration time of the outbreak, or the probability for the epidemic to die out [6]. Efficient mitigation strategies deduced from model results may stop an outbreak at early stages by reducing spreading parameters [6].
Globalization of trade and travel is one of the key factors driving the emergence of vectorborne diseases; heterogeneous structure also plays an important role on dynamics of infectious diseases [20]. Modeling the spatial spread of vector-borne diseases is a challenging task [3], but one possible approach is to consider a meta-population as a directed graph, or a network, with each vertex representing a subpopulation in a location, and links placed between two locations if there is a possibility of transmission, such as movement or proximity [5]. Network models are more widely used in epidemiology to understand the spread of infectious diseases through connected populations [27,39].
The basic reproduction number, R 0 defined as the number of secondary cases produced by an infected individual in a naive population [12] is an important threshold on epidemiology among many others, such as type reproduction number [32], target reproduction [34], and threshold index for epidemicity [18]. The basic reproduction number is an important metric, predicting whether a disease will spread or die out in deterministic population and communicable disease theory [2]. If R 0 > 1, one infectious individual generally produces more than one infection, leading to spread of an epidemic, whereas if R 0 < 1, one infectious individual generates less than one infection on average [10], and epidemic may die out [12]. The same trajectory can always be observed with deterministic models given the same initial conditions [21]. If it is possible for an epidemic to occur again, a real world epidemic does not allow us to observe that the same infection happens to the same person at the same time [21]. Moreover, deterministic models have the shortcoming that the number of infected individuals may go to less than one [23].
In comparison, Markov chain models are more realistic in the sense of only taking integer values instead of continuously varying quantities [23] and taking into account chances by approximating or mimicking the random or probabilistic factors. The last infectious individual may recover before the infection is transmitted to other susceptible individuals so that the disease may become extinct [23]. Consequently, an infection introduced to a completely susceptible population may not invade the system even if R 0 > 1 [23]. Threshold for the extinction of an infectious disease to occur and probability of disease extinction are of interests. Bienaymé-Galton-Watson branching processes are widely used to study extinction of diseases involving multi-type infections.
Lloyd [23] reviewed theory of branching processes and computed extinction probability using branching processes for Ross malaria model [33] taking into account stochasticity and heterogeneity. Pénisson [28] presented several statistical tools to study extinction of populations consisted of different types of individuals, and their behaviors before extinction and in the case of a very late extinction. Allen and Lahodny Jr [1] computed reproduction numbers for deterministic models, and extinction thresholds for corresponding continuous-time Markov chain (CTMC) models using continuous-time branching process, and derived their relationships. A CTMC model is a stochastic counterpart of a deterministic ordinary differential equation (ODE) model [1]. Lahodny Jr and Allen [22] estimated probability of disease extinction for a Susceptible-Infected-Susceptible (SIS) multipatch model and illustrated some differences between thresholds for deterministic models and stochastic models numerically. Allen and van den Driessche [2] established connections between extinction thresholds for continuous-time models and discrete-time models and illustrated the relations through numerical simulations. Although probability of disease extinction is defined as the probability for the number of infections to become zero when time goes to infinity, various numerical approximations for many types of models within finite time showed good agreement with predicted extinction probability using branching processes [1,2,22].
Deriving relationships between reproduction numbers and extinction thresholds is a complex task for vector-borne diseases transmitted on heterogeneous networks due to too many parameters and large size matrices. According to current knowledge, very little research has studied it. The objectives of this research are to relate the extinction threshold, E 0 in a stochastic setting and the reproduction number, R 0 in a deterministic setting for vector-host meta-population models, as well as gain understandings in how to increase extinction probability.
The contribution of our work is summarized as follows.
1. Relationships between extinction thresholds and the reproduction numbers are derived for network-based vector-host models under some assumptions.
2. Numerical simulations show that the relationships still exist after removing above assumptions.
3. Consistent trends of extinction probability varying with disease parameters are observed through extensive numerical simulations.
4. The key parameters in predicting uncertainty of the extinction threshold are identified using Latin Hypercube Sampling/Partial Rank Correlation Coefficient (LHS/PRCC).
5. The relationship between varying disease parameters and potential mitigation strategies is biologically interpreted.
This paper is organized as follows. Section 2 reviews the next generation matrix approach for computing R 0 and the branching process for deriving E 0 . Section 3 calculates R 0 for a deterministic vector-host model in which transmission dynamics of vectors are described by a Susceptible-Infected (SI) model and transmission dynamics of hosts are described by an SIS model. We relate E 0 of corresponding CTMC model and R 0 analytically. In Section 4, an analogue of results in Section 3 is obtained for a model in which transmission dynamics of vectors are described by a Susceptible-Exposed-Infected (SEI) model and transmission dynamics of hosts are described by a Susceptible-Exposed-Infected-Recovered (SEIR) model. Local transmission and trans-location transmission due to proximity for vector-borne diseases are both considered in Sections 3 and 4. In Section 5, the relationships derived in Sections 3 and 4 are numerically shown to hold without any assumptions for simplified malaria and Rift Valley fever meta-population models. The sensitivity test sorted out the key parameters in predicting uncertainty of extinction probability. Relationships between varying parameters and extinction probabilities are explored through extensive simulations for homogeneous populations and a two-node network. Section 6 provides a summary and discussion of mathematical derivations and simulation results.

Preliminary
The next generation matrix approach used to compute R 0 for compartmental models is reviewed here, followed by a review of the multitype branching process approximation used to derive E 0 for corresponding CTMC models.

Computation of R 0 using the next generation matrix approach
The next generation matrix approach is frequently used to compute R 0 . In this section, we quickly review this approach. For more details, we refer to [11,Chapter 5], [38]. For simplicity, let Y i , i = 1, · · · , m stand for compartments that are only related to infected and asymptomatically infected individuals. The original nonlinear system of ODEs including these compartments can be written as ∂Y i ∂t = F − V , where F = (F i ) and V = (V i ) represent new infections and transfer between compartments, respectively. Moreover, F i represents the rate at which new infections appear in compartment i, and represents the rate at which individuals transfer from (resp. into) compartment i. The Jacobian matrices F representing transmission, and V representing transition are defined as: where x 0 denotes disease free equilibrium (DFE), and x j is the number or proportion of infected individuals in compartment j, where j = 1, · · · , m. Matrix F is nonnegative and V is a nonsingular M-matrix. Matrix F V −1 is called the next generation matrix. The (i, k) entry of F V −1 indicates the expected number of new infections in compartment i produced by the infected individual originally introduced into compartment k, where i, k = 1, · · · , m.
The reproduction number, R 0 , is defined as the spectral radius of F V −1 , denoted by ρ(F V −1 ).

Deriving E 0 using branching process approximation
Calculating the probability of disease extinction is one of the most interesting applications of branching process. The branching process may lead to disease extinction or persistence. We are interested in the conditions under which a disease may become extinct and the probability for this event to occur. First, we review the approach of using branching process to compute extinction threshold and extinction probability for multi-type infections. We refer to [1,28] for the rest of this section. Let − → X (t) = (X 1 (t), · · · , X n (t)) T : t ∈ (0, ∞) be a set of discrete-valued vector random variables. Assume that individuals of type i produce individuals of type j and that the number of infected individuals produced by type i are independent of the number of infected individuals produced by other individuals of type i or type j for i, j = 1, · · · , n, i = j. Additionally, individuals of type i have identical probability generating function (pgf). Let {X ji } n j=1 be the offspring random variables for type i, where X ji is the number of infected individuals of type j produced by individuals of type i. The probability that one individual of type i produces x j infected individuals of type j is given as The offspring pgf array (g 1 , · · · , g n ) : [0, 1] n → [0, 1] n , is defined as Note that a trivial fixed point of (g 1 , · · · , g n ) always exists at 1 = (1, · · · , 1). We denote by M = [m ij ] n×n the expectation matrix of offspring distribution which is nonnegative, where m ij := ∂g i ∂w j | x=1 < ∞ represents the expected number of new infected individuals of type j produced by an individual of type i.
The extinction threshold, E 0 is defined as the spectral radius of the expectation matrix, denoted by ρ(M ).
Here, a function is called simple if it is linear with no constant term.

SI vector model and SIS host metapopulation model
In this section, a deterministic vector-host model in which disease transmission dynamics of vectors are described by an SI model, while transmission dynamics of hosts are described by an SIS model. The reproduction number and extinction threshold for corresponding CTMC model are analytically related.

The reproduction number
The model for vectors consists of compartment G representing susceptible vectors, and compartment J representing infected vectors. Disease dynamics of hosts are modeled by an SIS model.
The recruitment rate of vectors (resp. hosts) in node i is η i (resp. ψ i ) for all i = 1, · · · , n. The rate of new infections in vectors in node i produced by local hosts, and hosts in other nodes are β i G i I i /N i and n j=1,j =i ω ji G i I j /N j , respectively. The death rate of susceptible and infected vectors in node i are µG i and µJ i , respectively. The rate of host infection in node i produced by local vectors, and vectors in other nodes are α i S i J i /N i and n j=1,j =i σ ji S i J j /N i , respectively. The death rates of susceptible and infected hosts in node i are d i S i and d i I i , respectively. The rate of recovery for hosts in node i is γ i I i .
Since J i and I i , i = 1, · · · n are only compartments related to infected and asymptomatically infected, system of ODEs (4) can be rewritten as follows.
To calculate R 0 , we first prove the following lemma.

The threshold for extinction probability
In this section, we compute E 0 for corresponding CTMC of model (4). See Table 1 for state transitions and rates. The pgfs are:  (4) omitting node index i.
where j = 1, · · · , n, the index k = i − n for n + 1 ≤ i ≤ 2n, w i represents I V i = 1, I H i = 0, and u i represents The expectation matrix M is: where A, B are the same as those in (5), and Note that if both A and B are positive matrices, then the assumptions (B 0 ) and (B 1 ) in [28] hold for this model. Proof.
Proposition 2. The extinction threshold of model (4) satisfies that

The relationship between R 0 and E 0
To obtain a theoretical relationship between R 0 in (8) and E 0 , we assume that for constant numbers k 1 , k 2 ∈ [0, 1] throughout this section. The assumption can be interpreted biologically as: the probability of natural death is identical for vectors from each node, and the probability of natural death is identical for hosts from each node. The assumption shall be removed for numerical simulations in the next section.
Corollary 1. If the further assumption is made that k 1 = k 2 except assumption (11), then

The reproduction number
The following model extends the model in Section 3.1 by adding compartment Z for exposed vectors, and compartment E for exposed hosts. Other terms have identical meanings as corresponding ones in model (4). The rate at which exposed vectors and exposed hosts in node i transfer to infected compartments are ϕ i Z i and ε i E i , respectively.
Compartments related to infected and asymptomatically infected are Z i , E i , J i and I i , i = 1, · · · , n. The unique solution at DFE is (G 0 i , 0, 0, N 0 i , 0, 0, 0), where G 0 i and N 0 i are the same as those in Section 3.1. The above system of ODEs including these compartments can be rewritten as follows. d dt Z 1 · · · Z n E 1 · · · E n J 1 · · · J n I 1 · · · I n T = F − V .
By a direct calculation, Following Lemma 1, (13) is

The threshold for extinction probability
State transitions and rates for corresponding CTMC of model (13) are listed in Table 2. The pgfs are: where w i represents only Z i = 1, w i+n represents E i = 1, u i represents J i = 1, and u i+n represents I i = 1 for i = 1, · · · , n. The indexes k = i − n for 1 ≤ i ≤ n, p = i − 2n for n + 1 ≤ i ≤ 2n, and q = i − 3n for 3n + 1 ≤ i ≤ 4n. The expectation matrix M is: Similarly, the assumptions (B 0 ) and (B 1 ) in [28] ).

The relationship between R 0 and E 0
In this section, the assumption (11) holds and k 1 < k 2 . Under the assumption (11), by Lemma 2, Recall that, for any square matrices A, B with the same size, ρ(AB) = ρ(BA). By this property, By (14), (15) and (16), Hence, Similarly, we have the following theorem.

Theorem 2.
Under assumption (11), Corollary 2. If a further assumption is made that k 1 = k 2 besides assumption (11), then Proof. The proof is similar to that of Corollary 1.

Numerical results
We show numerically the general relations between R 0 and E 0 for two models on heterogeneous networks. Significant parameters in predicting the uncertainly in E 0 are found by Latin Hypercube Sampling/Partial Rank Correlation Coefficient (LHS/PRCC) analysis. Finally, the trends of parameters varying with extinction array is summarized.

Numerical results on relations between R 0 and E 0
Model (4) is applied to study thresholds for malaria transmission through numerical simulations. Five thousand realizations with parameters uniformly distributed within the ranges listed in Table 3 on a four-node network give rise to R 0 ranging from 0.7668 to 63.8111 and E 0 from 0.8965 to 1.9140. The ranges of R 0 and E 0 vary with the number of nodes of a network and the assumed ranges of vector (host) recruitment rates with ranges of other parameters fixed. The values of R 0 are sorted from small values to large valued in Figure 1(a) and 1(b), and E 0 are ranked from small values to large values in Figure 1(c) and 1(d). The largest value of E 0 is 0.9980 when all values of R 0 are smaller than 1 and R 0 ≤ E 0 , as shown in Figure 1(a). The smallest value of E 0 is 1.003 when all values of R 0 are greater than 1 and R 0 ≥ E 0 , as shown in Figure 1(b). The largest value of R 0 is 0.9947 when all values of E 0 are smaller than 1, as shown in Figure 1(c). The smallest value of R 0 is 1.006 when all values of E 0 are greater than 1, as shown in Figure 1(d). The value of E 0 is not monotonically increasing with the increase of R 0 , as shown Figure 1(a) and 1(b). Similarly, R 0 fluctuates as E 0 increases, as shown in Figure 1(c) and 1(d).
Model (13) is applied to numerical examine the relationship between R 0 and E 0 for Rift Valley fever. See Table 4 for descriptions and ranges of parameters. Five thousand realizations produce R 0 ranging between 0.2289 and 54.5086 and E 0 from 0.6757 to 1.9763. The values of R 0 are ordered from small to large magnitudes in Figure 2(a) and 2(b), and the values of E 0 are from small to large values in Figure 2(c) and 2(d). The largest value of E 0 is 1 when all values of R 0 are smaller than 1, and √ R 0 ≤ E 0 , as shown in Figure 2(a). The smallest value of E 0 is 1.005 when all values of R 0 are greater than 1, and √ R 0 ≥ E 0 , as shown in Figure 2(         largest value of R 0 is 0.9998 when all values of E 0 are smaller than 1, and √ R 0 ≥ E 0 , as shown in Figure 2(c). The smallest value of R 0 is 1.008 when all values of E 0 are greater than 1, and √ R 0 ≥ E 0 , as shown in Figure 2(d). When R 0 increases, E 0 does not always increase, as shown in Figure 2(a) and 2(b). Similarly, R 0 fluctuates as E 0 increases, as shown in Figure 2(c) and 2(d).

Sensitivity analysis
We employ Latin Hypercube Sampling/Partial Rank Correlation Coefficient (LHS/PRCC) analysis [25] to identify key parameters whose uncertainty contribute to predict uncertainty of E 0 for model (13) and rank the parameters by their significances. The parameters shown to be significant with large PRCC values (> 0.5) or small p-values (< 0.05) [15] by the sensitivity test with 5000 sets of parameter values, are listed in Table 5. The magnitude of PRCC value represent the contribution to the prediction for the imprecision of E 0 , and a negative sign indicates that the parameter is inversely proportional to the magnitude of E 0 . The closer PRCC value is to +1 or −1, the more greatly the parameter impacts the outcome of E 0 .

Trends of extinction array with varying parameters
Consistent trends of w * are observed by numerical simulations for homogeneous populations and a two-node network for Model (13). Table 6 lists three different values for each parameter and corresponding extinction array for homogeneous populations as an example. Table 7 shows the trends of extinction array by varying one parameter at a time, keeping other parameters fixed and E 0 > 1 for homogeneous populations and a two-node network. If at least one entry of extinction array increases and others remain constant, then we say that the array increases. The extinction array w * decreases with the increase of contact rates from local vectors and vectors in other nodes to local hosts, contact rates from local hosts and hosts in other nodes to local vectors, death rates of hosts, recruitment rates of vectors, and incubation rates of vectors and hosts, whereas increases with the increase of vector death rates, vector recovery rates, host recruitment rates.
Changing parameter  Increasing parameter (w * 1 , · · · , w * n , u * 1 , · · · , u * n ) Table 7: Summary of trends for extinction array changing with one parameter at a time, while keeping other parameters fixed and E 0 > 1 for model (13) for homogeneous populations and a two-node network throughout various simulations.
The first part was proven by Allen and van den Driessche under the assumption (16) in [2], i.e., (F − V ) T = W (M − I), where F and V are Jacobian matrices defined in (1), M is a mean matrix of offspring distribution defined in Section 2.2, I is the identity matrix, and W is a positive diagonal matrix with each entry w i representing the rate parameter at which lifespan of group i are exponentially distributed for i = 1, · · · n [28]. This assumption does not hold for both model (4) and model (13).
Consistent trends in the extinction array w * while changing one parameter through numerical simulations is helpful in deducing trends of extinction probability and possible interventions for vector-borne diseases. According to Equation (3), the probability of disease extinction is monotonically increasing (decreasing) with the increase (decrease) of the extinction array when the initial number of infection is fixed. The following biological interpretations on disease extinction or persistence are in terms of fixed initial number of infections. If contact rates from vectors to hosts (α, σ), or those from hosts to vectors (β, ω) increase, the probability for the disease to persist is higher. If death rates of hosts (d) increase, the number of vectors is relatively dominant. Consequently, the disease is more likely to persist. Similarly, growing vector recruitment rates (η) increase the probability for disease persistence. The higher incubation rates in vectors (ϕ) or that in hosts (ǫ) lead to faster vector or host infections, such that the disease is prone to persist. On the contrary, increasing death rates of vectors (µ) may reduce rates of host infection, and ultimately, may increase the likelihood of disease extinction. Increasing recovery rates of hosts (γ) may reduce the number of infections, such that the probability of disease extinction increases. Increasing recruitment rate of hosts (ψ) may reduce vector infection rates and increase probability of disease extinction.
The findings show that we may increase extinction probability of vector-borne diseases by properly controlling vector and host population size, and promptly detecting and applying treatment for hosts. Analytical and numerical results shed light on deriving relationships between R 0 and E 0 , as well as connections between varying parameters and increasing extinction probabilities for many other vector-borne diseases transmitted among heterogeneous works. In summary, the resulting mathematical derivations and numerical simulations facilitate understanding thresholds for the spread of vector-borne diseases, as well as providing novel insights into disease prevention, mitigation and control strategies.