As a submodel, an age-structured system is constructed to incorporate the possibility of disease transmission during travel, where age is the time elapsed since the start of the travel

Vol. 0, No. 0, pp. 000–000

Epidemic Spread and Variation of Peak Times in Connected Regions due to Travel-Related Infections—Dynamics of an Antigravity-Type Delay Differential

Model^∗

Di´ana H. Knipl^†, Gergely R¨ost^‡, and Jianhong Wu^§

Abstract. National boundaries have never prevented infectious diseases from reaching distant territories; how- ever, the speed at which an infectious agent can spread around the world via the global airline transportation network has significantly increased during recent decades. We introduce an SEAIR- based, antigravity model to investigate the spread of an infectious disease in two regions which are connected by transportation. As a submodel, an age-structured system is constructed to incorporate the possibility of disease transmission during travel, where age is the time elapsed since the start of the travel. The model is equivalent to a large system of differential equations with dynamically defined delayed feedback. After describing fundamental but biologically relevant properties of the system, we detail the calculation of the basic reproduction number and obtain disease transmission dynamics results in terms ofR0. We parametrize our model for influenza and use real demographic and air travel data for the numerical simulations. To understand the role of the different charac- teristics of the regions in the propagation of the disease, three distinct origin-destination pairs are considered. The model is also fitted to the first wave of the influenza A(H1N1) 2009 pandemic in Mexico and Canada. Our results highlight the importance of including travel time and disease dynamics during travel in the model: the invasion of disease-free regions is highly expedited by elevated transmission potential during transportation.

Key words. differential equations, transportation model, epidemic spread, influenza modeling AMS subject classifications. Primary, 34K05; Secondary, 92D30

DOI. 10.1137/130914127

1. Introduction. The global network of human transportation has played a paramount role in the spatial spread of infectious diseases. The high connectedness of distant territories by air travel makes it possible for a disease to invade regions far away from the source faster than ever. Some infectious diseases, such as tuberculosis, measles, and seasonal influenza, have been known to be transmissible during commercial flights. The importance of the global

∗Received by the editors March 22, 2013; accepted for publication (in revised form) by J. Sneyd July 17, 2013;

published electronically DATE. The work of the first and second authors was partly supported by the European Union and the European Social Fund through project FuturICT.hu (grant T´AMOP–4.2.2.C-11/1/KONV–2012–0013). The project was also subsidized by the European Union and co-financed by the European Social Fund.

http://www.siam.org/journals/siads/x-x/91412.html

†MTA–SZTE Analysis and Stochastics Research Group, Hungarian Academy of Sciences, Hungary (knipl@math.

u-szeged.hu). This author’s work was partially supported by the T´AMOP 4.2.4. A/2–11–1–2012–0001 “National Excellence Program — Elaborating and operating an inland student and researcher personal support system conver- gence program.”

‡Bolyai Institute, University of Szeged, Hungary (rost@math.u-szeged.hu). This author’s work was partially supported by the European Research Council StG Nr. 259559.

§Centre for Disease Modelling, Laboratory for Industrial and Applied Mathematics, York University, Toronto, ON, Canada (wujh@mathstat.yorku.ca). This author’s work was partially supported by Natural Sciences and Engineering Research Council of Canada and by the Canada Research Chair Program.

1

air travel network was highlighted in the 2002–2003 SARS outbreak (WHO [38]) and clearly contributed to the global spread of the 2009 pandemic influenza A(H1N1) (Khan et al. [15]).

Therefore, mathematically describing the spread of infectious diseases on the global human transportation network is of critical public health importance.

There are a few well-known studies which constructed and analyzed various meta-population models for disease spread in connecting regions (see Arino [2], Arino and van den Driessche [3], Ruan, Wang, and Levin [25], Wang and Zhao [37], and the references therein). These studies focus mainly on the impact of spatial dispersal of infected individuals from one region to another, and do not consider transportation as a platform of disease dynamics. However, during long-distance travel, such as intercontinental flights, a single infected individual may infect several other passengers (Wagner, Coburn, and Blower [36], European Center for Dis- ease Prevention an Control [12]), thus potentially inducing multiple generating infections in the destination region. It is therefore desirable to properly describe the spread of the disease via long-distance travel, in a way that incorporates into the models the transmission dynamics during the transition and travel.

