A General Multipatch Model of Ebola Dynamics


 A model for the transmission dynamics of Ebola virus in a multipatch network setting is studied. The model considers the contribution to the dynamics by people who are susceptible, infectious, isolated, deceased but still infectious and not yet buried, as well as the dynamics of the pathogen at interacting nodes or patches. Humans can move between patches carrying the disease to any patch in a region of n communities (patches). Both direct and indirect transmission are accounted for in this model. Matrix and graph-theoretic methods and some combinatorial identities are used to construct appropriate Lyapunov functions to establish global stability results for both the disease-free and the endemic equilibrium of the model. While the model is focused on Ebola, it can be adapted to the study of other disease epidemics, including COVID-19, currently affecting all countries in the world.


Introduction
The Ebola Virus Disease (EVD) is a rare but deadly disease that has mostly a ected countries in West Africa. It is caused by an infection with a group of viruses within the genus Ebolavirus, with the most common one being the Zaire ebolavirus (simply known as Ebola virus). According to the CDC -Center for Disease Control [6], Ebola virus was rst discovered in 1976 in the Democratic Republic of Congo (DRC), and since then it has caused outbreaks in several African countries. The 2014-2016 Ebola outbreak is the largest and started in Guinea. From there, it quickly spread to neighboring countries, especially Liberia and Sierra Leone, a ecting nearly 30,000 people and causing over 11,000 deaths. The more recent outbreak in DRC in 2018 reached about 3,000 cases, causing about 2,000 deaths. The numbers of deaths and countries a ected could probably have been much higher if the virus had started in countries that are fully connected to the world due to globalization. The corona virus disease (COVID- 19), currently causing an unprecedented pandemic, illustrates the importance of understanding the dynamics of a disease in an interconnected world (through a multipatch model). The model studied here can be adapted to the dynamics of COVID-19 and other disease epidemics.
Ebola transmission can be direct (person-to-person) or indirect (person-environment-person). One must remark that infected people who have been isolated can also transmit the disease. Due to funeral practices, deceased people can also transmit the disease before they are buried. Infected people shed the pathogen in the environment by leaving behind clothes, bedding or other objects contaminated with their body uids.
Most of the current mathematical models to study Ebola dynamics consider some combination of the following cases: only one community or a network of communities (the vast majority of them consider only one community); human mobility in a 2 or a multi-patch model; direct or indirect transmission (the majority addressing only direct transmission); inclusion (or omission) of the dynamics of pathogen concentration at each community [2,5,8,10,13,18,20].
Mathematical models of various disease epidemics typically focus on estimating the basic reproduction number, as a rst indicator or predictor of the conditions under which the disease will be endemic or die out. Then, the goal is to establish local and global stability properties of the equilibrium points (steady states), and giving an insight on possible measures to be taken to prevent or control the epidemics. See, e.g. [1,3,7,14,16,19,22]. For theoretical and numerical analyses related to chemotaxis models from mathematical biology described by classical PDEs, we refer the reader to [15,26].
One very important factor to consider in mathematical models of EVD and related diseases is the persistence of human mobility through di erent communities (patches or nodes) in the metacommunity. Local authorities are typically forced to impose lock down and even curfew measures to dramatically slow down human mobility and therefore the transmission of the disease (this is currently happening with the COVID-19 pandemic). Human mobility increases both direct (person-to-person) and indirect (exposure to contaminated environment) transmission, and it could dramatically increase the rate of spread to the point of collapsing existing health care infrastructure.
The mathematical model in this work considers a number n of communities (nodes or patches) connected by a human mobility network, to account for the fact that people move from node to node unintentionally carrying the disease or being exposed to it. Due to this mobility and connectivity in the network, both direct and indirect transmission happen simultaneously and the spread of the disease crosses city and country borders. The more connected a country is to the rest of the world (mostly because of globalization), the more likely the disease will spread across countries, close or far from the country where the diseases rst started.
In this kind of network connectivity models, the condition that one or more of the individual basic reproduction numbers R i , i = , . . . , n is greater than one, is not su cient nor necessary for an outbreak in the metacommunity to occur, and the disease can start and spread even if the basic reproduction number of each community is less than one [11].
In this article, we show that for the corresponding system of n di erential equations, if the general basic reproduction number R satis es R ≤ , then the disease-free equilibrium point is globally stable, and otherwise it loses stability, in which case an endemic equilibrium exists and is globally stable. These stability properties are proved using matrix and graph-theoretic methods and some combinatorial identities introduced in [17,21].

