Stability Criteria for a Multi-city Epidemic Model with Travel Delays and Infection during Travel

We present a compartmental SIR (susceptible–infected–recovered) model to describe the propagation of an infectious disease in a human population, when individuals travel between p different cities. The time needed for travel between any two locations is incorporated, and we assume that disease progression is possible during travel. The model is equivalent to an autonomous system of differential equations with multiple delays, and each delayed term is defined through a system of ordinary differential equations. We establish some necessary and sufficient conditions for the disease-free equilibrium of the model to be asymptotically stable.


Introduction
In recent years there has been an increasing interest in the mathematical modelling community for the spatial spread of infectious diseases.Historical examples like the 1918 influenza pandemic, as well as recent threats like the 2002-2003 SARS epidemic, the 2009 A(H1N1) influenza pandemic, the 2012 MERS coronavirus outbreak, and the 2015 extensive Ebola virus (EBOV) outbreak in West Africa exemplify that national boundaries or oceans have never prevented infectious diseases to reach distant territories.Differential equation-based models for spatial epidemic spread have been discussed for an array of infectious diseases including measles, influenza, malaria and SARS; most recent literature includes the work of Arino and coauthors [1,2], Gao and coauthors [7], Hsieh and coauthors [10], Li and coauthors [14], Ruan Email: d.knipl@ucl.ac.ukD. Knipl and coauthors [21], and Wang and coauthors [28,29].These studies are mainly concerned with the spatial dispersal of infected individuals in connected regions, and do not take it into account that long distance travel also provides a platform for disease dynamics.A recent report by the European Centre for Disease Prevention and Control [5] confirms that some infectious diseases, such as tuberculosis, measles and seasonal influenza, are transmissible during commercial flights, even on those with a duration of less than eight hours.There are some documented cases of influenza transmission occurring during train journeys, including when seven healthy people caught Tamiflu-resistant A(H1N1) influenza while travelling for 42 hours on a train in Vietnam (Mai et al. [17]).
In their disease transmission models, Cui and coauthors [3], Lui and Takeuchi [16], and Takeuchi and coauthors [25] examined the possibility that individuals may contract a disease while travelling.They used the SIS (susceptible-infected-susceptible) framework to describe disease progression, and implicitly used the assumption that transportation between regions takes place instantaneously.This approach ignores that the delay in the transportation can lead to a significant underestimation of the roles of global travel in the spread of some diseases of major public health concern, such as SARS and influenza.These diseases are rapidly progressive (typically, a few days), so even a short delay (a fraction of a day) can be significant.This consideration motivated the work of Liu, Wu and Zhou [15], where they introduced the time delay needed to complete the travel into the SIS-type epidemic model and considered the possibility of infections during this time.The global dynamics was described by Nakata [18] for two identical regions, and then by Nakata and Röst [19,20] for multiple regions with different characteristics and general travel networks.
For most diseases with a potential to pose a threat for a pandemic (like influenza, SARS and measles), the SIS-type model is inadequate.This was first noted by Knipl and Röst [12] when they suggested a general framework for SIR (susceptible-infected-recovered)-based disease progression with a general incidence term.This concept was further elaborated in an SEAIR (susceptible-exposed-asymptomatic infected-symptomatic infected-recovered)framework initiated by Knipl, Röst and Wu [13], closely mimicking the spatiotemporal evolution of the H1N1 pandemic in multiple regions connected by long distance travel.In these models the framework is required to include a subsystem structured by age, to incorporate the possibility of disease transmission during travel.In the subsystem, age represents the time elapsed since the start of the travel.Following the assumption that transmitting the infection is possible on-board, the model setup leads to a system of delay differential equations with delay representing travel time.The two systems, describing the dynamics in the regions and during transportation, are interconnected; initial values of the system for disease spread during travel depend on the state of the system in the regions, while the inflow term of arrivals to the regions after being in transportation for a fixed time arises as the solution of the subsystem for travel.If the subsystem can be solved analytically (this was the case for the age-structured SI model used in [15,18,20]), then the system for disease spread in the regions decouple from the subsystem.Recalling that initial data of the system during travel comes from the equations for the regions, the inflow term of travellers completing a trip appears as a delayed feedback term.On the other hand, in case of choosing SIR-type models as an epidemic building block when the subsystem does not admit a closed form solution, the delayed term in the system for the regions cannot be expressed explicitly (as exemplified in [12,13]), but is defined dynamically, via the solutions of another system.
The general form of initial value problems for nonautonomous functional differential equations with dynamically defined delayed feedback function were studied by Knipl [11].Fun-damental properties were obtained such as the usual existence, uniqueness and continuous dependence result for the solution, and some further, biologically relevant properties were also studied.Although the paper [11] provides mathematical framework for studying the SIR model in [12] and the SEAIR model in [13], it remains a major theoretical challenge to classify the global dynamics in such systems.In particular, no general method has been developed to our knowledge to investigate stability in delayed systems with dynamically defined delayed feedback function.To classify local stability of the disease-free steady state is extremely important as it characterizes whether the disease will become extinct or persist.This was attempted in [13] in the special case when the subsystem structured by travel time did obtain a closed form solution.In this paper, we present a general, delayed SIR model for the spread of an infectious disease when individuals travel between p cities.The model setup leads to a differential system with multiple delays, and transport-related infection is incorporated; more precisely the corresponding subsystems obtain SIR dynamics therefore cannot be solved analytically.We describe local stability properties of the disease-free equilibrium.The approach presented in the paper has the potential to be adopted to other delayed differential models where the delayed feedback term is defined via the solution of another system.