Cui, Takeuchi, and Saito [8] and Takeuchi, Liu, and Cui [31] modeled the possibility that individuals may contract a disease while traveling by a system of ordinary differential equations based on the standard SIS (susceptible-infected-susceptible) epidemic model. They discovered that the disease can persist in regions connected by human transportation even if the infection died out in all regions in the absence of travel. Liu, Wu, and Zhou [18]

noted that the previously proposed models [8, 19, 31] implicitly used the assumption that the transportation between regions occur instantaneously. For some diseases of major public health concern, such as SARS and influenza, the progress of the disease is so fast that even a short delay (a fraction of a day) can be significant. Based on such considerations, Liu, Wu, and Zhou [18] introduced the time needed to complete the travel into the SIS-type epidemic model, and also the possible infections during this time. Nakata [22] described the global dynamics of this system for two identical regions in terms of the basic reproduction number.

The model was later generalized by Nakata and R¨ost [23] to the case ofnregions with different characteristics and arbitrary travel networks. An SIR (susceptible-infected-recovered)-based model with a general incidence term was analyzed in Knipl and R¨ost [17] to describe the spread of infection in multiple regions with travel considered.

The purpose of this work is to formulate a model to properly describe the temporal evo- lution of an epidemic in regions connected by long-distance travel, such as intercontinental flights. The European Centre for Disease Prevention and Control (ECDC) developed a risk as- sessment guideline [12] for infectious diseases transmitted on aircrafts, like influenza. Existing studies confirmed that on-board transmission was possible in flights even with a duration of less than eight hours. For most diseases that pose a threat of a global pandemic, an SIS-type model is not adequate. For this reason, here we use the SEAIR model as a basic epidemic model building block in the regions and also during travel. The SIS model can be reduced to a logistic equation and then can be solved analytically. This property was heavily used in the analysis done in [18,22,23]. However, the lack of closed form solution causes substantial technical difficulties in the analysis of SEAIR-type models, as will be shown in this paper.

More significantly the aforementioned existing models did not distinguish local residents from temporary visitors in the model setup. In reality, the large part of travels are return

trips, and not only the number of visitors, but also the average time that visitors spend in the other region may significantly affect the speed of spatial spread of the disease. If visitors spend more time in a region which is a hotspot of the disease, they will more likely carry the disease back to their region of origin. In addition, visitors and local residents may have very different contact rates and mixing patterns, for example if the visitors are typically on holiday and stay in selected resorts and hotels. Hence in our model we use different compartments for residents and visitors to capture this phenomenon.

Many multiregional epidemic models, specially the gravity-type models, are based on the assumption that the speed of the spread of epidemics between regions is inversely proportional to the distance between those regions (see, for example, Tuite et al. [32] for the recent cholera outbreak in Haiti). However, in case of air travel, the travel behavior is different and can be just the opposite. First, the number of travelers does not depend directly on the distance between regions, but is determined by other more important factors, such as business and cultural relations or touristic attractivity. Second, the transmission rate of an infectious disease can be much higher than usual when a large number of passengers are sharing the same cabin, and the longer the flight (which means the larger distance is between regions), the greater the number of infections that can be expected (Wagner, Coburn, and Blower [36]).

Hence the air travel model we are proposing here is in principle “antigravity.”

The paper is organized as follows. In the next section we formulate an age-structured model (where age means time since the start of travel), which leads to a nonlinear system of functional differential equations. In section 3 we determine some fundamental properties of the model. Section 4 is concerned with the computation of local and global reproduction numbers. We parametrize our model for influenza in section 5, and then in section 6 we introduce three prototype origin-destination pairs (Canada-Mexico, Canada-China, Canada- UK) and run simulations using real air traffic and tourism data. In the last section we discuss our findings.