The Model
One of the main assumptions in this model is that people can move between the n nodes in the metacommunity, possibly carrying the disease and spreading it by direct person-to-person contact or by shedding the pathogen in the environment. This assumption is even more relevant in a pandemic case, as currently seen with COVID-19. While direct transmission is the main concern, indirect transmission is also a very important factor and it can happen by direct exposure to objects and the environment which have been contaminated by people living in the community or by travelers coming from any other community in the region. It is then important to consider the rate of change of pathogen concentration in the environment in each community.
In this model, we consider three groups of humans who can transmit the disease: Infected who are mostly freely moving around communities, infected people who have been isolated into some health facility or their own homes, and dead people who are still infectious but have not yet been buried due to funeral practices in some regions in Africa. For these groups, we take care of assigning di erent (direct) transmissibility rates as well as di erent rates at which they release the pathogen in the environment.
Following an approach by [11], we assume that individuals leave their original community or node i with probability m S for susceptibles and m I for infectives, reach their target node j with a probability Q ij and come back to node i. Clearly, ≤ m S , m I , Q ij ≤ . It is natural to assume that the graph Γ(Q) is strongly connected, meaning that people can potentially move from any node to any other node in the network. Clearly, the matrix Q describing the connections between communities is stochastic, with Q ii = for each i = , . . . , n. A possible explicit de nition of Q following a gravity model (relative attraction of communities -related to population size and distance between communities) is given in [11].
We let S i , I i , J i , D i , B i represent the population densities of susceptibles, infectious, isolated, deceased but still infectious and not yet buried and pathogen concentration respectively at the i-th community, with i = , ..., n. Assuming all parameters to be positive, and all variables nonnegative, the corresponding system of n di erential equations is: where Λ i is the constant recruitment into the i-th community, the function p(C i ) is given by p( and it represents the probability of direct transmission; τ and τ represent the transmissibility of isolated and dead but not buried respectively to susceptible humans; H i is a constant bound on the total population in the i-th community (see below); α i is the rate of direct transmission in node i, is the probability of indirect transmission due to exposure to pathogen B i in node i; β i is the rate of exposure to contaminated environment; µ is the natural human death rate; ϕ i = γ i + ξ i + µ represent the rates of recovery, disease-induced and natural death respectively of infectious humans in node i; similarly for ϕ i = γ i + ξ i + µ with respect to isolated humans; θ represents the fraction of infectious humans who are not isolated; r i is the rate at which pathogens are released by infectious individuals in the environment, in node i, with λ and λ representing the di erent contributions between infected humans; µ B is the decay rate of Ebola virus in node i. For some appropriate values of these parameters, see [2,11].
Thus, according to model (2.1), susceptible individuals can get infected via direct or indirect transmission in any given community (including their own) where people are infected or the environment is contaminated. If infected, susceptible people move to infectious and not isolated or infectious but isolated stages, from which they either recover or die of natural or induced-disease causes. Infected individuals are able to shed pathogens in their own community or the one they visit.
and using a di erential inequality in [12], . Thus, the feasible region for model (2.1) is: Observe that if we write system (2.1) as z ′ i = F i (z), with i = , . . . , n, then z i = clearly implies F i (z) ≥ , and hence the nonnegative cone (z i ≥ , i = , . . . , n) is invariant. This implies that the feasible region Γ in (2.2) is compact and (positively) invariant.
System (2.1) has a unique disease-free equilibrium point (DFE): The system also has an endemic equilibrium point, which will be discussed in Section 4.

Global Stability Analysis of the DFE
In this section we establish conditions for the epidemics to persist or die out. In addition to the role played by the problem parameters and by the density of each population (people and pathogen), the mobility of infected people throughout the network also plays a key role. We start by studying the stability properties of the disease-free equilibrium (2.3).
For global stability analysis of the DFE, system (2.1) is written as a compartmental model; that is, we split the variables into two compartments: a disease compartment x ∈ R n and a nondisease compartment y ∈ R n : In general, one denotes by F i the rate secondary infections increase the i-th disease compartment, and by V i the rate disease progression, death and recovery decrease the i-th compartment [25]. Also note that all the conditions in [21] and [25]: y) ≤ whenever x i = , and V i (x, y) ≥ , for i = , ..., n hold.
With this notation, model (2.1) can be split into disease/nondisease systems: where the nondisease system y ′ = G(x, y) is then de ned by the rst n di erential equations S ′ i in (2.1).
The basic reproduction number R of (2.1), de ned as the expected number of secondary cases produced by an infected individual in a completely susceptible population [9,24], is the spectral radius ρ of the Next Generation Matrix (NGM): R = ρ(FV − ). On the other hand, one can readily see that the matrix V − F is nonnegative and irreducible, guaranteeing that R is a positive eigenvalue of V − F (and of FV − ) and that there exists a positive left eigenvector ω of V − F corresponding to R [4]. This eigenvector ω will be used in the de nition of a Lyapunov function below, which in turn is used to prove the following stability result.
Observe that we have f (x, y) ≥ and that f ( , y) = . Then, by applying Theorem 2.1 in [21], we see that To prove global stability, rst consider R < . One can verify that the n-dimensional disease-free system associated to (2.1) has the unique equilibrium point y = ( Λ µ , · · · , Λn µ ), and that this equilibrium point is GAS in R n + . One can also see that f ( , y ) = in Γ. Then, applying the rst part of Theorem 2.2 in [21], we conclude that if R < , then the DFE (2.3) is GAS in Γ.
Then, Q ′ = implies f (x, y) = , as ω is strictly positive and all diagonal entries of V − are also strictly positive. Then, we get S i = S i and B i = = B i , for i = , . . . , n -see notation in (2.3). One can show that the remaining n-dimensional system in I i , J i and D i has the unique equilibrium (I i , J i , D i ) = ( , , ) ∈ R n , and that this equilibrium is GAS in R n + . But we already have S i = S i and B i = B i ; thus, the largest and only invariant set where Q ′ = is the DFE. Using LaSalle's invariance principle, one concludes that the DFE (2.3) is GAS in Γ, when R = .
Existence of an Endemic Equilibrium. In Theorem 1, we proved that if R ≤ , then the disease-free equilibrium point (DFE) is globally asymptotically stable in the given feasible region Γ. When this condition is reversed; that is, when we instead have R > , the DFE loses stability (see, e.g. Theorem 1 in [25]). As noted in the proof of Theorem 2.2 in [21], this instability implies uniform persistence of (2.1), and by Theorem D.3 in [23], the invariance of Γ and the uniform persistence imply the existence of an endemic equilibrium (EE). We will denote this endemic equilibrium point of (2.1) as In Theorem 3 below, we show that the condition R > also implies that the EE is globally asymptotically stable.

