ON A MODEL OF COVID-19 DYNAMICS

. A model of COVID-19 in an interconnected network of communities is studied. This model considers the dynamics of susceptible, asymp- tomatic and symptomatic individuals, deceased but not yet buried people, as well as the dynamics of the virus or pathogen at connected nodes or com- munities. People can move between communities carrying the virus to any node in the region of n communities (or patches). This model considers both virus direct (person to person) and indirect (contaminated environment to person) transmissions. Using either matrix and graph-theoretic methods and some combinatorial identities, appropriate Lyapunov functions are constructed to study global stability properties of both the disease-free and the endemic equilibrium of the corresponding system of 5 n diﬀerential equations.


1.
Introduction. Since the announcement of the novel virus SARS-CoV-2 affecting the population of Wuhan, China in December 2019, the disease caused by it, COVID-19, has extended from a small cluster in that region in China to almost every corner of the world. With globalization, countries are more interconnected than ever (with direct flights from Wuhan to the U.S. included), so human mobility is a factor that should be included in the mathematical modeling of this disease. Most current models of COVID-19 consider only one community or node (no human mobility), only direct or indirect transmission and do not include pathogen dynamics [10,13,18,26]. These and other models represent in part an extension of previous research on other epidemics including cholera, as well as Zika and Ebola virus diseases [2,6,9,12,14,17,19,21,23].
One of the main goals in this type of studies is to establish conditions under which the disease is expected to persist or die out, and to study possible control measures. To this end, it is important to study conditions for the existence and local or global stability of the steady states in the model.
The model considered here includes human mobility across n nodes (e.g. cities or countries) connected by a human mobility network (that is, people can move from node to node carrying the disease or being exposed to it), and both direct and indirect transmissions are considered simultaneously. In addition, the dynamics of concentration of the pathogen itself at different communities is also included in this model. Disease transmission can occur through person-to-person contact or through exposure to the virus in a local or visited community. The rate of exposure, contact and transmissibility vary with the parameters of the model. In practice, this is related to mitigation measures including travel restrictions, partial or complete isolation, social distancing, mask wearing, hand washing, and environment disinfection strategies. In this kind of network connectivity models, the condition that one or more of the individual basic reproduction numbers R (i) 0 , i = 1, . . . , n is greater than one, is neither sufficient nor necessary for an outbreak in the metacommunity to occur, and the disease can start even if the basic reproduction number of each community is less than one [9].
As with similar mathematical models, conditions under which the disease is expected to spread or die out are related to the stability properties of the disease-free equilibrium point, and to the value of the basic reproduction number. Here we prove that for the corresponding system of 5n differential equations, if the basic reproduction number R 0 satisfies R 0 ≤ 1, then the disease-free equilibrium point is globally stable, and otherwise it is unstable, in which case an endemic equilibrium exists and is globally asymptotically stable. These stability properties are proved using matrix and graph-theoretic methods introduced in [16,20].
2. The model. Researchers are still learning about the dynamics of COVID-19, which has been affecting regions and countries around the world, especially China, USA, Europe and Latin America. In this theoretical model, we assume that people can move with certain freedom between the n nodes in the metacommunity, possibly carrying the disease and spreading it by direct person-to-person contact or by shedding the virus in the environment on a given node or community. Indirect transmission is then possible by exposure to contaminated surfaces at any given node. As noted in [4,10,27], both asymptomatic and symptomatic carriers are able to spread the disease, although at different rates and probably at different timings of pathogen shedding. Therefore, we split infected people into asymptomatic and symptomatic variables, with distinct transmission probabilities, as well as distinct transition and shedding rates in this model. It is important to stress that in this work, asymptomatic refers to people carrying the virus but who never develop symptoms, while symptomatic refers to people also carrying the virus but showing mild to severe symptoms. It is also assumed that individuals leave their original community or node i with constant probabilities m S and m I for susceptibles and infectious (symptomatic and asymptomatic) respectively, reach their target node j with a probability Q ij and come back to node i. It is natural to assume that the graph Γ(Q) is strongly connected, meaning that people can move from any node to any other node in the network. Then, the matrix Q describing the connections between communities is stochastic, with Q ii = 0 for each i = 1, . . . , n. An explicit definition of Q following a gravity model (relative attraction of communities -related to population size and distance between communities) is given in [9].
We let S i , A i , I i , D i and V i represent the population densities of susceptibles, asymptomatic, symptomatic, disease-induced deceased but not yet buried, and virus concentration respectively at the i-th community, with i = 1, ..., n. With the virus currently moving its epicenter from the US and Europe to Latin America and probably later to Africa, it is important to consider the increased possibility of virus transmission by infected deceased people in the days before their burial (see e.g. [3,8]). With all parameters being positive, and all variables nonnegative, the corresponding system of 5n differential equations is: Qijαjp(Ej) Si where Λ i is the constant recruitment into the i-th community, and the function is the probability of direct transmission; τ and τ represent the transmissibility of asymptomatic and dead but not yet buried respectively to susceptible people; H i is a constant bound on the total population in the i-th community; α i is the rate of direct transmission in node i, g(Vi) = V i K+V i (where K is the half-saturation constant) is the probability of indirect transmission; that is, the probability of becoming infected due to exposure to the virus in node i; β i is the rate of exposure to virus-contaminated environment; µ is the natural human death rate; φ i = γ i + ξ i + µ represents the rates of recovery, disease-induced and natural death respectively of asymptomatic people in node i; similarly for φ i = γ i + ξ i + µ with respect to symptomatic people; θ represents the fraction of asymptomatic infections; p i is the rate at which virus is released by infected individuals into the environment (in the air or on surfaces) in node i, with λ and λ representing the different contributions of asymptomatic and deceased but not yet buried humans to virus concentration as compared to symptomatic ones; µ V is the rate of death of the virus (not anymore viable or capable of infection); 1/µ D is the mean time from COVID-19 death to burial. For some appropriate estimates of the values of some of these parameters, see [10,13,26].
Thus, according to model (1), susceptible individuals can get infected via direct or indirect transmission in any given community (including their own) where infected people or contaminated surfaces are present. If infected, they move to asymptomatic or symptomatic stage, from which they either recover (most of them) or die of natural or induced-disease causes. Infected individuals are able to shed the virus in their own community and in communities they visit. Denoting and using a differential inequality in Section I.6 in [11], one can show that mI Qji(λAj + Ij), we get Thus, the feasible region for model (1) is: Observe that if we write system (1) as z i = F i (z), with i = 1, . . . , 5n, then z i = 0 clearly implies F i (z) ≥ 0, and hence the nonnegative cone (z i ≥ 0, i = 1, . . . , 5n) is invariant. This means the feasible region Γ in (2) is compact and (positively) invariant. System (1) has a unique disease-free equilibrium point (DFE): The system also has an endemic equilibrium point, which will be discussed in Section 4.