2. Model description. We formulate a dynamical model describing the spread of an in- fectious disease within and between two regions, and also during travel from one region to the other. We divide the entire populations of the two regions into the disjoint classes S_j^m, E_j^m, A^m_j , I_j^m, R^m_j , j ∈ {1,2}, m ∈ {r, v}, where the letters S, E, A, I, and R represent the compartments of susceptible, exposed, asymptomatic infected, symptomatic infected, and recovered individuals, respectively. Lower index j∈ {1,2} specifies the current region, upper indexm∈ {r, v}denotes the residential status of the individual in the current region (resident versus visitor). For instance, S₁^v is the compartment of individuals who are susceptible to the disease and staying in region 1 as a visitor (hence, they originally belong to region 2); members of A^r₂ are those who are asymptomatic infected residents in region 2.

LetS_j^m(t),E_j^m(t),A^m_j (t),I_j^m(t),R^m_j (t),j∈ {1,2},m∈ {r, v}, be the number of individu- als belonging toS_j^m,E^m_j ,A^m_j ,I_j^m,R^m_j , respectively, at timet. The transmission rate between an infected individual with residential status m and a susceptible individual with residential status nin region j (j∈ {1,2}, m, n∈ {r, v}) is denoted byβ_j^m,n. Let F_j^r denote the force of infection of residents, and F_j^v the force of infection of visitors in region j. Model parameter μ_E denotes the inverse of the incubation period, μ_A and μ_I are the recovery rates of asymp- tomatic and symptomatic infected individuals. Let ρbe the reduction factor of infectiousness

of asymptomatic infected individuals (we assume they are capable of transmitting the disease, but generally with a lower rate than symptomatic infected individuals). Let p denote the probability that an infected individual develops symptoms, and let δ denote disease-induced mortality rate. We assume constant recruitment terms Λ_j, while d^r_j and d^v_j denote natural mortality rates of residents and visitors in region j. We denote the travel rate of residents between region j and region k by α_j, and the rate at which visitors of region j travel back to region k by γ_j; thus 1/γ_j is the average time visitors spend in region j. For the total population of residents, visitors, and all individuals currently being in region j at time t, we use the notation

(1)

N_j^r(t) =S_j^r(t) +E_j^r(t) +A^r_j(t) +I_j^r(t) +R_j^r(t), N_j^v(t) =S_j^v(t) +E_j^v(t) +A^v_j(t) +I_j^v(t) +R_j^v(t),

N_j(t) =N_j^r(t) +N_j^v(t).

We divide the population during travel into the classes s^m_j,k,e^m_j,k,a^m_j,k,i^m_j,k,r_j,k^m. Letterss, e,a, i, r denote susceptible, exposed, asymptomatic infected, symptomatically infected, and recovered travelers, respectively. Lower indices j, k ∈ {1,2}, j =k, indicate that individuals are traveling from region j to region k. Upper index m ∈ {r, v} determines individuals’

residential status in the region they have just left: for instance, an individual who is now in r_1,2^v is recovered, traveling from region 1 to region 2, and was a visitor in region 1, which means that the individual originally belongs to region 2.

Let τ >0 denote the average time required to complete a one-way trip. To describe the disease dynamics during travel, we defines^m_j,k(θ, t_∗),e^m_j,k(θ, t_∗),a^m_j,k(θ, t_∗),i^m_j,k(θ, t_∗),r^m_j,k(θ, t_∗), j, k ∈ {1,2}, j =k,m∈ {r, v}, as the density of individuals who started travel at timet^∗ and belong to classes s^m_j,k, e^m_j,k, a^m_j,k,i^m_j,k, r_j,k^m with respect to θ, whereθ ∈[0, τ] denotes the time elapsed since the beginning of the travel. Let