Global stability of endemic equilibrium
To study global stability of the EE (3.3), we use a graph-theoretic method as presented in [21]. First, we summarize some basic terminology and results from [17,21] on directed graphs, including a technique for the construction of a Lyapunov function.
A pair (i, j) is called an arc from vertex i to vertex j. Given a weighted digraph Γ(A) with p vertices, the p × p weight matrix A is de ned with a ij > equal to the weight of arc (j, i) if it exists, and a ij = otherwise. The Laplacian L of Γ(A) is de ned as The following theorem provides a graph-theoretic technique to construct a Lyapunov function D.
Theorem 2. [17,21] For a given open set E ⊂ R m , and a function f : E → R m , consider the systeṁ

3)
and assume that (i) There exist functions D i : E → R, G ij : E → R, and constants a ij ≥ such that With these tools at hand, we can give a result on global stability of the endemic equilibrium of system (2.1).
In the theorem below, we use the functions p(E j ) and g(B j ) as de ned in (2.1). We also use the following notations: and h i as in Section 3.
Proof. For i = , . . . , n, de ne Then, taking derivatives, and solving for Λ i , ϕ i and ϕ i at the equilibrium equations in (2.1), we have

Final Remarks
In this work, we have introduced a general multipatch model of Ebola dynamics. The model considers a metacommunity of n nodes (cities, countries, etc.) interconnected by a human mobility network. The EVD can spread through direct (person-to-person) or indirect (humans exposed to contaminated environment) transmission. As humans move from their original community or node i to community or node j, they can either get infected through contact with infected people (included those isolated and dead but not yet buried) or with contaminated environment in node j, or they can infect susceptible people by shedding pathogens in community j, and go back to their node i.
Using matrix and graph-theoretic techniques from [17,21], it was shown that when the basic reproduction number R satis es R ≤ , the disease-free equilibrium point is globally stable, whereas for R > it loses stability and an endemic equilibrium exists; it was then proved that such endemic equilibrium is unique and globally stable in the interior of a feasible region Γ.
Model (2.1) is general in the sense that it includes the following cases: If all α i = , then only indirect transmission is possible; if all β i = , then only direct transmission is possible; if θ = , λ = τ = and all ϕ i = , then isolated but still infectious people are not considered in the model; if λ = τ = and all ξ i = ξ i = , then dead individuals who are still infectious but not yet buried are not considered; if m S = m I = , then transmission is only possible within a given patch; nally, if n = , then the model reduces to one isolated community (and if n = , then it is simply a two-patch model).
Some extensions and generalizations of this model are possible for future research, including but not limited to: more general incidence functions for direct and indirect transmission, studying possible control strategies, and including exposed, recovered or other individuals as a separate groups. A promising new avenue of research is adapting these and other techniques to the study of the dynamics of the most recent corona virus disease COVID-19 that has caused a pandemic a ecting the whole world, with over 30 million cases and about 555,000 deaths in the US alone, as of the writing of this article. Mathematical models can be very helpful tools in understanding the dynamics of disease epidemics and providing with a rigorous analysis of those factors which are the most crucial into preventing a disease to spread beyond control.