3.
Global stability analysis of the DFE. In this section, we establish conditions under which the disease will persist or die out. In addition to the roles played by the problem parameters and by the density or concentration of each population (people and virus), 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 (3).
To study global stability of the DFE, we can write system (1) as a compartmental model, by splitting the variables into two compartments: a disease compartment x ∈ R 4n and a nondisease compartment y ∈ R n : In general, one denotes with each entry F i the rate secondary infections increase the i-th disease compartment, and with V i the rate disease progression, death and recovery decrease in the i-th compartment [25]. Also note that all the conditions in [20]: With this notation, model (1) can be split into disease/nondisease systems: Denote y 0 = ( Λ1 µ , · · · , Λn µ ). Then, following [20,24], we let Thus, the matrices F and V have the following form where each C i (i = 1, . . . , 6) is a positive, square matrix of order n, I represents the n × n identity matrix, and D =diag( ξ 1 , . . . , ξ n ). Thus, we have F ≥ 0, V is invertible, and also V −1 ≥ 0. The basic reproduction number R 0 of (1), defined as the expected number of secondary cases produced by an infected individual in a completely susceptible population [7,24], is the spectral radius ρ of the Next Generation Matrix (NGM): R 0 = ρ(F V −1 ). On the other hand, one can readily see that the matrix V −1 F is nonnegative and irreducible, guaranteeing that R 0 is a positive eigenvalue of V −1 F (and of F V −1 ) and that there exists a positive left eigenvector ω of V −1 F corresponding to R 0 [1]. This eigenvector ω will be used in the definition of a Lyapunov function below, which in turn is used to prove the following stability result.
Proof. Following the matrix-theoretic method in [20], we set y).
Observe that f (x, y) ≥ 0. Then, by applying Theorem 2.1 in [20], we see that is a Lyapunov function for (1) in Γ, when R 0 ≤ 1. To prove global stability, first consider R 0 < 1. One can verify that the n-dimensional disease-free system associated to (1) has a unique equilibrium point y 0 = ( Λ1 µ , · · · , Λn µ ), and that this equilibrium point is GAS in R n + . One can also see that f (0, y) = 0 in Γ. Then, applying the first part of Theorem 2.2 in [20], we conclude that if R 0 < 1, the DFE (3) is GAS in Γ. Now let R 0 = 1. Observe that Then, D = 0 implies f (x, y) = 0, as ω is strictly positive and all diagonal entries of V −1 are also strictly positive. Then, we get S i = S 0 i and V i = 0 = V 0 i , for i = 1, . . . , n -see notation in (3). One can show that the remaining 3n-dimensional system in A i , I i and D i has the unique equilibrium (A 0 i , I 0 i , D 0 i ) = (0, 0, 0) ∈ R 3n , and that this equilibrium is GAS in R 3n + . But we already have S i = S 0 i and V i = V 0 i ; thus, the largest and only invariant set where D = 0 is the DFE. Using LaSalle's invariance principle, one concludes that the DFE (3) is GAS in Γ, when R 0 = 1.
Existence of an Endemic Equilibrium. In Theorem 3.1, we proved that if R 0 ≤ 1, 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 0 > 1, the DFE loses stability (see e.g. Theorem 2 in [24], or Theorem 1 in [25]). As noted in the proof of Theorem 2.2 in [20], this instability implies uniform persistence of (1), and by Theorem D.3 in [22], the invariance of Γ and the uniform persistence imply the existence of an endemic equilibrium (EE). We will denote this endemic equilibrium point of (1) as In Theorem 4.2 below, we show that the condition R 0 > 1 also implies that the EE is globally asymptotically stable.