(2) n^m_j,k(θ, t_∗) =s^m_j,k(θ, t_∗) +e^m_j,k(θ, t_∗) +a^m_j,k(θ, t_∗) +i^m_j,k(θ, t_∗) +r^m_j,k(θ, t_∗), where j, k∈ {1,2}, j =k, m∈ {r, v}, and let

(3) n_j,k(θ, t_∗) =n^r_j,k(θ, t_∗) +n^v_j,k(θ, t_∗).

Thus, _θ1

θ2 n_j,k(θ, t−θ)dθ is the number of individuals who left region j in the time interval [t−θ₁, t−θ₂], where τ ≥ θ₁ ≥θ₂ ≥ 0. In particular, for θ₁ = τ and θ₂ = 0, this gives the total number of individuals who are in the travel transition from region j to region kat time t. We assume that infected individuals do not die during travel; hencen_j,k(θ, t_∗) =n_j,k(0, t_∗) for all θ∈[0, τ]. During the course of travel, infected individuals can transmit the disease at the rate β^T. We use the notation μ^T_E, μ^T_A, μ^T_I for the inverse of the incubation period and the recovery rates of asymptomatic and symptomatic infected individuals during travel. Let F_j,k^T denote the force of infection during travel from region j to region k. Thens^m_j,k(τ, t−τ), e^m_j,k(τ, t−τ), a^m_j,k(τ, t−τ), i^m_j,k(τ, t−τ), r^m_j,k(τ, t−τ) gives the inflow of individuals arriving from region jto compartments S_kⁿ,E_kⁿ,Aⁿ_k,I_kⁿ,Rⁿ_k,j, k ∈ {1,2},j=k,m, n∈ {r, v}, m=n, respectively, at time t.

All variables and model parameters are listed in Tables 1 and 2. The flow chart of the model is depicted in Figure 1. Based on the assumptions formulated above, we obtain the following system of differential equations for disease transmission in the two regions:

Table 1

Model variables (j,k ∈ {1,2},j=k). In the figure, “density” means the density with respect to the age since the start of travel.

Variables

F_j^r Force of infection of residents in regionj Fj^v Force of infection of visitors in regionj

F_j,k^T Force of infection during travel from region jto regionk Sj^r,E^rj,A^rj,Ij^r,R^rj Susceptible, exposed, asymptomatic, symptomatic

infected, or recovered residents in regionj Sj^v,Ej^v,A^vj,Ij^v,R^vj Susceptible, exposed, asymptomatic, symptomatic

or infected, recovered visitors in regionj Nj^r,Nj^v,Nj Total population size of residents, visitors, and

all individuals in region j

s^r_j,k,e^r_j,k,a^r_j,k,i^r_j,k,r^r_j,k Density of susceptible, exposed, asymptomatic, symptomatic infected, or recovered individuals during the travel fromj tok (traveling to visitk) s^vj,k,e^vj,k,a^vj,k,i^vj,k,r^vj,k Density of susceptible, exposed, asymptomatic,

symptomatic infected, or recovered individuals during the travel fromjtok(returning tokfrom visitingj) n^r_j,k,n^v_j,k,nj,k Total density of residents, visitors, and

all individuals during the travel fromjtok

Table 2

Key model parameters (j,k ∈ {1,2},j=k).

Key model parameters Λ_j Recruitment rate in regionj

d^r_j,d^v_j Natural death rate of residents and visitors of regionj δ disease-induced death rate

β^m,n_j Transmission rate between an infected individual with residential statusmand a susceptible individual with residential statusnin regionj(m, n∈ {r, v}) β^T Transmission rate during the travel

αj Traveling rate of residents of regionjto regionk γj Inverse of duration of visitors’ stay in regionj

τ Duration of travel between the regions p Probability of developing symptoms

ρ Reduction of infectiousness of asymptotic infecteds μ_E,μ^T_E Reciprocal of the length of the incubation period

in the regions and during the travel μA,μ^T_A Recovery rate of asymptomatic infecteds

in the regions and during the travel μI,μ^T_I Recovery rate of symptomatic infecteds

in the regions and during the travel

E

A

I