Model formulation
We describe the spread of an infectious disease in a population of individuals who travel between different locations.Considering p locations (cities), we divide the entire populations into the disjoint classes I j , S j , and R j , where the letters I, S, and R represent the compartments of infected, susceptible, and recovered individuals respectively, and lower index j (1 ≤ j ≤ p) specifies the current city.Let I j (t), S j (t), and R j (t) be the number of individuals belonging to I j , S j , and R j respectively, at time t.We assume constant recruitment terms A j , while d j denotes natural mortality in city j.Disease transmission is modelled by standard incidence and the transmission rate between an infected individual and a susceptible individual in city j is denoted by β j .Infected individuals in city j recover at rate µ j .For the total population currently being in city j at time t, we use the notation N j (t) = I j (t) + S j (t) + R j (t).
We denote by m k,j the movement rate from city j to city k (1 ≤ j, k ≤ p, j = k), furthermore, we let m j,j = 0 for j = 1, . . ., p.The total travel outflow rate from city j is given by M j = ∑ p k=1 m k,j , and the movement matrix M is defined as M = (m j,k ) p j,k=1 .We assume that M is an irreducible matrix.We introduce the terms I τ j,k , S τ j,k , and R τ j,k for the inflow of infected, susceptible, and recovered individuals respectively, into city j from k (j = k).These terms will be precisely described later.Based on the assumptions above, the following system is formulated for the dynamics of an infectious disease in p cities: (2.1)

