Electronic Journal of Differential Equations, Vol. 2020 (2020), No. 62, pp. 1–29. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu STABILITY AND SENSITIVITY ANALYSIS OF THE EPIDEMIOLOGICAL MODEL BE-CODIS PREDICTING THE SPREAD OF HUMAN DISEASES BETWEEN COUNTRIES BENJAMIN IVORRA, DIÈNE NGOM, ANGEL M. RAMOS Communicated by Jesus Ildefonso Diaz Abstract. The Ebola virus disease is a lethal human and primate disease that requires a particular attention from the international health authorities due to important recent outbreaks in some Western African countries and isolated cases in Europe and North-America. Regarding the emergency of this situation, various decision tools, such as mathematical models, were de- veloped to assist the authorities to focus their efforts in important factors to eradicate Ebola. In a previous work, we proposed an original deterministic spatial-temporal model, called Be-CoDiS (Between-Countries Disease Spread), to study the evolution of human diseases within and between countries by tak- ing into consideration the movement of people between geographical areas. This model was validated by considering numerical experiments regarding the 2014-16 West African Ebola Virus Disease epidemic. In this article, we perform a stability analysis of Be-CoDiS. Our first objective is to study the equilibrium states of simplified versions of this model, limited to the cases of one or two countries, and determine their basic reproduction ratios. Then, we perform a sensitivity analysis of those basic reproduction ratios regarding the model pa- rameters. Finally, we validate the results by considering numerical experiments based on data from the 2014-16 West African Ebola Virus Disease epidemic. 1. Introduction Modelling and simulation are important decision tools that can be used to control human and animal diseases [1, 22, 31, 35]. However, since each disease exhibits its own biological characteristics, the models need to be adapted to each specific case to be able to handle real situations [6, 45]. In a previous work [23], we presented a spatial-temporal epidemiological model, called Be-CoDiS (Between-Countries Disease Spread), for the study of the spread of human diseases between and within countries. This model is an adaptation of a previous one, called Be-FAST (Between Farm Animal Spatial Transmission), which simulates the spread of animal diseases between and within farms [22, 31, 29, 30, 28]. Be-CoDiS is based on the combination of a deterministic Individual-Based model (where countries are considered as individuals) [11], simulating the between-country 2010 Mathematics Subject Classification. 92D30, 92C60, 34D20, 34K28. Key words and phrases. Epidemiological modelling; deterministic models; stability analysis; sensitivity analysis; ebola virus disease; basic reproduction ratio. c©2020 Texas State University. Submitted April 24, 2020. Published June 18, 2020. 1 2 B. IVORRA, D. NGOM, A. M. RAMOS EJDE-2020/62 interactions (here, migratory flux, also considered in other works such as [8, 25]) and disease spread, with a deterministic compartmental model [6, 17] (a system of ordi- nary differential equations), simulating the within-country disease spread. At the end of a simulation, Be-CoDiS returns outputs referring to outbreaks characteristics (for instance, the epidemic magnitude, the risk of disease introduction or spread per country, etc.). The main characteristic of this approach is the consideration of the following effects at the same time: migratory flux between countries, control measure effects and dynamic model parameters fitted to each country. Then, as a second part of that work, Be-CoDiS was validated by considering the case of the 2014-16 West African Ebola Virus Disease (EVD) epidemic [15, 7, 16, 44]. EVD is a human and primate virus disease that causes a high mortality rate (between 50% and 90%) [14, 33]. During the period from December 2013 to March 2016, several important outbreaks were reported in West Africa (Guinea, Liberia, Sierra Leone and Nigeria). Furthermore, 16 isolated cases were detected in Mali, Senegal, the USA, the United Kingdom, Italy and Spain. The outbreak was considered over on March 29th 2016. It is estimated that around 28616 people were infected during those outbreaks and 11310 deaths were reported [43]. Starting with data from De- cember 2013, Be-CoDiS predicted (see [23]) a total of 28475 infected people, 11797 deaths and that the epidemic would end on April 19th, 2016. Once the model (available at: https://www.ucm.es/momat/epidemics) was de- veloped and validated, it would be good idea to study its mathematical properties. Its properties can be important for extracting conclusions that can be used in fur- ther analysis and/or developments of the model. In this article, we perform a stability analysis of two simplified versions of Be- CoDiS. To this aim, we first analyze the equilibrium states of the model by consid- ering only one country. More precisely, by applying the so-called next generation matrix method [36], we estimate an analytical expression of the disease basic repro- duction ratio [3, 12, 13], denoted by R0 according to the model parameters. The basic reproduction ratio associated to a disease free equilibrium, is used in epidemi- ology to determine the behavior of an epidemic. It is defined as the average number of new infections caused by one infected individual in a population in the conditions of a disease free equilibrium, over the course of its infectious period [1, 6]. We note that the mathematical definition of R0 used in this paper is specific to deterministic finite dimensional systems such as those considered here [36, 46]. It is generally expected that if R0 > 1 then the epidemic becomes endemic, whereas if R0 < 1 then the epidemic tends to a disease free equilibrium [1, 3]. What is expected for R0 = 1 is not always clear. Here, by considering particular techniques and results from [2, 34], we show that, if R0 ≤ 1, then the epidemic converges to a disease free equilibrium; otherwise, the epidemic does not disappear. We highlight the fact that the proposed approach is based on specific results not commonly used in the literature. Then, we extend this study to the case of two countries, when one country sends infected people to other country. Finally, we validate and illustrate the theoretical results obtained here, with numerical experiments based on data from the 2014- 16 West African Ebola virus epidemic and perform a sensitivity analysis of the estimated basic reproduction ratio, with respect to the model parameters. This work is organized as follows. In Section 2, we recall the formulation of the Be-CoDiS model presented in [23]. In Section 3, we study the equilibrium states of EJDE-2020/62 EPIDEMIOLOGICAL MODEL BE-CODIS 3 simplified versions of this model for one and two countries. In Section 4, considering data from the 2014-16 West African Ebola virus epidemic, we validate and illustrate the theoretical results with numerical experiments and perform a sensitivity analysis of the basic reproduction ratio with respect to the model parameters. Finally, in Section 5, we explain our conclusions. 2. Be-CoDiS model formulation We consider a disease with the following states for people (see [23, 26, 32, 33, 44]): • Susceptible (denoted by S): A person that is not infected by the disease pathogen. • Infected (denoted by E): A person that is infected by the disease pathogen but cannot infect other people and has no visible clinical signs (e.g., fever, haemorrhages, etc.). After an incubation period, the person passes to the infectious state. • Infectious (denoted by I): A person that can infect other people and start developing clinical signs. The mean duration of a person in this state is called the infectious period. After this period, infectious people are detected by sanitary authorities and hospitalized. • Hospitalized (denoted by H): A person that is hospitalized and can still infect other people. At the end of this state, the person can pass either to the Recovered state or to the Dead state. We point out that state H does not include hospitalized people which cannot infect other people any more. This last category of people is included in the Recovered state explained below. • Dead (denoted by D): A person that did not survived the disease. The cadavers of infected people can infect other people until they are buried. After a fixed average time, the body is buried. • Buried (denoted by B): A person that is dead because of the disease. His/her cadaver is buried and can no longer infect other people. • Recovered (denoted by R): A person that survived the disease, is no longer infectious and develop a natural immunity to the disease pathogen. After an infected person is hospitalized, the authorities may apply various control measures to control the disease spread (see [15, 19]): • Isolation: Infected people are isolated from contact with other people. Only sanitary professionals are in contact with them. However, depending on the considered disease, contamination of those professionals may also occur (see [15]). Isolated people receive an adequate medical treatment that reduces the disease fatality rate. • Quarantine: Movement of people in the area of origin of an infected person is restricted and controlled (e.g., quick sanitary check-points at the airports) to avoid potentia infected people spreading the disease. • Tracing: The objective of tracing is to identify potential infectious contacts which may have infected a person or spread the disease to other people. • Increase of sanitary resources: The number of operational beds and san- itary personal available to detect and treat affected people is increased. When necessary, the funerals of infected cadavers are controlled by sani- tary personal in order to reduce the contacts between the dead bodies and susceptible people. 4 B. IVORRA, D. NGOM, A. M. RAMOS EJDE-2020/62 Considering those general disease and control measures, the Be-CoDiS model is used to evaluate the spread of a human disease within and between countries during a fixed time interval. At the beginning of the simulation, the model parameters are set by the user. At the initial time (t = 0), only susceptible people live in the countries that are free of the disease, whereas the number of people in states S, E, I, H, R, D and B of the infected countries are set to their corresponding values. Then, during the time interval [0, Tmax], with Tmax ∈ IN being the maximum number of simulation days, within-country and between-country daily spread procedures, are applied. If at the end of a simulation day t, the number of people in state E, I, H and D is lower than a fixed threshold (smaller than 1), the simulation is stopped. Else, the simulation ends when t = Tmax. Furthermore, the control measures are also implemented and they can be activated or deactivated, when starting the model, in order to quantify their effectiveness to reduce the magnitude and duration of an epidemic. The spread of a disease within a particular country is modelled by using a de- terministic compartmental model (see [6]). For the sake of simplicity, we assume that, at each time, the population inside a country is homogeneously distributed and constant. Thus, the spatial distribution of the epidemic inside a country can be omitted. We also suppose that new births are susceptible people and the birth rate is, for the sake of simplicity, equal to the death rate (due to the disease or other causes). The disease spread between countries is modelled by using a spatial deterministic Individual Based model (see [11]), based on a population level transmission model. We consider that the flow of people between countries i and j at time t (i.e., people travelling per day from i to j at time t), is the only way to introduce the disease from country i, infected by the disease, to country j. To do so, we consider the matrix (τi,j) NCO i,j=1, where τi,j ∈ [0, 1] is the rate of transfer (day−1) of people from country i to country j, which is expressed in % of population in i per day. Furthermore, we assume that only people in states S and E can travel (we could also include here people in R, if it is suitable to do that), as other categories are not in condition to perform trips due to the clinical signs or to quarantine. Moreover, as a result of control measures in countries i and j, we assume that those rates can vary in time and are multiplied by a function denoted by mtr,i,j(t). Under those assumptions, the evolution of Si(t), Ei(t), Ii(t), Hi(t), Ri(t), Di(t) and Bi(t), denoting the number of susceptible, infected, infectious, hospitalized, recovered, dead and buried people in country i at time t, respectively, is modelled EJDE-2020/62 EPIDEMIOLOGICAL MODEL BE-CODIS 5 by the system of ordinary differential equations [23] dSi dt (t) = − Si(t) ( mI,i(t)βI,iIi(t) +mH,i(t)βH,iHi(t) +mD,i(t)βD,iDi(t) ) NPi(t) − µm,iSi(t) + µn,i ( Si(t) + Ei(t) + Ii(t) +Hi(t) +Ri(t) ) + ∑ i 6=j mtr,j,i(t)τj,iSj(t)− ∑ i 6=j mtr,i,j(t)τi,jSi(t), dEi dt (t) = Si(t) ( mI,i(t)βI,iIi(t) +mH,i(t)βH,iHi(t) +mD,i(t)βD,iDi(t) ) NPi(t) − µm,iEi(t) + ∑ i6=j mtr,j,i(t)τj,iXεfit(Ej(t)) − ∑ i 6=j mtr,i,j(t)τi,jXεfit(Ei(t))− γEXεfit(Ei(t)), dIi dt (t) = γEXεfit(Ei(t))− (µm,i + γI,i(t))Ii(t), dHi dt (t) = γI,i(t)Ii(t)− ( µm,i + (1− ωi(t))γHR,i(t) + ωi(t)γHD,i(t) ) Hi(t), dRi dt (t) = (1− ωi(t))γHR,i(t)Hi(t)− µm,iRi(t), dDi dt (t) = ωi(t)γHD,i(t)Hi(t)− γDDi(t), dBi dt (t) = γDDi(t), (2.1) where • i ∈ {1, . . . , NCO} is the index of each country. • NCO ∈ N is the number of countries. • NPi(t) = Si(t)+Ei(t)+Ii(t)+Hi(t)+Ri(t)+Di(t) is the number of people (alive and also died because of the disease) in country i at time t. • µn,i ∈ [0, 1] is the birth rate (day−1) in country i: the number of births per day and per capita. • µm,i ∈ [0, 1] is the mortality rate (day−1) in country i: the number of deaths per day and per capita (or, equivalently, the inverse of the mean life expectancy (day) of a person). • βI,i ∈ R+ is the disease effective contact rate (day−1) of a person in state I in country i: the mean number of effective contacts (i.e., contacts sufficient to transmit the disease) of a person in state I per day before applying control measures. • βH,i ∈ R+ is the disease effective contact rate (day−1) of a person in state H in country i. • βD,i ∈ R+ is the disease effective contact rate (day−1) of a person in state D in country i. • γE(i, t), γI,i(t) and γD(i, t) ∈ (0,+∞) denote the transition rate (day−1) from state E, I or D to state I, H or B, respectively: the number of people per day and per capita passing from one state to the other (or, equivalently, the inverse of the mean duration of one of those people in state E, I or D, respectively). We note that those parameters are time and country 6 B. IVORRA, D. NGOM, A. M. RAMOS EJDE-2020/62 dependent, since, due to the applied control measures in country i, their value could vary in time. • γHR,i(t) and γHD,i(t) denote the transition rate (day−1) from state H to state R (for people surviving the disease) and D (for people not surviving the disease), respectively. Here, we assume that the values of γHR,i(t) and γHD,i(t) are very close. Indeed, focusing on EVD, it has been observed that a person in state H passes to state R or S in around 4.5 days [26]. If values of γHR,i(t) and γHD,i(t) are different, the flow from I to H should be split into two sub flows, one from I to HR (hospitalized people surviving the disease), and other one from I to HD (hospitalized people not surviving the disease), with corresponding new flows from HR to R and from HD to D. • ωi(t) ∈ [0, 1] is the disease fatality rate in country i at time t: the proportion of people who do not survive the disease. • mI,i(t), mH,i(t), mD,i(t) ∈ [0, 1] (%) are functions representing the effi- ciency of the control measures applied to non-hospitalized infectious people, hospitalized people and infected cadavers respectively, in country i at time t to eradicate the outbreaks. Focusing on the application of the control measures, we multiply the disease contact rates (i.e., βI,i, βH,i and βD,i) by decreasing functions simulating the reduction of the number of effec- tive contacts as the control measures efficiency is improved. Here, we have considered the functions (see [27]): mI,i(t) = mH,i(t) = mD,i(t) = exp ( − κi max(t− λi, 0) ) , (2.2) where κi in [0,+∞) (day−1) simulates the efficiency of the control mea- sures (greater value implies lower value of disease contact rates) and λi in R ∪ {+∞} (day) denotes the first day of application of those control mea- sures. • Xεfit(x) = x if x ≥ εfit, Xεfit(x) = 2x− εfit if (εfit/2) ≤ x ≤ εfit, and 0 else- where, with εfit ≥ 0 being a small tolerance parameter. This function, with εfit small enough, is a filter used to avoid artificial spread of the epidemic due to negligible values of x. System (2.1) is completed with initial data Si(0), Ei(0), Ii(0), Hi(0), Ri(0), Di(0) and Bi(0) given in [0,+∞); for i ∈ {1, . . . , NCO}. Model (2.1) is summarized in Figure 1. Remark 2.1. We note that the Be-CoDiS model proposed here is not only limited to the study of the EVD, but can also tackle other diseases, such as the Middle East respiratory syndrome-related coronavirus (MERS) or the Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [9, 24, 21], by adapting the compartmental model and its parameters. 3. Analytical behavior of the Be-CoDiS model Here, we are interested in studying the equilibrium states and in estimating the basic reproduction ratio of simplified versions of the Be-CoDiS model presented in Section 2. First, we focus on the case of one country with an emigration flow of susceptible or infected people and an immigration flow of susceptible people. Then, EJDE-2020/62 EPIDEMIOLOGICAL MODEL BE-CODIS 7 γ I E n m D γ(1−ω) HR D m µ (i,j)τ RHIES µ γHI m m HI tr B γ β + β + β D m D γ HDω Migratory flow with other countries Natural mortality/natality dynamic Control measures Figure 1. Summary of the complete version of the Be-CoDiS model. we extend the study to the case of two countries, with one country sending people to the other one. 3.1. Simplified model for one country. Here, we are interested in studying the behavior of the epidemic inside one single country (i.e., we omit to study the spread in other countries). Additionally, we assume that this country send infected people to other countries, but does not receive infected people from other countries. Those assumptions are suitable if we consider that this country is the index case (i.e., the country of origin of the disease) and the spread of the disease to other countries is not expected. For the sake of simplicity, we assume that the population size in the considered country is constant and equal to N ∈ N (i.e., emigration or death flows are compensated by immigration or birth flows entering the susceptible state). This hypothesis is reasonable as, due to the size of the population in a country (generally larger than a million of people) and the time scale of the study (generally lower than five years) considered here, the global variation of the population size during a simulation is negligible [20]. To avoid asymptotic endemic solutions with the whole population concentrated in state B, the flow from state D to state B is replaced by a flow from state D to state S and state B is omitted. We note that this change in the model satisfies the hypothesis mentioned previously that deaths are replaced by immigration or births in the susceptible state in order to keep constant the size of the population. Furthermore, to simplify notations, we consider that S, E, I, H, R and D denote the ratio of people in each state in the considered country (rather than the total number of people). Additionally, we assume that the model coefficients are constant and no control measures are applied. As no other country is considered, the filter Xεfit is omitted. A diagram of this model for one country is shown in Figure 2. Under these assumptions, the evolution of the epidemic, is modelled by dS dt (t) = −S(t) ( βII(t) + βHH(t) + βDD(t) ) + τE(t) + µ ( E(t) + I(t) +H(t) +R(t) ) + θD(t), (3.1) dE dt (t) = S(t) ( βII(t) + βHH(t) + βDD(t) ) − (µ+ δ + τ)E(t), (3.2) dI dt (t) = δE(t)− (µ+ γ)I(t), (3.3) 8 B. IVORRA, D. NGOM, A. M. RAMOS EJDE-2020/62 τ τ Θ Θ λ αγβ + β + β D RIES µ H I H D δ Figure 2. Simplified model with one country considered in Sec- tion 3.1. dH dt (t) = γI(t)− ( µ+ λ+ α ) H(t), (3.4) labelmod dR dt (t) = αH(t)− µR(t), dD dt (t) = λH(t)− θD(t), (3.5) where • µ ∈ [0, 1] is the mortality rate (day−1). • βI ∈ R+ is the disease effective contact rate (day−1.person−1) of people in state I. • βH ∈ R+ is the disease effective contact rate (day−1.person−1) of people in state H. • βD ∈ R+ is the disease effective contact rate (day−1.person−1) of people in state D. • δ, γ and θ denote the transition rates (day−1) from state E to I, I to H and D to S, respectively. • λ is the disease fatality percentage times the transition rate from state H to state D. • α is the disease survival percentage (1 minus the disease fatality percentage) times the transition rate from state H to state R. • τ ∈ [0, 1] is the daily rate (%) of the movement of people in states S or E. We point out that, to keep constant the population, people leaving the country for unit of time (i.e., τ(S +E)) and people being buried for unit of time (i.e., θD) are added to the susceptible state. The main parameters used in this work and their corresponding range of values used in Section 4 are summarized in Table 1. For convenience, we write the solutions of (3.5) as vectors (E(t), I(t), H(t), D(t), S(t), R(t)) ∈ [0, 1]6, for all t ≥ 0. We also consider Ω = {(E, I,H,D, S,R) ∈ [0, 1]6 : E + I +H +D + S +R = 1}. Theorem 3.1. The set Ω is positively invariant for system (3.5) (i.e., if (E(0), I(0), . . . , R(0)) ∈ Ω, then (E(t), I(t), . . . , R(t)) ∈ Ω, for all t > 0). Proof. First, we note that system (3.5) is positive (i.e., if (E(0), . . . , R(0)) ∈ [0,+∞)6, then (E(t), . . . , R(t)) ∈ [0,+∞)6, for all t > 0). Indeed, if (E(t), . . . , R(t)) ∈ [0,+∞)6 and S(t) = 0, then dS(t) dt ≥ 0, which guarantees that S cannot become EJDE-2020/62 EPIDEMIOLOGICAL MODEL BE-CODIS 9 Table 1. Summary of the main notations used in this work. A brief description (Description) and the range (i.e., minimum and maximum values) of the values (Range of Value) used in Section 4, with i ∈ {1, 2}. Notation Range of Value Description βIi [0.0494,0.2671] Disease effective contact rate (day−1person−1) of people in state I in country i βHi [0.0107,0.0200] Disease effective contact rate (day−1person−1) of people in state H in country i βDi [0.0494,0.2671] Disease effective contact rate (day−1person−1) of people in state D in country i δi [0.0120,0.0230] Transition rate (day−1) from state E to state I in country i γi [0.2000,0.5000] Transition rate (day−1) from state I to state H in country i αi [0.1050,0.1480] Disease survival proportion times transition rate (day−1) from state H to state R λi [0.0328,0.1282] Disease fatality rate times transition rate (day−1) from state H to state D θi [0.5000,1.0000] Transition rate (day−1) from state D to state S in country i µi [0.012,0.023] Natural mortality rate (day−1) in country i τi [0,2.4]·10−5 Daily rate (% day−1) of the movement of peo- ple leaving country i Ni [10,20]·106 Number of people in country i Si, Ei, Ii, Hi, Ri, Di [0,1] proportion (%) of people in state S, E, I, H, R, D in country i at time t negative. This property is also true for the other disease states, ensuring the posi- tivity of the considered system. Additionally, since dE dt + dI dt + dH dt + dD dt + dS dt + dR dt = 0, we have that E(t)+I(t)+ H(t) +D(t) + S(t) + R(t) = E(0) + I(0) +H(0) +D(0) + S(0) + R(0) = 1 for all t ≥ 0. Thus, we deduce that Ω is positively invariant for system (3.5). � For the study of the stability properties of system (3.5), we will use the basic reproduction ratio R0, which is the average number of secondary cases produced by one infected individual during its entire infectious period in an otherwise uninfected population [1, 18]. From a mathematical point of view, the value of R0 associated to the epidemi- ological compartmental model (3.5) can be computed as the spectral radius of the so-called next generation matrix (see [36] for more details) as explained below. More precisely, let us consider a general compartmental model for infectious disease transmission, defined (following the notation used in [34]) by Ẋ = F(X,Y )− V(X,Y ) Ẏ = g(X,Y ), (3.6) where vectors X = (x1, . . . , xn)T ∈ Rn and Y = (y1, . . . , ym)T ∈ Rm represent the populations in infected and non-infected states, respectively; vector F = (F1, 10 B. IVORRA, D. NGOM, A. M. RAMOS EJDE-2020/62 F2, . . . ,Fn)T , with Fi being the rate of appearance of news infections in com- partment i; vector V = (V1,V2, . . . ,Vn)T , with Vi being the rate of transfer of individuals out of (for positive values) or into (for negative values) compartment by all other means; and vector g = (g1, . . . , gm)T represents the transition terms for non-infected states. Following [34, 37], we assume that: (A1) All functions Fi, Vi and gi are C2([0, 1]n+m;R). (A2) F(0, Y ) = V(0, Y ) = 0 (if the population is free of disease then it will remain free of disease). (A3) Fi(X,Y ) ≥ 0 for all i ∈ {1, . . . , n} if X,Y satisfy xi ≥ 0 and yj ≥ 0 for all (i, j) ∈ {1, . . . , n} × {1, . . . ,m}. (A4) Given i ∈ {1, . . . , n}, Vi(X,Y ) ≤ 0 if vector X satisfies xi = 0; Given j ∈ {1, . . . ,m}, gj(X,Y ) ≥ 0 if vector Y satisfies yj = 0, (if a compartment is empty, then there can be no transfer of individuals out of the compartment). (A5) ∑n i=1 Vi (X,Y ) ≥ 0 for all X,Y such that xi ≥ 0 and yj ≥ 0, for all (i, j) ∈ {1, . . . , n} × {1, . . . ,m} (the total outflow from all infected compartments is non negative). (A6) The disease-free system Ẏ = g(0, Y ) has a unique equilibrium Yf ∈ ΩY = {Y ∈ [0, 1]m : Y1 + · · ·+ Ym = 1} which is globally asymptotically stable in ΩY . Therefore, using assumptions (A1)–(A6), Pf = (Xf , Yf), with Xf = (0, . . . , 0), is the unique admissible disease free equilibrium for system (3.6); we refer to this point as the disease-free equilibrium. Furthermore, there will be a unique basic reproduction ratio R0, the one associated to Pf , which can be computed as follows: Let F = [∂Fi ∂xj (Pf) ]n i,j=1 and V = [∂Vi ∂xj (Pf) ]n i,j=1 . (3.7) From assumptions (A1)–(A6) we have (see [37, page 174]) that, if V is nonsingu- lar and F and V −1 are non-negative matrices (i.e., square matrices all of whose elements are nonnegative), then the basic reproduction ratio associated with Pf is given by R0 = ρ(FV −1), the spectral radius of matrix FV −1 (see [36], page 33), which is the so called next generation matrix. Taking into consideration this result, we introduce the following formulation of system (3.5). Let P = (X,Y )T , with X = (E, I,H,D)T and Y = (S,R)T . system (3.5) can be rewritten as system (3.6), where F(X,Y ) = F̃(S)X, V(X,Y ) = ṼX and g(X,Y ) = g̃(S) ( X Y ) , with F̃(S) =  0 βIS βHS βDS 0 0 0 0 0 0 0 0 0 0 0 0  , Ṽ =  (µ+ δ + τ) 0 0 0 −δ (µ+ γ) 0 0 0 −γ (µ+ λ+ α) 0 0 0 −λ θ  , g̃(S) = ( µ+ τ µ− βIS µ− βHS θ − βDS 0 µ 0 0 α 0 0 −µ ) . EJDE-2020/62 EPIDEMIOLOGICAL MODEL BE-CODIS 11 Theorem 3.2. The basic reproduction ratio of system (3.5) is R0 = δ(αθβI + γλβD + γθβH + λθβI + µθβI) (µ+ δ + τ)(µ+ γ)(µ+ λ+ α)θ . (3.8) Proof. To compute the basic reproduction ratio of system (3.5), we apply the Next Generation Matrix methodology briefly described above (see [36]). We first rewrite (as done above) system (3.5) as (3.6) and notice that the assumptions (A1)–(A5) are satisfied. To check (A6), we consider the system Ẏ = g(0, Y ) which is given by dS dt (t) = µR(t), dR dt (t) = −µR(t). (3.9) In ΩY = {(S,R)T ∈ [0, 1]2 : S +R = 1}, the solutions of system (3.9) satisfy dS dt (t) = µR(t) = µ(1− S(t)). Thus, (S(t), R(t))T = (1 − R(0)e−µt, R(0)e−µt) is the unique solution of (3.9) in ΩY , once an initial value R(0) ∈ [0, 1] is given, and Yf = (1, 0)T is its unique equilibrium point in ΩY , which is globally asymptotically stable. Therefore, all the assumptions (A1)–(A6), detailed above, are satisfied. Matrices F = [∂Fi ∂xj (Pf)] 4 i,j=1 and V = [∂Vi∂xj (Pf)] 4 i,j=1 are given by F =  0 βI βH βD 0 0 0 0 0 0 0 0 0 0 0 0  , V = Ṽ =  (µ+ δ + τ) 0 0 0 −δ (µ+ γ) 0 0 0 −γ (µ+ λ+ α) 0 0 0 −λ θ  , (3.10) respectively. Furthermore, V is nonsingular, V −1 =  (µ+ δ + τ)−1 0 0 0 δ (µ+γ)(µ+δ+τ) (µ+ γ)−1 0 0 γδ (µ+λ+α)(µ+γ)(µ+δ+τ) γ (µ+λ+α)(µ+γ) (µ+ λ+ α)−1 0 λ γ δ (µ+λ+α)(µ+γ)(µ+δ+τ)θ λ γ (µ+λ+α)(µ+γ)θ λ (µ+λ+α)θ θ−1  , and FV −1 =  δ(αθβI+γλβD+γθβH+λθβI+µθβI) (µ+δ+τ)(µ+γ)(µ+λ+α)θ βIθ(µ+λ+α)+βHγθ+βDλ γ (µ+λ+α)(µ+γ)θ βHθ+βDλ (µ+λ+α)θ βD θ 0 0 0 0 0 0 0 0 0 0 0 0  It is clear that the eigenvalues of FV −1 are δ(αθβI + γλβD + γθβH + λθβI + µθβI) (µ+ δ + τ)(µ+ γ)(µ+ λ+ α)θ (simple) and 0 (whose algebraic multiplicity is 3). Therefore, as claimed previously, R0 = ρ(FV −1) = δ(αθβI + γλβD + γθβH + λθβI + µθβI) (µ+ δ + τ)(µ+ γ)(µ+ λ+ α)θ , which is the value given in (3.8). � 12 B. IVORRA, D. NGOM, A. M. RAMOS EJDE-2020/62 Using the basic reproduction ratio R0 given in (3.8), we have the following sta- bility results. Theorem 3.3. With the notation used in (3.6), if (A1)–(A6) are satisfied, we have the following results: (1) The (unique) disease free equilibrium Pf = (0, 0, 0, 0, 1, 0) of system (3.5) is globally asymptotically stable in Ω if R0 ≤ 1, and unstable if R0 > 1. (2) If R0 > 1, system (3.5) has, at least, one endemic disease equilibrium Pe = (Xe, Ye) ∈ Ω, with Xe = (Ee, Ie, He, De) and Ye = (Se, Re) given by Se = 1 R0 , Ee = θµ(µ+ γ)(µ+ α+ λ)φ, Ie = δθµ(µ+ α+ λ)φ, He = δθγµφ, Re = δθαγφ, De = δγλµφ φ = 1 (δγλ(µ− θ) + (µ+ δ + τ)(µ+ γ)(µ+ λ+ α)θ) ( 1− 1 R0 ) . Before proving Theorem 3.3 we recall some results that will be used in the proof. We consider a general autonomous differential equation ẋ(t) = f(x(t)), t > 0, (3.11) such that x ∈ Rn and f ∈ C0(Rn,Rn). Let us denote, for any function V ∈ C1(U,R) with U ⊂ Rn, V̇ (x) = ∇V (x)Tf(x). (3.12) Theorem 3.4. Let G ⊂ Rn be a compact and positively invariant set for system (3.11). Let L ∈ C1(G,Rn) such that L̇(x) ≤ 0, for all x ∈ G. We consider the following sets: S1 = {x ∈ G : L(x) = minw∈G L(w)}; S2 = {x ∈ G : L̇(x) = 0} and S3 the largest invariant set in S2 for system (3.11). If S1 = S3, then, starting from any point in G, system (3.11) converges asymptotically to S3. For a proof of the above theorem, see, for instance, [2, page 346]. Proof of Theorem 3.3. First, we determine the equilibrium states of system (3.6) by solving (F̃(S)− Ṽ)X = 0 and g̃(S)P = 0. After some computation (with Maple 16), it can be proved that Pf ∈ Ω and Pe, given in the statement of Theorem 3.3, are equilibrium points. If R0 > 1, Pe ∈ Ω, else R0 ≤ 1, Pe /∈ Ω (we also note that, if R0 > 1 and tends to 1, then Pe tends to Pf). Let us assume that R0 ≤ 1: To build a function L as defined in the statement of Theorem 3.4, we use a method developed in [34] to determine a Lyapunov function for the disease free points of system (3.6) (i.e., points such that E = I = H = D = 0). With that aim, the first line of system (3.6) is rewritten as Ẋ = (F − V )X − f(X,Y ), where F and V are defined in (3.7) and f(X,Y ) = (F −V )X−F(X,Y )+V(X,Y ). Thus, in our case, from (3.10), we have that f(X,Y ) = (F − V )X − F̃(S)X + ṼX = (F − F̃(S))X. Matrices FV −1 and V −1F have the same eigenvalues (we remind that the prod- ucts of AB and BA of two arbitrary matrices A ∈Mn×m and B ∈Mm×n have the same set of eigenvalues since, if u is an eigenvector of AB, then Bu is an eigenvector EJDE-2020/62 EPIDEMIOLOGICAL MODEL BE-CODIS 13 of BA of the same eigenvalue), which are (as seen previously) R0 (simple) and 0. Furthermore, it can be easily proved (for instance with symbolic computations done with Maple) that w = (0, βI , βH , βD) is a left eigenvector of V −1F corresponding to the eigenvalue R0. Let Lf : R4 × R2 → R given by Lf(X,Y ) = wV −1X, which leads, by simple calculations, to Lf(X,Y ) = ( βIδ (µ+ γ)(µ+ δ + τ) + βHγ δ (µ+ δ + τ)(µ+ λ+ α)(µ+ γ) + βDλ γ δ (µ+ λ+ α)(µ+ γ)(µ+ δ + τ)θ ) E + ( βI (µ+ γ) + βHγ (µ+ γ)(µ+ λ+ α) + βDλ γ (µ+ λ+ α)(µ+ γ)θ ) I + ( βH (µ+ λ+ α) + βDλ (µ+ λ+ α)θ ) H + βD θ D. We note that Lf is non negative in Ω and, according to (3.12), L̇f(X,Y ) = wV −1((F − V )X − f(X,Y )) = (R0 − 1)wX − wV −1f(X,Y ). Additionally, since R0 ≤ 1 and the coordinates of V −1 and f(X,Y ) are non negative for all (X,Y ) ∈ Ω then L̇f(X,Y ) ≤ 0 for all (X,Y ) ∈ Ω. Let Ωf = {(E, I,H,D, S,R) ∈ Ω : E = I = H = D = 0}. If x ∈ Ωf , then Lf(x) = 0, and if x ∈ Ω\Ωf then Lf(x) > 0. Moreover, we note that L̇f(X,Y ) = 0 if and only if (R0 − 1)wX = 0, and wV −1f(X,Y ) = ( − βIδ (µ+ γ)(µ+ δ + τ) − βHγ δ (µ+ δ + τ)(µ+ λ+ α)(µ+ γ) − βDλ γ δ (µ+ λ+ α)(µ+ γ)(µ+ δ + τ)θ )( 1− S )( βII + βHH + βDD ) = 0 From the first condition, we have that I = H = D = 0 if R0 < 1 (no condition if R0 = 1). From the second condition, we have that S = 1 (and E = I = H = D = R = 0) or I = H = D = 0. Using the notation introduced in Theorem 3.4, we have that • S1 = {(X,Y ) ∈ Ω : Lf(X,Y ) = minW∈Ω Lf(W ) = 0} = Ωf . • S2 = {(X,Y ) ∈ Ω : L̇f(X,Y ) = 0} = {(X,Y ) ∈ Ω : I = H = D = 0}. • S3 = Ωf . Indeed, on the one hand, starting from any point in S2\Ωf , we have that E(0) > 0 and I(0) = 0, thus, due to the third line of system (3.5), dI dt (0) > 0 which implies that the trajectory gets out of S2. On the other hand, the solution of system (3.5) starting from any point in Ωf , say, (0, 0, 0, 0, 1−R(0), R(0)), with R(0) ∈ [0, 1], is given by (E(t), I(t), H(t), D(t), S(t), R(t)) = (0, 0, 0, 0, 1−R(t), R(0) exp(−µt)) ∈ Ωf for any t ≥ 0. Thus, Ωf is is positively invariant for system (3.5) and the largest invariant set in S2 is S3 = Ωf . Since all the hypothesis of Theorem 3.4 are fulfilled, we deduce that system (3.5), starting from any point in Ω, converges asymptotically to Ωf . 14 B. IVORRA, D. NGOM, A. M. RAMOS EJDE-2020/62 In particular, this implies that limt→+∞H(t) = 0. From the fifth line of system (3.5) we have that dR dt (t) = αH(t)− µR(t), whose solution is given by R(t) = R(0)e−t + e−t ∫ t 0 H(x)exdx. Let K(t) = ∫ t 0 H(x)ex. Since H(x)ex ≥ 0 for any x ∈ R, K(t) is a non-decreasing and non-negative function. Thus, we have two cases: (1) limt→+∞K(t) ∈ R which implies that limt→+∞R(t) = 0. (2) limt→+∞K(t) = +∞. Then, applying L’Hôpital’s rule, we obtain that lim t→+∞ R(t) = lim t→+∞ e−t ∫ t 0 H(x)exdx = lim t→+∞ −e−tH(t)et = lim t→+∞ −H(t) = 0. This implies limt→+∞R(t) = 0 and limt→+∞ S(t) = 1. We conclude that Pf is globally asymptotically stable in Ω for system (3.5). Let us assume that R0 > 1. System (3.6) satisfies the hypothesis in [36, Theorem 2 ] and, thus, the equilibrium point Pf is unstable when R0 > 1. � 3.2. Simplified model for 2 countries. In this section, we study the epidemi- ological behavior of two countries. One of them (denoted by Country 2) receiving infected people from the other country (denoted by Country 1). We suppose no extra control measures are taken (so that mi,j = 1) and take into account the same assumptions and notations (but indexed by i = 1, 2 according to the country) as those introduced in Section 3.1. For the sake of simplicity, we assume that each country has constant (in time) population, N1 and N2, respectively, and movement of people from Country 1 to Country 2. This is a reasonable assumption since, for the typical duration of a EVD epidemic, the population of a country does not vary significantly. Remark 3.5. The same formulation results for the case of two countries keeping constant populations, with movements from Country 1 to Country 2 and viceversa, with no movements from the infected state in Country 2 (E2) to Country 1. As usual (see [23]) only people in the susceptible (S) or infected (E) states are supposed to travel. The proportion of susceptible and infected people travelling, at time t, from Country 1 to Country 2 per unit time is τ1S1(t) and τ1E1(t), respec- tively. To keep constant populations in both countries, this rate of people is also reintroduced in the susceptible state of Country 1 and removed (after normalizing with the population of Country 2) from the susceptible state of the second coun- try. Furthermore, to avoid removing too many people from the susceptible state in Country 2 in the extreme case (and clearly unrealistic) of having E1(t) much bigger than S2(t), we set a limit K (big enough) for the ratio, so that when E1(t) S2(t) ≥ K, only τ̃1KS2(t) of τ̃1E1(t) is removed from the susceptible state S2 per unit time and the rest, τ̃1(E1(t) −KS2(t)), is removed from the infected state E2, per unit time. Here, τ̃1 = τ1 N1 N2 represents the proportion of persons going from Country 1 to Country 2 per unit of time and relative to the population size N2. More precisely, τ1E1(t) is leaving the infected state E1 and entering in the suscep- tible states S1 per unit time, in Country 1. In Country 2, τ̃1E1(t) is entering in the infected state E2(t) per unit time. Furthermore, to keep constant the population N2, EJDE-2020/62 EPIDEMIOLOGICAL MODEL BE-CODIS 15 • τ̃1χK(E1, S2) is removed from the susceptible state S2 per unit time, where χK : [0, 1]2 → R is the continuous function defined by χK(x, y) = min{x,Ky} = { x, if x ≤ Ky Ky, if x ≥ Ky . • The following quantity is removed from the infected state E2 per unit time: τ̃1 max{E1 −KS2, 0} = { 0, if E1 ≤ KS2 τ̃1(E1 −KS2), if E1 ≥ KS2 We observe that the final balance of the proportion of infected people entering the infected state in Country 2 is τ̃1(E1 −max{E1 −KS2, 0}) = τ̃1 min{E1,KS2} = τ̃1χK(E1, S2). Furthermore, following an idea in [23] to avoid unrealistic spread of the epidemic because of unrealistic negligible values of movement of people in the state E from one country to another, we consider only the contribution of infected individuals from Country 1 in Country 2, when the number of infected individuals in Country 1 (given by N1E1) is greater than a given threshold N1ε > 0. To do that, we set ε ≥ 0 (small enough) and change the function χK by the continuous function χK,ε : [0, 1]2 → R defined by χK,ε(x, y) =  x, if x ≤ Ky and x ≥ ε, 2x2 ε − x, if x ≤ Ky and ε ≥ x ≥ ε/2, Ky, if x ≥ Ky and x ≥ ε 2Ky ε x−Ky, if x ≥ Ky and ε ≥ x ≥ ε/2, 0, if x ≤ ε/2. (3.13) This function is a filter used to avoid artificial spread of the epidemic due to negligible values of E1(t) and to keep nonnegativity of the function S2(t) and constant population of the second country. We point out that, when ε = 0, χK,0(E1(t), S2(t)) = χK(E1(t), S2(t)). A particular graphical representation of this function is depicted in Figure 3. A representation of the distribution of the analytical formulation of χK,ε in the plane OXY is presented in Figure 4. 16 B. IVORRA, D. NGOM, A. M. RAMOS EJDE-2020/62 Figure 3. Representation of the function χ1.5,0.6(x, y). P2(x,y) P1(x,y) Ky x 0 0 0 εε/2 y=x/K y x 1 1 Figure 4. Distribution of the analytical formulation of χK,ε in the plane OXY . Here, P1(x, y) = 2Ky ε x − Ky and P2(x, y) = 2x2 ε − x. Taking into account these hypothesis, we now consider the two systems dS1 dt (t) = −S1(t) ( βI,1I1(t) + βH,1H1(t) + βD,1D1(t) ) + τ1E1(t) + µ1 ( E1(t) + I1(t) +H1(t) +R1(t) ) + θ1D1(t), dE1 dt (t) = S1(t) ( βI,1I1(t) + βH,1H1(t) + βD,1D1(t) ) − (µ1 + δ1 + τ1)E1(t), dI1 dt (t) = δ1E1(t)− (µ1 + γ1)I1(t), dH1 dt (t) = γ1I1(t)− ( µ1 + λ1 + α1 ) H1(t), dR1 dt (t) = α1H1(t)− µ1R1(t), dD1 dt (t) = λ1H1(t)− θ1D1(t), (3.14) EJDE-2020/62 EPIDEMIOLOGICAL MODEL BE-CODIS 17 and dS2 dt (t) = −S2(t) ( βI,2I2(t) + βH,2H2(t) + βD,2D2(t) ) − τ̃1XK,ε(E1(t), S2(t)) + τ2E2(t) + µ2 ( E2(t) + I2(t) +H2(t) +R2(t) ) + θ2D2(t), dE2 dt (t) = S2(t) ( βI,2I2(t) + βH,2H2(t) + βD,2D2(t) ) − (µ2 + δ2 + τ2)E2(t) + τ̃1XK,ε(E1(t), S2(t)), dI2 dt (t) = δ2E2(t)− (µ2 + γ2)I2(t), dH2 dt (t) = γ2I2(t)− ( µ2 + λ2 + α2 ) H2(t), dR2 dt (t) = α2H2(t)− µ2R2(t), dD2 dt (t) = λ2H2(t)− θ2D2(t), (3.15) where all the constants involved are positive and denoting similar things as those of systems (2.1) and (3.5). A diagram summarizing this model is presented in Figure 5. 1 τ X ∼ Κ,ε Κ,ε ∼ Xτ 1 2Θ 1Θ 2 2 2 222 2 1 1 11 1 1I1 D2H2I2 D1H1 2 22222 1 11H111 µ S E I R D S E I H R D β + β + β δ γ α λ Θ µ γ α λ Θ β + β + β δ τ τ1 Figure 5. Simplified model with two countries considered in Sec- tion 3.2. Let Ω2 = Ω × Ω = {(E1, I1, H1, D1, S1, R1, E2, I2, H2, D2, S2, R2) ∈ [0, 1]12 : E1 + I1 +H1 +D1 + S1 +R1 = 1 and E2 + I2 +H2 +D2 + S2 +R2 = 1}. Theorem 3.6. The set Ω2 is positively invariant for system (3.14)-(3.15). Proof. Let (E1(0), I1(0), . . . , R2(0)) ∈ Ω2. Since (E1(0), . . . , R1(0)) ∈ Ω and sys- tem (3.14) has the same structure as system (3.5), by Theorem 3.1, the solution (E1(t), I1(t), H1(t), D1(t), S1(t), R1(t)), which is governed by system (3.14), re- mains in Ω for all t ∈ R. 18 B. IVORRA, D. NGOM, A. M. RAMOS EJDE-2020/62 Moreover, (E2(0), . . . , R2(0)) ∈ Ω implies that (E2(t), . . . , R2(t)) ∈ Ω for all t ≥ 0. Indeed, if S2(t) = 0, then dS2 dt (t) ≥ 0, which guarantees that S2 cannot become negative. Same reasoning applies for E2, I2, H2, D2 and R2. Additionally, since dE2 dt (t) + dI2 dt (t) + dH2 dt (t) + dD2 dt (t) + dS2 dt (t) + dR2 dt (t) = 0, we have that E2(t)+I2(t)+H2(t)+D2(t)+S2(t)+R2(t) = E2(0)+I2(0)+H2(0)+ D2(0)+S2(0)+R2(0) = 1 for all t ≥ 0. Thus, (E2(t), I2(t), H2(t), D2(t), S2(t), R2(t)) remains in Ω for all t ∈ R. Thus, if (E1(0), I1(0), . . . , R2(0)) ∈ Ω2, then for all t ∈ R (E1(t), I1(t), . . . , R2(t)) ∈ Ω2. � Following the process in Theorem 3.2, we consider R0,i = δi(αiθiβI,i + γiλiβD,i + γiθiβH,i + λiθiβI,i + µiθiβI,i) (µi + δi + τi)(µi + γi)(µi + λi + αi)θi , (3.16) with i = 1, 2, which is helpful for studying the stability of system (3.14)-(3.15). Theorem 3.7. System (3.14)–(3.15) admits a disease free equilibrium Pf = (Pf,1, Pf,2), with Pf,1 = Pf,2 = (0, 0, 0, 0, 1, 0). Furthermore, the following results hold: (1) If R0,1 ≤ 1 and R0,2 ≤ 1, then Pf is globally asymptotically stable. (2) If R0,1 ≤ 1 and R0,2 > 1, then the solution (E1(t), I1(t), H1(t), D1(t), S1(t), R1(t)) of sub-system (3.14) with any initial data (E1(0), . . . , R(0)) ∈ Ω, tends to the disease free state Pf,1 = (0, 0, 0, 0, 1, 0) as t → ∞, and there exists an endemic equilibrium Pe,2 = (Ee,2, Ie,2, He,2, De,2, Se,2, Re,2) ∈ Ω for sub-system (3.15). (3) If R0,1 > 1, then Pf,1 is an unstable disease free equilibrium and there exists an endemic equilibrium Pe,1 = (Ee,1, Ie,1, He,1, De,1, Se,1, Re,1) for sub-system (3.14). Additionally, if (E1, . . . .R1, E2, . . . , R2) is a solution of (3.14)-(3.15) with (E1(0), . . . , R1(0), E2(0), . . . , R2(0)) ∈ Ω2 and KS2(t) > E1(t) > ε for all t > 0 (which is a reasonable assumption in real cases), then E2(t) does not converge to 0 as t→∞. Here, for i = 1, 2, we have Se,i = 1 R0,i , Ee,i = θiµi(µi + γi)(µi + αi + λi)φi, Ie,i = δiθiµi(µi + αi + λi) φi, He,i = δiθiγiµiφi, Re,i = δiθiαiγiφi and De,i = δiγiλiµiφi, with φi = 1 (δiγiλi(µi − θi) + (µi + δi + τi)(µi + γi)(µi + λi + αi)θi) ( 1− 1 R0,i ) . Proof. It is obvious that Pf is always a disease free equilibrium. Let us prove the other points of the theorem. (1) and (2): We first assume that R0,1 ≤ 1. Since sub-system (3.14) is inde- pendent of sub-system (3.15) and similar to system (3.5), as proven in Theorem 3.3, (E1(t), I1(t), H1(t), D1(t), S1(t), R1(t)) converges to the disease free equilibrium Pf,1. This implies that there exists a time tε > 0, such that E1(t) < ε 2 , for all t > tε. Thus, for all t > tε, XK,ε(E1(t), S2(t)) = 0 and the first and second equations of sub-system (3.15) are given by dS2 dt (t) = −S2(t) ( βI,2I2(t) + βH,2H2(t) + βD,2D2(t) ) EJDE-2020/62 EPIDEMIOLOGICAL MODEL BE-CODIS 19 + τ2E2(t) + µ2 ( E2(t) + I2(t) +H2(t) +R2(t) ) + θ2D2(t), dE2 dt (t) = S2(t) ( βI,2I2(t) + βH,2H2(t) + βD,2D2(t) ) − (µ2 + δ2 + τ2)E2(t). Therefore, sub-system (3.15) is equivalent to system (3.5) for t > tε and, as proven in Theorem 3.3: • If R0,2 ≤ 1 , then (E2(t), I2(t), H2(t), D2(t), S2(t), R2(t)) converges to the disease free equilibrium Pf,2 as t→∞. • If R0,2 > 1, then the point Pe,2 given in the statement of Theorem 3.7 is an endemic equilibrium for sub-system (3.15). (3) Assuming R0,1 > 1, from Theorem 3.3, we deduce that Pf,1 is an unstable disease free equilibrium, the point Pe,1 given in the statement of Theorem 3.7 is an endemic equilibrium point of sub-system (3.14) and Ee,1 > 0. Let us assume that for all t > 0, KS2(t) > E1(t) > ε. By reductio ad absurdum, if limt→+∞E2(t) = 0, then there exists t1 > 0 such that for all t > t1, E2(t) < τ̃1ε 2(µ2+δ2+τ2) . Additionally, from the second equation of sub-system (3.15), dE2(t) dt ≥ −(µ2 + δ2 + τ2)E2(t) + τ̃1XK,ε(E1(t), S2(t)) > − τ̃1ε 2 + τ̃1E1(t) > τ̃1ε 2 for all t > t1. This implies that limt→+∞E2(t) = +∞, which is not possible because, as said previously 0 ≤ E2(t) ≤ 1 for all t ≥ 0 (see Theorem 3.6). Thus, E2(t) does not converge to 0 as t→∞. � Remark 3.8. From Theorem 3.7, we can define a basic reproduction ratio for the disease described by System (3.14)-(3.15) by considering R0 = max(R0,1, R0,2). In- deed, if R0 ≤ 1, System (3.14)-(3.15) converges globally and asymptotically to the disease free equilibrium (0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0), else, under reasonable hypoth- esis, it does not converges to this disease free state. 4. Application to the 2014-2016 West African EVD epidemics In this section, to validate and illustrate the interest of the theoretical results obtained previously, we present some numerical experiments based on data from the 2014-2016 West African EVD epidemics. To do so, in subsection 4.1, we perform a sensitivity analysis of the basic reproduction ratio studied in section 3.1, regarding the model parameters. This sensitivity analysis will be used later, in Section 5, to propose strategies to allocate the resources for fighting EVD. Next, in section 4.2, to exhibit the stability results highlighted in Theorem 3.7, we present the evolution of the epidemic between two countries by considering several sets of parameters. 4.1. Sensitivity analysis of the basic reproduction ratio. In Table 2, we show the maximum and minimum values of the parameters of System (3.5) reported in the literature for the 2014-2016 West African EVD case [7, 15, 16, 23, 32]. Considering those values, we study the impact of variations in each model pa- rameter on the value of the basic reproduction ratio given in Theorem 3.3 and rewritten as a function R0(P ) = δ(αθβI + γλβD + γθβH + λθβI + µθβI) (µ+ δ + τ)(µ+ γ)(µ+ λ+ α)θ , (4.1) with P = (βI , βD, βH , δ, α, θ, γ, λ, µ, τ) = (P1, P2, . . . , P10) ∈ (0,+∞)9 × [0,+∞). 20 B. IVORRA, D. NGOM, A. M. RAMOS EJDE-2020/62 Table 2. Minimum and maximum values of the parameters of system (3.5) for the 2014-2016 West African EVD case. Parameters Unit Mininmum Maximum µ (day−1) 0.0120 0.0230 τ (day−1) 0 2.4×10−5 βI (day−1) 0.0494 0.2671 βH (day−1) 0.0020 0.0107 βD (day−1) 0.0494 0.2671 δ (day−1) 0.0476 0.5000 θ (day−1) 0.5000 1.0000 γ (day−1) 0.2000 0.5000 λ (day−1) 0.0328 0.1272 α (day−1) 0.0148 0.1050 We note that, because of the complexity of the nonlinear functionR0 : (0,+∞)9× [0,+∞) → R, performing an analytical study of its gradient is quite complicated. Thus, we decided to use a numerical approach. Furthermore, for the sake of sim- plicity, we decided to exclude form this study the possible relationships between parameters (e.g., α and λ are related). To do so, we apply the following algorithm which is described for a general function with Nparam parameters: Step 1 Let Nparam ∈ N be the number of parameters (here, Nparam = 10 for the function R0 given in (4.1)). We set Nscen and Npoints, the number of random scenarios and the number of points inside the interval of values of parameters, respectively. Step 2 For s = 1, 2, . . . , Nscen, Step 2.1 We randomly generate a set of parameters, denoted by Ps = (Ps,1, . . . , Ps,Nparam ), considering a uniform distribution in the interval of values reported in Table 2. We denote by Ps,p and P s,p the minimum and maximum value of parameter Ps,p, with p = 1, . . . , Nparam, respec- tively. Step 2.2 We compute R̄s,0 = R0(Ps). Step 2.3 For p = 1, . . . , Nparam, Step 2.3.1 For k = 1, . . . , Npoints, Step 2.3.1.1 We set Ps,p,k = Ps,p × (1− k−1 Npoints−1 ) + P s,p × ( k−1 Npoints−1 ). Step 2.3.1.2 We set P = Ps and replace its pth component Ps,p by Ps,p,k. Step 2.3.1.3 We compute R̃s,p,k,0 = R0(P ). Step 2.3.1.4 We compute the relative proportion error between R̃s,p,k,0 and R̄s,0 as Err(s, p, k) = |R̃s,p,k,0 − R̄s,0| R̄s,0 . End For Step 2.3.2 Considering the values R̃s,p,k,0, we compute C̃r(s, p), the Pearson correlation coefficient (see, e.g., [5]) between vectors( R̃s,p,1,0, . . . , R̃s,p,Npoints,0 ) and ( Ps,p,1, . . . , Ps,p,Npoints ) . End For EJDE-2020/62 EPIDEMIOLOGICAL MODEL BE-CODIS 21 End For Step 3 For p = 1, . . . , Nparam and k = 1, . . . , Npoints, we compute the mean relative error given by Err(p, k) = 1 Nscen Nscen∑ s=1 Err(s, p, k). End For For p = 1, . . . , Nparam, we compute the mean and maximum values of {Err(p, k)}Npoints k=1 . End For Step 4 For p = 1, . . . , Nparam, compute the mean correlation coefficient of the value of R0 with respect to the p-th parameter, given by Cr(p) = 1 Nscen Nscen∑ s=1 C̃r(s, p). End For Here, we have considered Nscen = 106 and Npoints = 100. Thus, the number of R0 evaluations performed with this algorithm is 109 and the computation time on a 3.6 Ghz I7 Intel Computer with 32 Gb of RAM is 573 seconds. The results obtained with the algorithm detailed above, when studying the sensi- tivity analysis of R0 with respect to each parameters of system (3.5) and considering data from the 2014-2016 West African EVD epidemic, are reported in Table 3. We observe that parameter τ has a negligible influence on R0. Additionally, βI and γ are the most sensitive parameters with a mean error greater than 20% and reaching errors up to 374% for the worst scenarios. All other parameters have a moderated impact on the basic reproduction ratio with a mean error lower than 8%, but may produce differences up to 162% for extreme variation cases. Regarding the Pearson correlation coefficient, we note that increasing the values of λ and δ should increases the value of R0. Furthermore, R0 exhibits a linear increasing dependency regard- ing βI , βH or βD that is why the corresponding Pearson correlation coefficient is 1; furthermore, in this case, computing the partial derivatives of R0 with respect to those parameters is trivial. On the other hand, increasing all other parameters should decreases the value of R0. The values of the basic reproduction ratios obtained with the algorithm de- scribed above were included in the interval [0.0957, 1.7424], with a mean value of 0.5833. Those results show that, for the considered range of parameters, there exist scenarios for which the EVD epidemic may remain endemic in the considered pop- ulation and, thus, the application of control measure should be applied to contain the disease spread. 4.2. Disease evolution between 2 countries. We now focus on the case of system (3.14)–(3.15), with Country 1 potentially sending infected people to Country 2. To study some representative numerical examples, we consider two set of pa- rameters, denoted by Set 1 and Set 2, detailed in Table 4, corresponding to basic reproduction ratios of 0.3491 and 1.3910, respectively. Furthermore, we assume 22 B. IVORRA, D. NGOM, A. M. RAMOS EJDE-2020/62 Table 3. Mean and maximum values of the relative error Err obtained with the algorithm in Section 4.1, when studying the sensitivity analysis of R0 with respect to each parameters of system (3.5) and considering data from the 2014-2016 West African EVD epidemic. We also report in the last column the value of Cr, the mean Pearson correlation coefficient of the value of R0 with respect to each parameter. Parameter Unit Mean Err Maximum Err Cr µ (day−1) 3.5 42 -0.99 τ (% day−1) 10−3 0.03 -0.99 βI (day−1person−1) 40 387 +1.00 βH (day−1person−1) 4.4 135 +1.00 βD (day−1person−1) 7.7 139 +1.00 δ (day−1) 5.8 41 +0.86 θ (day−1) 3.7 56 -0.98 γ (day−1) 21.3 132 -0.97 λ (day−1) 4.7 162 +0.96 α (day−1) 6.0 121 -0.95 that the population sizes are N1 = 2× 107 and N2 = 107 in Country 1 and Coun- try 2, respectively (corresponding to countries of approximatively the size of Sierra Leone and Guinea). The initial conditions are set to S1(0) = 0.999 (equivalent to 1.998 × 107 people in this particular case), E1(0) = 0.001 (equivalent to 20000 people in this particular case), S2(0) = 1 (equivalent to 107 people in this particular case) and all other ratios set to 0. Additionally, ε = 1/N1 to consider emigration flow from Country 1 to Country 2 only in the case that it exists at least one infected individual in Country 1. The model is discretized by considering the explicit Euler scheme with a step size of 1 hour (in previous numerical experiments, this scheme and step size have given a good ratio between precision and computational time, see [23]). The simulation is stopped after a maximum number of 3650 days; or if the change in state S from one iteration to other is lower than 10−9 for both countries; or if the proportion of contaminated people (e.g., people either in the state E, I, H or D) in each country is lower than the inverse of the population size. Taking into account those parameters and numerical methods, we perform the following four experiments: • Country 1 with Set 1 and Country 2 with Set 1 (Exp11): The proportion of contaminated people in both countries is presented in Figure 6. In this case, this proportion is decreasing in Country 1. In Country 2, the maximum ratio of contaminated people is 1.3 × 10−7 (equivalent to 2 people) and is reached after 8.9 days. The initial outbreak in Country 2 is due to the transportation of infected people from Country 1 occurring during the first 77.5 days of the simulation. The simulation stops after 102.7 days because of the low proportion of contaminated people in both countries. • Country 1 with Set 1 and Country 2 with Set 2 (Exp12): The evolution of the proportions of affected and healthy (i.e., people either in EJDE-2020/62 EPIDEMIOLOGICAL MODEL BE-CODIS 23 Table 4. Values of the parameters in Set 1 and Set 2 used in during the experiments presented in Section 4.2. The basic repro- duction ratio (R0) generated by those values is also reported. Parameter Unit Set 1 Set 2 µ (day−1) 0.0197 0.0120 τ (%day−1) 2×10−5 2.4×10−5 βI (day−1) 0.1147 0.2671 βH (day−1) 0.0046 0.0107 βD (day−1) 0.1147 0.2671 δ (day−1) 0.3643 0.0476 θ (day−1) 0.8500 0.5000 γ (day−1) 0.4100 0.2000 λ (day−1) 0.0564 0.1272 α (day−1) 0.0693 0.0148 R0 1.3910 0.3491 20 40 60 80 100 Time (days) 10-15 10-10 10-5 100 R a ti o Exp11 Country 1 Country 2 Figure 6. Evolution (starting from day 2) of the proportion (in logarithmic scale) of contaminated people in Countries 1 and 2 simulated for the experiment Exp11 in Section 4.2. the state S or R) people are depicted in Figure 7. We can see in Figure 7 that the ratio of affected people decreases in Country 1. On the opposite, in Country 2 the epidemic starts due to the movement of infected people from Country 1 during 77.5 days and, then, reaches an endemic equilibrium with 23% of affected people. The simulation stops after 1238 because of the stabilization of the numerical solutions. We note that, at the end of the simulation, the final value of (E2, . . . , R2) is numerically close to the value of Pe,2 reported in Theorem 3.7. • Country 1 with Set 2 and Country 2 with Set 1 (Exp21): The evolution of the proportions of affected and healthy people are shown in 24 B. IVORRA, D. NGOM, A. M. RAMOS EJDE-2020/62 0 200 400 600 800 1000 1200 Time (days) 10-40 10-20 100 R a ti o Country 1 Healthy Affected 0 200 400 600 800 1000 1200 Time (days) 0 0.2 0.4 0.6 0.8 1 R a ti o Country 2 Healthy Affected Figure 7. Evolution of the proportions of affected and healthy people in Countries 1 (in logarithmic scale) and 2 simulated for experiment Exp1 in Section 4.2. Figure 8. We can see that the epidemic reaches an endemic equilibrium of 23% of affected people in Country 1. For Country 2, due to the continuous movement of infected people coming from Country 1, the epidemic starts and remains endemic with an equilibrium of 10−4% of affected people in the population. The simulation stops after 1149 days due to the stabilization of the numerical solutions. We note that, as remarked in Theorem 3.7, despite the fact that the basic reproduction ratio of country 2 is lower than 1, the emigration of people from Country 1 does not allow Country 2 to approach a disease free state. Again, we note that, the final value of (E1, . . . , R1) is numerically close to the value of Pe,1 reported in Theorem 3.7. 0 200 400 600 800 1000 Time (days) 0.2 0.4 0.6 0.8 R a ti o Country 1 Healthy Affected 0 200 400 600 800 1000 Time (days) 10-8 10-6 10-4 10-2 100 R a ti o Country 2 Healthy Affected Figure 8. Evolution of the proportions of affected and healthy people in Countries 1 and 2 simulated for experiment Exp21 in Section 4.2. In country 2, the evolution is shown from day 2. • Country 1 with Set 2 and Country 2 with Set 2 (Exp22): In Figure 9, we report the proportions of affected people in both countries. Endemic states of 23.28% and 23.36% of affected people are reached in Countries 1 and 2, respectively. The epidemic in Country 2 suffers a delay, regarding Country 1, because of the time required to move infected people from Country 1 to Country 2. The simulation stop after 1436 days because EJDE-2020/62 EPIDEMIOLOGICAL MODEL BE-CODIS 25 of the stabilization of the numerical solutions. We observe that the final value of (E1, . . . , R1) and (E2, . . . , R2) are numerically close to the value of Pe,1 and Pe,2 reported in Theorem 3.7. 0 200 400 600 800 1000 1200 1400 0 0.05 0.1 0.15 0.2 0.25 Time (days) R a ti o Exp22 Country 1 Country 2 Figure 9. Evolution of the proportions of affected people in Coun- tries 1 and 2 simulated for experiment Exp22 in Section 4.2. We note that all values reported previously are obtained numerically and corre- spond to approximations of the results detailed in Theorem 3.7. Furthermore, we point out the fact that, when R0,i > 1 in country i ∈ {1, 2}, we have obtained in this country the numerical convergence of the solution to the endemic equilibrium point reported in Theorem 3.7. This seems to indicate that if R0,i > 1, although this result has not been proven theoretically and was only illustrated numerically with particular sets of parameters, this endemic equilibrium point is asymptotically stable. 5. Discussion and conclusions In this article, we performed an analysis of the equilibrium states of simplified versions of the Be-CoDiS model proposed in [23]. This model aims to study the spread of human diseases between countries. In Section 3.1, we have estimated a basic reproduction ratio (denoted by R0) of a version of the model for one country. In particular, we have obtained in Theorem 3.3 an analytical expression of R0 according to the model parameters. Additionally we have proven that if R0 ≤ 1, then the disease free equilibrium is globally and asymptotically stable which is a desirable biological situation because the epidemic will disappear. When R0 > 1, we show that the disease free equilibrium is unstable. This implies that the epidemic will not disappears in the considered population. Then, starting from this R0 expression and data from the 2014-16 West African Ebola epidemic, we have performed a sensitivity analysis of the basic reproduction ratio regarding the model parameters. We point out that due to biological reasons, one generally does not have control on parameters µ (the mortality rate) and δ 26 B. IVORRA, D. NGOM, A. M. RAMOS EJDE-2020/62 (transition from E to I). Taking into account this observation, due to the control measures applied by the authorities in order to eradicate the EVD spread (i.e., Isolation, Quarantine, Tracing and Increase of sanitary resources, see [38, 43, 15]), other model parameters can be regulated according to the technical limitations of those control measures. In particular, this sensitivity analysis seems to indicate that decreasing the time of detection of infectious people (1/γ, the inverse of the transition rate from I to H) and the contact rate with infectious people (βI) are the most efficient way to reduce the epidemic evolution. During the 2014-16 EVD epidemic, both variables were controlled, for instance, by monitoring the popu- lation in areas of EVD risk with health-care workers, by performing information campaigns about the disease and by isolating suspicious cases [10, 7, 42]. For ex- ample, considering the case of Guinea, it was estimated that βI and γ were reduced by the control measures, from 0.1944 and 0.2000, in December 2013, to 0.0871 and 0.3333, in October 2015, respectively [23]. Additionally, controlling contact with hospitalized people (βH) and dead body (βD), should have an impact on the EVD magnitude, although lower than reducing βI and 1/γ. In particular, it was observed that, during the first months of this EVD epidemic, around 20% of the infections were due to contacts with dead bodies [42, 39]. Additionally, the reported number of health workers infected due to contacts with hospitalized people was around 815 in May 2015, which correspond to 4% of the total number of EVD cases [40]. For these variables, control measures, such as the increase of sanitary conditions in hos- pitals and the supervision of funerals, allowed to reduce those risk factors. It was estimated that, those contact rates were both reduced by two from the beginning to the end of the epidemic [23]. The increase of sanitary resources in hospitals also allowed to increase the value of α (transition from H to R), for instance, in Guinea from 0.0847 to 0.1250 [23]. Regarding θ and λ, both parameters were controlled by reducing the duration of the funerals and the death rate (e.g., by improving the healthcare system). In particular, for Guinea, θ passes from 0.5 to 1 and λ from 0.2381 to 0.1707 [23]. We note that the classification of the importance of the model parameters in EVD control proposed here is consistent with the response plan proposed by the international community to fight the EDV outbreaks [41]. All those results seem to validate the interest of using system (3.5) and its R0 value to identify the most important factors of an epidemic evolution. Next, in Section 3.2, we have described the behavior of the epidemic evolution when two countries are connected by an emigration flow. From Theorem 3.7, we conclude that when R0,1 ≤ 1 (where R0,1 is computed with formula (3.16)) in Country 1, the evolution of the disease in Country 2 only depends on the value of R0,2. More precisely, if R0,2 ≤ 1, the epidemic disappears in Country 2, whereas if R0,2 > 1 it may become endemic in Country 2. On the opposite, when R0,1 > 1, under some reasonable assumptions, the epidemic may remain active in Country 2, even if R0,2 ≤ 1. This behavior was illustrated in Section 4.2 by performing four particular numerical experiments with several sets of parameters estimated from the 2014-16 EVD epidemic. The numerical results shown here are consistent with those found theoretically. Additionally, those numerical results seem to indicate that if R0,i > 1 for Country i ∈ {1, 2}, the epidemic in this country should converge to the endemic equilibrium point Pe,i defined in Theorem 3.7 (property which is not proven theoretically). Those outcomes tend to show the necessity to control the emigration EJDE-2020/62 EPIDEMIOLOGICAL MODEL BE-CODIS 27 flows from countries with serious epidemic scenarios. This recommendation was also proposed in the literature for the case of the 2014-16 EVD epidemic [4]. Acknowledgments. This work was carried out with financial support from the Spanish “Ministry of Economy and Competitiveness” under projects MTM2015- 64865-P, the MOMAT UCM research Group (Ref. 910480), and the “Institut de Recherche et de Développement” through projects UMI 2019-IRD and UMMUSCO- UGB. References [1] M. Anderson,; Population biology of infectious diseases: Part 1, Nature 280 (1979), 361–367. [2] N. P. Bhatia, G. P. Szegö,; Dynamical systems: Stability theory and applications, Lecture Notes in Mathematics, No. 35, Springer Verlag, Berlin, 1967. [3] D. Bichara, A. Iggidr, G. Sallet; Global analysis of multi-strains sis, sir and msir epidemic models, Journal of Applied Mathematics and Computing 44 (2014), no. 1, 273–292. [4] I.I . Bogoch, M. I. Creatore, M. S. Cetron, J.S . Brownstein, N. Pesik, J. Miniota, T. Tam, W. Hu, A. Nicolucci, S. Ahmed, J. W. Yoon, I. Berry, S.I. Hay, A. Anema, A. J. Tatem, D. MacFadden, M. German, K. Khan; Assessment of the potential for international dissem- ination of ebola virus via commercial air travel during the 2014 west african outbreak, The Lancet 9962 (2016), 29–35. [5] S. Boslaugh; Statistics in a nutshell, 2nd edition, O’Reilly Media, Incorporated, 2012. [6] F. Brauer, C. Castillo-Chávez; Mathematical models in population biology and epidemiology, Texts in applied mathematics, Springer, 2001. [7] C. D. C.; Ebola disease, Centers for Disease Control (2014). [8] T. Chen, Z. Sun, B. Wu; Stability of multi-group models with cross-dispersal based on graph theory, Applied Mathematical Modelling 47 (2017), 745 – 754. [9] G. Chowell, F. Abdirizak, S. Lee, J. Lee, E. Jung, H. Nishiura, C. Viboud; Transmission char- acteristics of mers and sars in the healthcare setting: a comparative study, BMC Medicine 13 (2015), no. 1, 1–12. [10] G. Chowell, H. Nishiura; Transmission dynamics and control of ebola virus disease (evd): a review, BMC Medecine 12 (2014), no. 196. [11] D. L. DeAngelis, L. J. Gross; Individual-based models and approaches in ecology: Populations, communities, and ecosystems, Chapman & Hall, 1992. [12] O. Diekmann, H. Heesterbeek, T. Britton; Mathematical tools for understanding infectious disease dynamics:, Princeton Series in Theoretical and Computational Biology, Princeton University Press, 2012. [13] O. Diekmann, J. A. P. Heesterbeek, J. A. J. Metz; On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations, Journal of Mathematical Biology 28 (1990), no. 4, 365–382. [14] R. T. D. Emond, B. Evans, E. T. W. Bowen, G. Lloyd; A case of ebola virus infection, Br Med J. 2 (1977), 541–544. [15] D. Fisman, E. Khoo, A. Tuite; Early epidemic dynamics of the west african 2014 ebola outbreak: Estimates derived with a simple two-parameter model, PLOS Currents Outbreaks 1 (2014). [16] M. F. C. Gomes, A. Pastore y Piontti, L. Rossi, D. Chao, I. Longini, M. E. Halloran, A. Vespignani; Assessing the international spreading risk associated with the 2014 west african ebola outbreak, PLOS Currents Outbreaks 1 (2014). [17] Y. He, S. Gao, D. Xie; An sir epidemic model with time-varying pulse control schemes and saturated infectious force, Applied Mathematical Modelling 37 (2013), no. 16, 8131 – 8140. [18] J. A. P. Heesterbeek, K. Dietz; The concept of ro in epidemic theory, Statistica Neerlandica 50 (1996), no. 1, 89–110. [19] N. Hernandez-Ceron, Z. Feng, C. Castillo-Chavez, Discrete epidemic models with arbitrary stage distributions and applications to disease control, Bulletin of Mathematical Biology 75 (2013), no. 10, 1716–1746 (English). [20] H. W. Hethcote; The mathematics of infectious diseases, SIAM Review 42 (2000), no. 4, 599–653. 28 B. IVORRA, D. NGOM, A. M. RAMOS EJDE-2020/62 [21] B. Ivorra, M. R. Ferrández, M. Vela-Pérez, A. M. Ramos; Mathematical modeling of the spread of the coronavirus disease 2019 (covid-19) taking into account the undetected infections. the case of china, Communications in Nonlinear Science and Numerical Simulation 88 (2020), 105303. [22] B. Ivorra, B. Mart́ınez-López, J. M. Sánchez-Vizcáıno, A. M. Ramos; Mathematical formula- tion and validation of the Be-FAST model for classical swine fever virus spread between and within farms, Annals of Operations Research 219 (2014), no. 1, 25–47 (English). [23] B. Ivorra, D. Ngom, A. M. Ramos; Be-codis: A mathematical model to predict the risk of human diseases spread between countries—validation and application to the 2014–2015 ebola virus disease epidemic, Bulletin of Mathematical Biology 77 (2015), no. 9, 1668–1704. [24] B. Ivorra, A. M. Ramos; Application of the be-codis mathematical model to forecast the international spread of the 2019-20 wuhan coronavirus outbreak, Technical report from Re- sarchGate. (2020), 1–13. [25] E. Kokomo, Y. Emvudu; Mathematical analysis and numerical simulation of an age- structured model of cholera with vaccination and demographic movements, Nonlinear Anal- ysis: Real World Applications 45 (2019), 142 –156. [26] J. Legrand, R. F. Grais, P. Y. Boelle, A. J. Valleron, A. Flahault; Understanding the dynamics of ebola epidemics., Med. Hypotheses. 135 (2007), no. 4, 610–621. [27] P. E. Lekone, B. F. Finkenstädt; Statistical inference in a stochastic epidemic seir model with control intervention: Ebola as a case study, Biometrics 62 (2006), no. 4, 1170–1177. [28] B. Mart́ınez-López, B. Ivorra, E. Fernández-Carrión, A. M. Perez, A. Medel-Herrero, F. Sánchez-Vizcáıno, C. Gortázar, A.M. Ramos, J. M. Sánchez-Vizcáıno; A multi-analysis approach for space-time and economic evaluation of risks related with livestock diseases: The example of fmd in Peru, Preventive Veterinary Medicine 114 (2014), no. 1, 47–63. [29] B. Mart́ınez-López, B. Ivorra, D. Ngom, A. M. Ramos, J. M. Sánchez-Vizcáıno; A novel spatial and stochastic model to evaluate the within and between farm transmission of classical swine fever virus. II. Validation of the model, Veterinary Microbiology 155 (2012), no. 1, 21 – 32. [30] B. Mart́ınez-López, B. Ivorra, A. M. Ramos, E. Fernández-Carrión, T. Alexandrov, J. M. Sánchez-Vizcáıno; Evaluation of the risk of classical swine fever spread from backyard pigs to other domestic pigs by using the spatial stochastic disease spread model Be-FAST: The example of Bulgaria, Veterinary Microbiology 165 (2013), no. 1–2, 79 – 85. [31] B. Mart́ınez-López, B. Ivorra, A. M. Ramos, J. M. Sánchez-Vizcáıno; A novel spatial and stochastic model to evaluate the within- and between-farm transmission of classical swine fever virus. I. General concepts and description of the model, Veterinary Microbiology 147 (2011), no. 3–4, 300 – 309. [32] M. I. Meltzer, C. Y. Atkins, S. Santibanez, B. Knust, B. W. Petersen, E. D. Ervin, S. T. Nichol, I. K. Damon, M. L. Washington;Estimating the future number of cases in the ebola epidemic - Liberia and Sierra Leone, 2014–2015, Centers for Disease Control. MMWR / Early Release 63 (2014). [33] C. J. Peters, J. W. Peters; An introduction to ebola: The virus and the disease, Journal of Infectious Diseases 179 (1999), no. Supplement 1, ix–xvi. [34] Z. Shuai, P. Van den Driessche; Global stability of infectious disease models using lyapunov functions, SIAM Journal on Applied Mathematics 73 (2013), no. 4, 1513–1532. [35] H. R. Thieme; Mathematics in population biology, Mathematical Biology Series, Princeton University Press, 2003. [36] P. van den Driessche, J. Watmoughm; Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Mathematical Biosciences 180 (2002), no. 1-2, 29–48. [37] P. van den Driessche, J. Watmough; Further notes on the basic reproduction number, pp. 159– 178, Springer Berlin Heidelberg, Berlin, Heidelberg, 2008. [38] W. H. O.; Fact sheet no103: Ebola virus disease, World Health Organization (2014). [39] W. H. O.; New who safe and dignified burial protocol - key to reducing ebola transmission, World Health Organization (2014). [40] W. H. O.; Health worker ebola infections in guinea, liberia and sierra leone, World Health Organization (2015). [41] W. H. O.; Who strategic response plan 2015: West africa ebola outbreak, World Health Organization (2015), ISBN: 9789241508698. EJDE-2020/62 EPIDEMIOLOGICAL MODEL BE-CODIS 29 [42] W. H. O.; Ebola publications: Case management, infection prevention and control, World Health Organization (2016). [43] W. H. O.; Global alert and response: Ebola virus disease, World Health Organization (2016). [44] W. H. O. Response Team; Ebola virus disease in west africa - the first 9 months of the epidemic and forward projections., New England Journal of Medicine 371 (2014), 1481–1495. [45] D. Yan and H. Cao,; The global dynamics for an age-structured tuberculosis transmission model with the exponential progression rate, Applied Mathematical Modelling 75 (2019), 769 – 786. [46] L. Zhu, G. Guan, Y. Li; Nonlinear dynamical analysis and control strategies of a network- based sis epidemic model with time delay, Applied Mathematical Modelling 70 (2019), 512 – 531. Benjamin Ivorra Instituto de Matemática Interdisciplinar (IMI), Department of Applied Mathematics and Mathematical Analysis, and MOMAT Research Group, Complutense University of Madrid, Plaza de Ciencias, 3, 28040, Madrid, Spain Email address: ivorra@ucm.es Diène Ngom Département de Mathématiques, Université Assane Seck de Ziguinchor, UMI 2019-IRD & UMMUSCO-UGB, BP 5233, Ziguinchor, Senegal Email address: dngom@univ-zig.sn Angel M. Ramos Instituto de Matemática Interdisciplinar (IMI), Department of Applied Mathematics and Mathematical Analysis, and MOMAT Research Group, Complutense University of Madrid, Plaza de Ciencias, 3, 28040, Madrid, Spain Email address: angel@mat.ucm.es 1. Introduction 2. Be-CoDiS model formulation 3. Analytical behavior of the Be-CoDiS model 3.1. Simplified model for one country 3.2. Simplified model for 2 countries 4. Application to the 2014-2016 West African EVD epidemics 4.1. Sensitivity analysis of the basic reproduction ratio 4.2. Disease evolution between 2 countries 5. Discussion and conclusions Acknowledgments References