S

R

S

E

I

A

R

S

E

A

I

R

S

E

A

I

R

Travel

e

r

s

s

a

i

i

a

e

r

s

e

i

a

r

e

a

i

r

s

Region 1 Region 2

Figure 1. Color-coded flow chart of disease transmission and travel dynamics. The disease transmission in the two regions is shown in two different columns; the disease progresses vertically from top to bottom. Classes having the same origins are marked by the same colors. Red corresponds to the classes originating from region 1; blue represents classes of region2. Arrows colored with the same colors indicate how the disease progresses.

Green dashed-dotted arrows represent individuals that are traveling. Green solid arrows show the dynamics of the pandemic during the course of the travel. The description of the variables can be found in Table1.

(L)

⎧⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎨

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎩

S˙_j^r(t) = Λ_j−S_j^r(t)F_j^r(t)−(d^r_j +α_j)S_j^r(t) +s^v_k,j(τ, t−τ), E˙_j^r(t) =S_j^r(t)F_j^r(t)−(d^r_j +μ_E +α_j)E_j^r(t) +e^v_k,j(τ, t−τ), A˙^r_j(t) = (1−p)μ_EE_j−(d^r_j +α_j+μ_A)A^r_j(t) +a^v_k,j(τ, t−τ),

I˙_j^r(t) =pμ_EE_j−(d^r_j +α_j+δ+μ_I)I_j^r(t) +i^v_k,j(τ, t−τ), R˙^r_j(t) =μ_II_j^r(t) +μ_AA^r_j(t)−(d^r_j +α_j)R^r_j(t) +r_k,j^v (τ, t−τ), S˙_j^v(t) =−S_j^v(t)F_j^v(t)−(d^v_j +γ_j)S_j^v(t) +s^r_k,j(τ, t−τ), E˙_j^v(t) =S_j^v(t)F_j^v(t)−(d^v_j +μ_E+γ_j)E_j^v(t) +e^r_k,j(τ, t−τ), A˙^v_j(t) = (1−p)μ_EE_j^v(t)−(d^v_j +γ_j +μ_A)A^v_j(t) +a^r_k,j(τ, t−τ),

I˙_j^v(t) =pμ_EE_j^v(t)−(d^v_j +γ_j +δ+μ_I)I_j^v(t) +i^r_k,j(τ, t−τ), R˙^v_j(t) =μ_II_j^v(t) +μ_AA^v_j(t)−(d^v_j +γ_j)R^v_j(t) +r_k,j^r (τ, t−τ), where

F_j^r(t) = 1 N_j(t)

β_j^rr(I_j^r(t) +ρA^r_j(t)) +β_j^vr(I_j^v(t) +ρA^v_j(t)) , F_j^v(t) = 1

N_j(t)

β_j^rv(I_j^r(t) +ρA^r_j(t)) +β_j^vv(I_j^v(t) +ρA^v_j(t)) .

For each given t_∗, the following system (T) describes the evolution of the densities during the travel initiated at time t_∗:

(T)

⎧⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎨

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎩ d

dθs^r_j,k(θ, t_∗) =−s^r_j,k(θ, t_∗)F_j,k^T (θ, t_∗), d

dθe^r_j,k(θ, t_∗) =s^r_j,k(θ, t_∗)F_j,k^T (θ, t_∗)−μ^T_Ee^r_j,k(θ, t_∗), d

dθa^r_j,k(θ, t_∗) = (1−p)μ^T_Ee^r_j,k(θ, t_∗)−μ^T_Aa^r_j,k(θ, t_∗), d

dθi^r_j,k(θ, t_∗) =pμ^T_Ee^r_j,k(θ, t_∗)−μ^T_Ii^r_j,k(θ, t_∗), d

dθr_j,k^r (θ, t_∗) =μ^T_Aa^r_j,k(θ, t_∗) +μ^T_Ii^r_j,k(θ, t_∗), d

dθs^v_j,k(θ, t_∗) =−s^v_j,k(θ, t_∗)F_j,k^T (θ, t_∗), d