D. Knipl
The time needed to complete travel from city k to j is denoted by τ j,k .It is assumed that τ j,k > 0 whenever m j,k > 0; however, for convenience we let τ j,k = 0 for all j, k such that m j,k = 0, that is, there is no movement from city k to j.
A principal assumption of our modelling approach is that we do not neglect the possibility of infection and recovery during the travel between the cities.To describe these processes during a travel to city j started from city k at any time t * , we divide the population during the travel into the classes i j,k , s j,k , and r j,k for infected, susceptible, and recovered travellers, respectively.Then i j,k (θ), s j,k (θ), and r j,k (θ) give the density of infected, susceptible, and recovered individuals respectively, who started travel at time t * , with respect to θ that denotes the time elapsed since the beginning of the travel.With and by denoting by β T j,k and µ T j,k the transmission rate and the recovery rate respectively, we consider the ODE system for the dynamics of the disease during a travel from city k to j: (2.2) This system can be obtained in a compact form as where y j,k = (i j,k , s j,k , r j,k ) and g j,k : R 3 → R 3 is defined by the right hand sides in (2.2).
Standard results from the theory of ordinary differential equations guarantee that for any initial data y * ∈ R 3 there is a unique solution to (2.3) on [0, ∞).We denote the solution by ỹj,k (θ) = ỹj,k (θ; 0, y * ).
Let h j,k : R 3 → R 3 , h j,k (v) = m j,k v, and let I j,k , S j,k , R j,k : R 3 → R such that (2.4) gives the number of infected, susceptible, and recovered individuals who leave city k at time t * to travel to city j.Therefore ỹj,k (θ; 0, h j,k (v)) describes the density of individuals in the three disease classes θ unites of time after the beginning of their travel from city k to j, that started at t * .In particular, for θ = τ j,k and using the notations in (2.4), (I j,k (v), S j,k (v), R j,k (v)) gives the density of infected, susceptible, and recovered individuals at the end of the travel.
We are now in the position to define the inflow terms I τ j,k , S τ j,k , R τ j,k in system (2.1).The individuals who arrive at time t to city j from k are precisely those who started their travel at time t − τ j,k ; hence, letting t * = t − τ j,k we derive by the above arguments that the density of infected, susceptible, and recovered individuals arriving at time t to city j from k is given as It is therefore meaningful to define the inflow terms as Note that these terms depend on some past states of the system and therefore they are delayed terms; however they are not expressed explicitly by the variables of the system (2.1) but are defined dynamically, through the solution of another system (2.2).
Let σ j = max{τ k,j : 1 ≤ k ≤ p} for j = 1, . . ., p.We define where for σ > 0 we denote by C([−σ, 0], R) the space of continuous functions on [−σ, 0], which is a Banach space with the usual uniform norm |φ| = sup{φ(θ) : −σ ≤ θ ≤ 0}.The generic element of C τ is denoted by φ = (φ 1 , . . ., φ 3p ), and C τ is a Banach space with the norm |φ| = ∑ 3p n=1 |φ n |.The system (2.1) in a compact form reads as where x = I 1 , . . ., I p , S 1 , . . ., S p , R 1 , . . ., R p , and f : . . ,f Rp ) where f I j , f Sj , f Rj : C τ → R are defined as the right hand side of the corresponding equation in system (2.1), for 1 ≤ j ≤ p.For the segment of the solution we use the usual notation x t where x t = (x t 1 , . . ., x t 3p ), and x t j = x j (t + θ), x t j+p = x j+p (t + θ) and The feasible phase space in our model is defined as the nonnegative cone C τ + of C τ .The general form of initial value problems for nonautonomous functional differential equations with dynamically defined delayed feedback function were studied in [11].Fundamental properties were obtained such as the usual existence, uniqueness and continuous dependence result for the solution, and some further, biologically relevant properties were also studied.

Stability analysis
In Appendix A we show that the system (2.1) has a unique disease-free equilibrium (DFE) Ê, where Ê ∈ C τ is the constant function equal to (0, . . ., 0, S1 , . . ., Sp , 0, . . ., 0) ∈ R 3p for all values of its argument and For 1 ≤ k ≤ p we denote by Êk the constant function equal to (0, Sk , 0) for all values of its argument.To study the stability of Ê, we investigate the linear variational system corresponding to Ê: for 1 ≤ j ≤ p, where we use the usual notation a • b for the dot product of two vectors a and b.In Appendix B we show that the gradients on the right hand side are obtained as The stability of the trivial solution of the linear variational system (3.1) is determined by the characteristic equation ∆(λ) = 0, which arises by seeking solutions of the form X(t) = e λt u where u ∈ R 3p .If the stability modulus, defined as is negative then the trivial solution of the linear variational system is locally asymptotically stable (LAS) whereas it is unstable if s(F) > 0. These stability properties extend to the stability of the DFE Ê in the system (2.1) by the principle of linearised stability.
Using the equations (3.2) and (3.3) it is possible to rewrite the linear variational system (3.1) as d dt and with all other elements of F j,k being zero.It is useful to note that for any j and k it holds that (F j,k ) u,v = 0 whenever u ≤ p and v > p.This implies that the first p equations (the infection subsystem) in the linear variational system are independent from the other equations.Moreover, we also note that for any j and k it holds that (F j,k ) u,v = 0 whenever u, v are such that either p + 1 ≤ u ≤ 2p and v > 2p, or u > 2p and p + 1 ≤ v ≤ 2p holds.It follows that the two systems, formed by equations p + 1, . . ., 2p (the susceptible subsystem) and by equations 2p + 1, . . ., 3p (the recovered subsystem) respectively, are also independent from one another.These remarks imply that the characteristic equation factorizes as the product of the characteristic equations of the three subsystems, hence the stability in (3.1) is determined by the stability properties in the three subsystems.In particular, both the susceptible and recovered subsystems can be obtained in the form and the infection subsystem reads We define the following ODE system, that can be associated with (3.5) by ignoring the delays: We present now one of our most important results, which states that the stability analysis of the DFE in the nonlinear model (2.1) with dynamically defined delayed feedback function reduces to studying the non-delayed infection subsystem of the linear variational system.The principal result of Section 5, Chapter 5 in [24] is that the stability of an equilibrium of a cooperative and irreducible system of delay differential equations is the same as that of an associated system of cooperative ordinary differential equations.To study these properties for system (3.5),we recall some definitions and results from [24] in Appendix C. Then we prove in Appendix D that the system (3.5) is cooperative and irreducible.
To study the stability of the linear system (3.6), we define and obtain the matrix of the linear system as T − Σ, where we used the short hand notation c j,k = e τ j,k (β T j,k −µ T j,k ) (j = k) and c j,j = 0.The stability of the zero solution is determined by the eigenvalues of the coefficient matrix; the solution is LAS if s(T − Σ) < 0 and it is unstable D. Knipl if s(T − Σ) > 0, where s(A) denotes the maximum real part of all eigenvalues of any square matrix A.
We say that a square matrix A is a non-singular M-matrix if it has the Z-sign pattern (all entries are non-positive expect possibly those in the diagonal) and A −1 ≥ 0 holds (several definitions exist for M-matrices, see [6,Theorem 5.1]).It is easy to see that T is a nonnegative matrix and Σ is a non-singular M-matrix.The proof of the next result follows by similar arguments as those in the proof of [26,Theorem 2], where we denote by ρ(A) the dominant eigenvalue of any square matrix A.
In the mathematical epidemiology literature, it is common to relate stability properties with the next-generation matrix (NGM) and the basic reproduction number (R 0 ) (for description of the next-generation method and computation of this matrix and number, see [4] and the references therein).If we define reproduction (the emergence of new infections) by the matrix T and transition between the infection classes by Σ then the NGM of system (3.6)arises as The dominant eigenvalue of the NGM is denoted by R 0 (that is, R 0 := ρ(K)), that is of particular importance as by Theorem 3.1 and Proposition 3.2 it is a threshold quantity for the stability of the DFE in the model (2.1).However, due to the complexity of the NGM it is not possible to derive the formula of R 0 and therefore we cannot explicitly obtain the stability condition.
The next section is devoted to addressing this issue as we consider two special cases for the movement network: first we investigate the case of two cities, then we obtain more general results as we assume that p cities are organized in a ring-like structure.