Global stability of endemic equilibrium.
To study global stability of the EE (6), we use a graph-theoretic method as presented in [20]. First, we summarize some basic terminology and results from [16,20] on directed graphs, including a technique for the construction of a Lyapunov function. Given a weighted digraph Γ(A) with p vertices, the p × p weight matrix A is defined with a ij > 0 equal to the weight of arc (j, i) if it exists, and a ij = 0 otherwise. The Laplacian L of Γ(A) is defined as Let c i be the cofactor of l ii . If Γ(A) is strongly connected, then c i > 0, for all i = 1, . . . , p. The following combinatorial identities are useful in finding explicit expressions for c i : If a ij > 0 and the out-degree of vertex j satisfies d + (j) = 1, for some i, j, then If a ij > 0 and the in-degree of vertex i satisfies d − (i) = 1, for some i, j, then The following theorem provides a graph-theoretic technique to construct a Lyapunov function D.
Theorem 4.1. [16,20] For a given open set E ⊂ R m , and a function f : E → R m , consider the systemż = f (z), (9) and assume that (i) There exist functions Q i : E → R, G ij : E → R, and constants a ij ≥ 0 such that where S(C) denotes the set of all arcs in C.
With these tools at hand, we can give a result on global stability of the endemic equilibrium of system (1). In the theorem below, we use the functions p(E j ) and g(V j ) as defined in (1). We also use the following notations: and h i as in Section 3.
Proof. For i = 1, . . . , n, define Then, taking derivatives, and solving for Λ i , φ i and φ i at the equilibrium equations in (1), we have With corresponding functions G ij within both sets of brackets above. Similarly, for i = 1, . . . , n, defining: we get: With the coefficients a ij > 0 found above, we can construct the digraph Γ(A) shown in Figure 1 so that along any cycle C, we get (s,r)∈ S(C) Then, we can let for some constants c i ≥ 0, i = 1, . . . , 2n. Also observe that the in-degree of node n + i satisfies d − (n + i) = 1, so that from (8) we get c n+i = n j=1 c j a j,n+i /a n+i,i .
On the other hand, note that Gi,n+j + Gn+j,j = Gij . Then, following the same approach used in [20], we get ci (aij + ai,n+j) ln One can let a ij = a ij + a i,n+j , and consider the corresponding digraph Γ ci aij ln In addition, one can show that the endemic equilibrium (S * i , A * i , I * i , D * i , V * i ) is the largest (and only) invariant set for which Q = 0. Then, using LaSalle's invariance principle we conclude that such equilibrium is globally asymptotically stable in int(Γ), and thus unique.

Final remarks.
A general theoretical model of COVID-19 dynamics has been proposed. The model considers the connection and interaction between people from n different communities. The disease can spread through direct (person-toperson) or indirect (humans exposed to contaminated environment) transmission. While most infection transmission might happen within a given community, as people move from their original community i to community j, they can either get infected through contact with (asymptomatic or symptomatic) infected people or contaminated surfaces in node j, or they can pass the disease to susceptible people and shed pathogens in the environment in node j, and go back to their node i. The possibility of people getting infected when handling disease-induced dead bodies, and the dynamics of pathogen concentration are also included in the model. A simple simulation of a 2-node system is shown above, assuming 200K and 400K people in the communities respectively. For the particular case of the U. S., after about five months of the first case reported, the number of asymptomatic people seems to be much larger of that of symptomatic, and they seem to be the group most responsible for recent spikes in the number of infections and spreading of the disease. This is probably explained in part because of the inconsistent protocols about reopening the country, with strategies varying very much from state to state and from county to county, ignoring the role that connectivity and human mobility across counties and states play in the persistence of the disease.
Using matrix and graph-theoretic techniques from [16,20], it was shown that when the basic reproduction number R 0 satisfies R 0 ≤ 1, the disease-free equilibrium point is globally stable, whereas for R 0 > 1 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 (1)  This work builds on previous research on modeling and analyzing mathematical models of COVID-19 and other disease epidemics, including those in [2,5,9,8,10,13,15,16,21,26]. Some extensions and generalizations of this model are possible for future research, such as splitting infectious people into tested (or detected) and non-tested (not detected), isolated and not isolated, adding the exposed stage into the dynamics, and possibly also splitting symptomatic into those with mild and those with severe symptoms.