dθe^v_j,k(θ, t_∗) =s^v_j,k(θ, t_∗)F_j,k^T (θ, t_∗)−μ^T_Ee^v_j,k(θ, t_∗), d

dθa^v_j,k(θ, t_∗) = (1−p)μ^T_Ee^v_j,k(θ, t_∗)−μ^T_Aa^v_j,k(θ, t_∗), d

dθi^v_j,k(θ, t_∗) =pμ^T_Ee^v_j,k(θ, t_∗)−μ^T_Ii^v_j,k(θ, t_∗), d

dθr_j,k^v (θ, t_∗) =μ^T_Aa^v_j,k(θ, t_∗) +μ^T_Ii^v_j,k(θ, t_∗),

where j, k∈ {1,2}, j =k, and F_j,k^T (θ, t_∗) = β^T

n_j,k(θ, t_∗)(i^r_j,k(θ, t_∗) +i^v_j,k(θ, t_∗) +ρ(a^r_j,k(θ, t_∗) +a^v_j,k(θ, t_∗))), n_j,k(θ, t_∗) =α_j(S_j^r(t_∗) +E_j^r(t_∗) +A^r_j(t_∗) +I_j^r(t_∗) +R^r_j(t_∗))

+ γ_j(S_j^v(t_∗) +E_j^v(t_∗) +A^v_j(t_∗) +I_j^v(t_∗) +R^v_j(t_∗))

=α_jN_j^r(t_∗) +γ_jN_j^v(t_∗).

For θ = 0, the densities are determined by the rates at which individuals start their travels from one region to the other at time t_∗. Hence, the initial values for system (T) atθ= 0 are given by

(IVT)

⎧⎪

⎪⎪

⎪⎪

⎪⎪

⎨

⎪⎪

⎪⎪

⎪⎪

⎪⎩

s^r_j,k(0, t_∗) =α_jS_j^r(t_∗), s^v_j,k(0, t_∗) =γ_jS₁^v(t_∗), e^r_j,k(0, t_∗) =α_jE_j^r(t_∗), e^v_j,k(0, t_∗) =γ_jE_j^v(t_∗), a^r_j,k(0, t_∗) =α_jA^r_j(t_∗), a^v_j,k(0, t_∗) =γ_jA^v_j(t_∗), i^r_j,k(0, t_∗) =α_jI_j^r(t_∗), i^v_j,k(0, t_∗) =γ_jI_j^v(t_∗), r^r_j,k(0, t_∗) =α_jR_j^r(t_∗), r^v_j,k(0, t_∗) =γ_jR_j^v(t_∗) forj, k ∈ {1,2}, j=k.

Now we turn our attention to the termss^m_j,k(τ, t−τ),e^m_j,k(τ, t−τ),a^m_j,k(τ, t−τ),i^m_j,k(τ, t−τ), r^m_j,k(τ, t−τ) in system (L), which are the densities of individuals arriving to classes S_kⁿ, E_kⁿ, Aⁿ_k,I_kⁿ,Rⁿ_k,j, k∈ {1,2},j =k,m, n∈ {r, v}, m=n, respectively, at timet, upon completing a one-way trip from regionj. At timet, these terms are determined by the solution of system (T) with initial values (IVT) for t_∗ =t−τ at θ=τ:

(i) individuals who enter region kat timet are those who left regionj at timet−τ; (ii) residents of region j become visitors of region k and vice versa (m=n) upon com- pleting a one-way trip;

(iii) an individual may move to a different compartment during travel; for example, a susceptible resident who travels from region j may arrive as an infected visitor to region k (j, k∈ {1,2}, j=k), as given by the dynamics of system (T).

Next we specify initial values for system (L) att= 0. Since travel takesτ units of time to complete, arrivals to region j are determined by the state of region k(j, k∈ {1,2}, j=k) at t−τ, via the solution of systems (T) and (IVT). Thus, we set up initial functions as follows:

(IVL)

⎧⎪