Two cities
For p = 2, the system (3.6)reduces to We define and obtain the coefficient matrix of the linear system as T 2 − Σ 2 .We derive that and observe that where the relation is meant entry-wise.We introduce the notation where R j and P j are the basic reproduction numbers in city j in the absence of movement and in the absence of inflow into city j, respectively.With these definitions, it holds that P 1 < R 1 and P 2 < R 2 .We arrive at the following result.
Theorem 4.1.If max { P 1 , P 2 } ≥ 1 then the zero solution of (4.1) is unstable for all choices of delays.
It remains to investigate stability in the case when P 1 < 1 and P 2 < 1.In this case, ρ(K 0 2 ) < 1 holds so it is meaningful to define that can be explicitly calculated as A threshold for stability in the system (4.1) is given by T 2 = 1.
Proof.We use Theorem 2.1 from [22,23] to establish the relationship between the dominant eigenvalue of

D. Knipl
We use T 2 to give stability conditions in terms of the reproduction numbers and movement rates.We assume that β T 1,2 − µ T 1,2 > 0 and β T 2,1 − µ T 2,1 > 0; such conditions are reasonable as evidence ( [5,27]) supports that the transmission rate of an infectious disease can be much higher than usual when a large number of passengers are sharing the same cabin during travel.Furthermore, we consider enables us to give critical values for the delays.We define and note that τ c > 0 if and only if the following condition holds: Theorem 4.3.Assume that P 1 < 1 and P 2 < 1 hold.If R 1 > 1 and R 2 > 1 then the zero solution of (4.1) is unstable for all choices of delays.If R 1 < 1 and R 2 < 1 then the zero solution is LAS for The assumption β T − µ T > 0 implies that the square root of the left hand side above is a lower bound for T 2 , therefore T 2 > 1 for all τ 1,2 and τ 2,1 , that is, the zero solution is unstable in the case when R 1 > 1 and R 2 > 1.
If R 1 < 1 and R 2 < 1 then by similar arguments as in the previous case, we derive that The last inequality is equivalent to T 2 < 1 for τ 1,2 = 0, τ 2,1 = 0, that is, the zero solution is stable in the model without delays.T 2 is monotonically increasing in the delays, and the condition T 2 < 1 for stability is equivalent to It can also be shown that , that is, τ c gives stability threshold for the delays.Note that In the two cases when R 1 < 1 and R 2 > 1, and when R 1 > 1 and R 2 < 1, the derivation of the stability condition T 2 < 1 ⇔ τ 1,2 + τ 2,1 < τ c goes analogously as above; however τ c might be non-positive in which case the zero solution is unstable for all choices of delays.A necessary and sufficient condition for τ c > 0 is that (c) holds, and in this case the zero solution is stable for τ 1,2 + τ 2,1 < τ c but it loses its stability when the sum of the delays exceeds τ c .