⎪⎪

⎪⎪

⎪⎪

⎨

⎪⎪

⎪⎪

⎪⎪

⎪⎩

S_j^r(u) =ϕ^r_S,j(u), S_j^v(u) =ϕ^v_S,j(u), E_j^r(u) =ϕ^r_E,j(u), E_j^v(u) =ϕ^v_E,j(u), A^r_j(u) =ϕ^r_A,j(u), A^v_j(u) =ϕ^v_A,j(u), I_j^r(u) =ϕ^r_I,j(u), I_j^v(u) =ϕ^v_I,j(u), R^r_j(u) =ϕ^r_R,j(u), R^v_j(u) =ϕ^v_R,j(u),

where u ∈ [−τ,0], and each ϕ^m_K,j is a continuous function for j ∈ {1,2}, m ∈ {v, r}, K ∈ {S, E, A, I, R}.

Note that systems (L) and (T) are interconnected; in order to determine the dynamics of the model, simultaneous solution of both is required. Considering the fact that disease transmission is possible during travel, the solution of system (T) at (τ, t−τ) is required for allt≥0 to find the solution of (L). However, in order to obtain the solution of (T) at (τ, t−τ), it is necessary to use the solution of (L) att−τ, because (T) takes the initial conditions from (L). Hence, in order to describe the disease transmission in the regions, the solution of another differential equation system is required at each time t, which has initial values depending on the earlier state of the system on the regions. Thus (L) is a delay differential system, where the delayed feedback is determined by a solution of a parallel system of ordinary differential equations. In previous papers with travel delay, such as [18, 22, 23], the authors used an SIS-type system during travel, which was analytically solvable; thus it was possible to express the delayed feedback explicitly. Unlike the SIS model, the SEAIR model is not analytically solvable; therefore here we have to deal with a system of functional differential equations, where the delay term is given only implicitly via a solution of a nonlinear system of ordinary differential equations.

3. Basic properties of the model. In this section, we show that our model is equivalent to a system of nonlinear functional differential equations where the delay term is defined dynam- ically, via the solution of another system of differential equations. Then we also investigate some biologically relevant properties of the system. Set

X_j^r=

⎛

⎜⎜

⎜⎜

⎜⎝ S_j^r E_j^r A^r_j I_j^r R^r_j

⎞

⎟⎟

⎟⎟

⎟⎠

, X_j^v =

⎛

⎜⎜

⎜⎜

⎜⎝ S_j^v E_j^v A^v_j I_j^v R^v_j

⎞

⎟⎟

⎟⎟

⎟⎠

, x^r_j,k =

⎛

⎜⎜

⎜⎜

⎜⎝ s^r_j,k e^r_j,k a^r_j,k i^r_j,k r^r_j,k

⎞

⎟⎟

⎟⎟

⎟⎠

, x^v_j,k =

⎛

⎜⎜

⎜⎜

⎜⎝ s^v_j,k e^v_j,k a^v_j,k i^v_j,k r_j,k^v

⎞

⎟⎟

⎟⎟

⎟⎠ ,

where j, k∈ {1,2}, j =k, and set

X=

⎛

⎜⎜

⎝ X₁^r X₁^v X₂^r X₂^v

⎞

⎟⎟

⎠, x=

⎛

⎜⎜

⎝ x^v_2,1 x^r_2,1 x^v_1,2 x^r_1,2

⎞

⎟⎟

⎠,

so that X∈R²⁰and x∈R²⁰.For a given t_∗ we define the system (T^∗) as (T^∗)

∂

∂θx(θ, t_∗) =f(x(θ, t_∗)), x(0, t_∗) =g(X(t_∗)), where t_∗, θ∈R+, f, g:R²⁰ R²⁰,

g_i(y) =

⎧⎪

⎪⎪

⎪⎨

⎪⎪

⎪⎪

⎩

γ₂y_i if i= 1, . . . ,5, α₂y_i if i= 6, . . . ,10, γ₁y_i if i= 11, . . . ,15, α₁y_i if i= 16, . . . ,20,

and f_i(x) equals the right-hand side of the equation for x_i in system (T). For instance, f₇(x) =β^T x₆

₁₀

j=1x_j(x₄+x₉+ρ(x₃+x₈))−μ^T_Ex₇.

Define G:R²⁰ R²⁰, whereG_i(X) is given by the right-hand side of the equation of X_i in (L) without the inflow from travel. For instance,

G₁₆(X) =− X_{16 20}

j=11X_j(β₂^rv(X₁₄+ρX₁₃) +β₂^vv(X₁₉+ρX₁₈))−(d^v₂ +γ₂)X₁₆.

Let ˆx(θ, t_∗;Y) denote the solution of the initial value problem (T^∗) for t_∗ with initial value ˆ

x(0, t_∗) = g(Y), where Y ∈ R²⁰, and let H(Y) := ˆx(τ, t−τ;Y), H : R²⁰ R²⁰. Then our system (L) can be written in closed form as a system of functional differential equations, (L^∗) X(t) =˙ F(X(t), X(t−τ)),

where F(X(t), X(t−τ)) = G(X(t)) + H(X(t−τ)), F : R²⁰×R²⁰ R²⁰. Clearly (T^∗) is also a compact form of (T). To study the dynamics of (L^∗), we define our phase space as the nonnegative cone C+ =C([−τ,0],R²⁰₊) of the Banach space of continuous functions from [−τ,0] to R²⁰, equipped with the supremum norm. For each Φ ∈ C₊, standard arguments guarantee that there exists a unique solution of system (L^∗) with initial valuesX(u) = Φ(u), u∈[−τ,0] (see [16,17]). Using the notation of (IVL), we have Φ = (Φ^r₁,Φ^v₁,Φ^r₂,Φ^v₂)^T, where Φ^r_j = (ϕ^r_S,j, ϕ^r_E,j, ϕ^r_A,j, ϕ^r_I,j, ϕ^r_R,j)^T, Φ^v_j = (ϕ^v_S,j, ϕ^v_E,j, ϕ^v_A,j, ϕ^v_I,j, ϕ^v_R,j)^T,j∈ {1,2}.

Proposition 1. For any Φ∈ C₊, the solution of system (L^∗) is nonnegative.

Proof. It is straightforward to see that (T^∗) preserves nonnegativity; thus H(Y) ≥ 0 if Y ≥ 0. Since our function F(y, z) : R²⁰₊ ×R²⁰₊ R²⁰ and Fy(y, z) are continuous on R²⁰ and for every i= 1, . . .20, for every y, z ∈ R²⁰₊, y_i = 0 implies Fi(y, z) ≥0, all the conditions of Theorem 3.4 in [27] hold. This implies that for nonnegative initial data the corresponding solution of system (L^∗) remains nonnegative.

We define the disease-free subspace C₊^df as C₊^df =

Φ|Φ = (ϕ^r_S,j,ˆ0,ˆ0,ˆ0, ϕ^r_R,j, ϕ^v_S,j,ˆ0,ˆ0,ˆ0, ϕ^v_R,j)^T

⊂ C₊, where ˆ0 denotes the constant 0 function. If Φ∈ C^df₊, then

E_j^r(t) =E_j^v(t) =A^r_j(t) =A^v_j(t) =I_j^r(t) =I_j^v(t)≡0 for all t≥0, and hence the disease-free subspace is positively invariant.

Proposition 2. In the disease-free subspace C₊^df there exists a unique positive equilibrium of system (L^∗) which is globally asymptotically stable inC₊^df.

Proof. Using the definition ofN_j^r andN_j^v (see Table1) in section2, for these variables we derive the following differential equation system:

(4)

N˙_j^r(t) = Λ_j−(d^r_j +α_j)N_j^r(t) +γ_kN_k^v(t−τ), N˙_j^v(t) =−(d^v_j +γ_j)N_j^v(t) +α_kN_k^r(t−τ),