Ring of cities
For a movement network where all movement rates are zero except Using the short hand notations c 1 = e τ 1,p (β T 1,p −µ T 1,p ) and c j = e τ j,j−1 (β T j,j−1 −µ T j,j−1 ) for j = 2, . . ., p, we define and obtain the matrix of the linear system as T p − Σ p , where T p is a nonnegative matrix and Σ p is a non-singular M-matrix.Therefore it is meaningful to define , and similarly as in the two-city case, a relationship between s(T p − Σ p ) and ρ(K p ) can be established by Proposition 3.2, as follows: the zero solution of (4.
and observe that The proof of the following result is omitted due to its similarity to Theorem 4.1.
Theorem 4.4.If max 1≤j≤r P j ≥ 1 then the zero solution of (4.2) is unstable for all choices of delays.
If P j < 1 for all j = 1, 2, . . ., p then ρ(K 0 p ) < 1 clearly holds, so it is meaningful to define ), that can be explicitly calculated as By the assertion of Proposition 3.2 and [22, 23, Theorem 2.1] we obtain that the zero solution is LAS if and only if T p < 1, and the stability threshold T p = 1 can be used to derive explicit conditions for the delays.To this end, it is assumed that β where τ gives the sum of all delays, and note that τ p c > 0 is equivalent with the following condition: ( (c*) Theorem 4.5.If R j > 1 for all j then the zero solution of (4.2) is unstable for all choices of delays.If R j < 1 for all j then the zero solution is LAS for τ < τ p c whereas it is unstable for τ > τ p c .If there are j and k such that R j > 1 and R k < 1 then the zero solution is unstable for all choices of delays if the condition (c*) does not hold, it is also unstable if (c*) holds and τ > τ Proof.If R j > 1 for all j then it follows by that implies instability for all choices of delays by T p > 1.One can easily show that the above inequality is reversed if R j < 1 for all j, and it follows that τ p c > 0. Straightforward calculations yield that the zero solution is LAS in the model without delays because T p < 1 holds, and stability is preserved in the delayed system as long as τ < τ p c , where τ p c is the critical threshold for the total sum of delays, derived from T p = 1.
The condition (c*) is equivalent to T p < 1 in the special case when all delays are zero, moreover, T p is increasing in the delays hence the condition (c*) is necessary for stability in the case when there are j and k such that R j > 1 and R k < 1.The stability threshold τ p c for τ is calculated using T p = 1, and (c*) guarantees that τ p c > 0. The proof is complete.

Discussion on the impact of travel on stability
In this last section we return to system (4.1) to comment on the sensitivity of stability in the movement rates.The following computations are useful:

D. Knipl
We learnt from Theorem 4.3 that if R 1 < 1 and R 2 < 1 then for all choices of the movement rates there is a stability threshold τ c > 0. The above calculations reveal that the derivative of τ c is negative with respect to both movement parameters, which means that increasing the movement rates shrinks the stability region and therefore it can destabilize the zero solution.
Stability is also possible if the condition (c) holds and either R 1 < 1 and R 2 > 1 or R 1 > 1 and R 2 < 1.We investigate now the former case, and find that τ c increases in m 1,2 whereas it decreases in m 2,1 , thereby the two movement rates have opposite impact on τ c .The first conclusion that we can draw is that encouraging movement from the endemic city 2 to the non-endemic city 1 expands the stability region and helps stabilizing the DFE.On the other hand, increasing the movement rate in the other direction can lead to stability loss; there is a critical value for m 2,1 , given by such that τ c becomes negative when m 2,1 exceeds m u 2,1 (note that the condition (c) ceases to hold).An important observation to make is that m u 2,1 depends on the other movement rate m 1,2 .Lastly, we note that the conditions P 1 < 1 and P 2 < 1 in Theorem 4.3 give lower bounds for m 2,1 and m 1,2 , respectively, as Summarizing, in the case when R 1 < 1 and R 2 > 1 hold, stability is possible only if the movement rates are chosen such that m The above analysis revealed some non-trivial dependence of the stability region on the movement rates when R 1 < 1 and R 2 > 1.To further investigate this question, we let m 1,2 = m 2,1 =: m and study the dependence of τ c on the unified movement rate m.In this special case the function H arises as We assume that R 1 < 1, R 2 > 1 and the condition (c) hold, or equivalently µ and m is chosen such that H(m) > 1.To satisfy the last inequality, it follows by arises for m, that is weaker than the one derived above from H(m) > 1. Straightforward calculations yield and we obtain that where we used that R 1 and R 2 are chosen such that (µ It is easy to see that m u > m H > 0, therefore we arrive at the following conclusions.If R 1 and R 2 are chosen such that (µ 1 + d 1 − β 1 + µ 2 + d 2 − β 2 ) > 0 does not hold then the zero solution is unstable for all movement rates and delays.If the last inequality is reversed, then with words, for all movement rates m such that m > m H stability is possible for small delays, however increasing m expands and then shrinks the stability region.

Appendix B Calculation of the gradients in the system (3.2)
To derive the gradients on the right hand side of the system (3.2),we recall the definition I j,k (v), S j,k (v), R j,k (v) = ỹj,k (τ j,k ; 0, h j,k (v)) from (2.4), and it follows that Dg j,k ( ỹj,k (ν;0,h j,k (v))) dν Dh j,k (v).
Proof.Theorem 3.3 in Chapter I in [8] states that as g j,k has continuous first derivative, the solution ỹj,k (θ; 0, y * ) of system (2.3) is continuously differentiable with respect to θ and y * on its domain of definition.The matrix Dg j,k ( ỹj,k (ν;0,h j,k (v))) dν Dh j,k (v).
If h j,k (v) is an equilibrium of the system (2.3) then ỹj,k (ν; 0, h j,k (v)) = h j,k (v) holds for all ν ∈ [0, τ j,k ], and therefore the above formula reduces to ∂ ∂v ỹj,k (τ j,k ; 0, h j,k (v)) = e τ j,k Dg j,k (h j,k (v)) Dh j,k (v).where z : R → R n , L : R → L(C r , R n ) is continuous.For the segment of the solution we use the usual notation z t where z t = (z t 1 , . . ., z t n ), z t i = z i (t + θ) for θ ∈ [−r i , 0].Let L i (t)φ denote the i-th component of L(t)φ.We introduce the following definition for the property (K).
Lemma C.2. (K) holds if and only if there exists a i (t) ∈ R for 1 ≤ i ≤ n and positive Borel measures η ij (t) for 1 ≤ i, j ≤ n on [−r j , 0] such that L(t) has a representation in terms of Borel measures, as and η ii (t){0} = 0.Moreover, if (K) holds then the representation is unique and a i (t) and η ij (t) are continuous functions of t.
For 1 ≤ i ≤ n let e i denote the standard basis vectors for R n , and let êi ∈ C r be the constant function equal to e i for all values of its argument.We need two further conditions.Definition C.3.(R) For each j for which r j > 0, there exists i such that for all t, η ij (t)([−r j , −r j + )) > 0 for all small > 0.
If (K) holds then by Lemma C.2 we obtain that A(L) ij (t) = L i (t) êj = a i (t)(e j ) i + η ij (t)([−r j , 0]) = a i (t) + η ii (t)([−r i , 0]) if i = j, The condition (K) implies that if a component of a solution of the differential equation gets turned on, then it stays on.We need the condition (R) to ensure that some component actually gets turned on if the initial condition is nontrivial.(I) is a kind of irreducibility assumption; it implies that once one component gets turned on then all components eventually get turned on.
A general autonomous functional differential equation (FDE) is denoted by where f : U → R n is continuously differentiable and U is an open subset of C r .For ψ ∈ U, D f (ψ) gives a bounded linear operator, that is, D f (ψ)φ ∈ R n for φ ∈ C r , and the corresponding linear system arises as d dt z(t) = D f (ψ)z t .

Theorem 3 . 1 .
The DFE of the model (2.1) is LAS if and only if the zero solution of system (3.6) is LAS.Proof.By the above arguments, the stability analysis of the DFE in the model (2.1) reduces to studying stability in the systems (3.4) and(3.5).It follows from[20, Theorem 4.2]  that the zero solution is LAS in system (3.5)(see system(4.4)in[20] and[9, Theorem 1]).
it is LAS if (c*) holds and τ < τ p c .