key: cord-0166691-dasa87e5
authors: Basnarkov, Lasko
title: Epidemic spreading model of COVID-19
date: 2020-05-24
journal: nan
DOI: nan
sha: 38a574ba63da77d503b1b652f3430c5d96cbcf2e
doc_id: 166691
cord_uid: dasa87e5

We study Susceptible-Exposed-Asymptomatic-Infectious-Recovered (SEAIR) epidemic spreading model inspired by two characteristics of the infectiousness of COVID-19: delayed start and its appearance before onset of symptoms, or even with total absence of them. The model is theoretically analyzed in continuous-time compartmental version and discrete-time version on random regular graphs and complex networks. We show analytically that there are relationships between the epidemic thresholds and the equations for the susceptible populations at the endemic equilibrium in all three versions, which hold when the epidemic is weak. We provide theoretical arguments that eigenvector centrality of a node approximately determines its risk to become infected.

Understanding epidemic spreading of contagious diseases and effectiveness of various countermeasures is of high interest for public health and the society in general, with important contributions provided by epidemiologists, mathematicians and physicists as well. Although earliest theoretical work in mathematical epidemiology dates back to Daniel Bernoulli [1] , the development of the modern approach started in the beginning of the past century [2] [3] [4] . In the last two decades, since the emergence of the complex networks theory, epidemic modeling has gained novel insights. By modeling the contacts between the individuals with complex networks, association was found between the epidemic threshold with the network properties like the degree distribution or eigenvalue of the adjacency matrix [5] [6] [7] [8] [9] [10] . Moreover, the epidemic spreading has grown as a concept that extends its original design for modeling diffusion of infectious disease to sharing ideas, rumors, or computer viruses [9] .

In the classical approach, the individuals are conveniently grouped in compartments or classes. The mathematical models for epidemic spreading are systems of differential equations for the evolution of the size of those compartments. In such models, there is assumption of homogeneous mixing, which means that the pathogen can spread between each pair of individuals with equal probability. In a more realistic modeling each individual is considered as node in certain network of contacts, where infection can spread only among neighbors in that network. This is particularly relevant, because real world networks besides being random, among others they possess properties like small-world phenomenon [11] , or scale-free distribution of the node degrees [12] . One approach for studying the disease spreading on networks is the heterogeneous mean-field [6, 13] in which all nodes with the * lasko.basnarkov@finki.ukim.mk same degree are assumed to be statistically equivalent. The quenched mean-field technique is applied in even more realistic scenario, where each node is treated separately [14, 15] . Among major contributions in the field of epidemic spreading on networks are also the studies of spreading on multilayer networks [16] , disease localization [17] and the assumption of non-exponential distribution of periods between consecutive events [18, 19] .

Various infections are characterized with different stages in the course of development of the disease in the host, starting from contracting the pathogen to the healing. Depending on the disease under study, several compartments, or classes are defined in order to differentiate between the stages. The most frequently used compartments are the Susceptible (S), Exposed (E), Infectious (I) and Recovered (R) [9, 20] . The meaning of the compartmental symbols usually is: S -healthy individuals subjected to infection, E -infected which do not transmit the disease yet, I -infected and infectious and R -cured which cannot become infected again. The most popular models are the SIS and SIR, which are sufficiently simple to provide mathematical tractability, and powerful enough to capture the features of epidemics of many contagious diseases [20] .

The ongoing COVID-19 is the largest pandemic in modern history. The understanding of the virus and its influence on the infected individuals are still in progress. However, it was uncovered that it features a slow progress that might take a long period from inception to cure and some infected individuals can be unnoticed because they might not have any symptoms while spreading the disease [21] . Some studies even report that the infectiousness is strongest before and at the onset of appearance of symptoms [22] . In order to properly describe such contagion one needs a model which incorporates delay effects and also has different stages with possibly different virulence. Since its emergence, numerous mathematical models of the COVID-19 with different complexity have appeared (a short and not exhaustive list of references is [23] [24] [25] [26] ). In this theoretical work we con-sider the known Susceptible-Exposed-Infected-Recovered (SEIR) model augmented with another state, the asymptomatic state (A), which precedes the infectious. The chosen SEAIR model has a built-in delay which means that the infected person do not start to spread the disease immediately, allows to address the different contagiousness of the disease in different phases and captures the presence of undetected spreaders. Also, in our opinion it is simple enough to allow for theoretical study of different properties. We study the SEAIR model with two approaches: the classical, which uses differential equations, and the statistical physics-based, by considering discretetime evolution model on complex networks. We study the latter model on random regular graph and on complex network separately. For the three versions we obtain the epidemic threshold, equation for the fraction of susceptible individuals at the end of epidemic and study linear stability of the disease-free and endemic equilibria. The results for the random regular graph hold when the contagiousness is weak, while for the complex network it is also needed that epidemic is weak in a sense that small fraction of population is affected by the epidemic. We furthermore study the roles of the leading eigenvalue and the principal eigenvector of the adjacency matrix in the epidemic spreading. It is already known that leading eigenvalue determines the epidemic threshold and thus determines the stability [14, 27] . We show that the principal eigenvector and its associated eigenvector centrality have important role in estimation of the risk of infection. Finally, the techniques which we apply for analysis of the Jacobian matrix of the evolution might be relevant in studies of linear stability of coupled multidimensional dynamical systems.

The paper is organized as follows. In Section II we introduce and analyze the SEAIR compartmental model. In the following Section III is studied the discrete-time version, while in Section IV we present some results of numerical experiments. The paper ends with the Conclusion.

Let the variables S, E, A, I and R denote the fraction of individuals which are respectively susceptible, exposed, asymptomatic, infectious and recovered. We assume that the exposed state corresponds to the incubation period when host has the pathogen, but he or she cannot infect the others. In the asymptomatic state the individual spreads the disease, possibly with higher virulence without being aware about having the virus. The person can even recover without ever noticing that he, or she had the disease. Certain percentage of carriers of the virus will show symptoms and we classify them as infectious. Let the infecting rate of the asymptomatic persons is α, while of infectious ones is β. The rate at which exposed individuals become asymptomatic is γ. The growth of the fraction of infectious hosts, which have symptoms is determined with rate σ. For simplicity, we assume that healing of both the asymptomatic and infectious persons is modeled with same rate µ. We emphasize that µ does not exactly correspond to the time of complete healing, but the period in which a person can infect the others. With these assumptions one has the following SEAIR compartmental model

We have neglected the births and deaths in the population and one can easily verify that the total number of persons in all states is constant S(t)+E(t)+A(t)+I(t)+ R(t) = 1.

One trivial solution of the system (1) is the diseasefree state S = 1, when the pathogen is absent. If some virus is introduced, an epidemic can occur. Then there is an endemic equilibrium which corresponds to the situation when the fraction of susceptibles is not sufficient for further spread of the disease. When epidemic occurs, the number of unaffected people can be obtained by standard technique [28] . If we sum the top four equations in (1) the following relationship will hold

In situations when the epidemic starts with negligibly small number of virus bearers, by taking that at the finish of the epidemic the fractions of individuals with the pathogen is zero, after integration of the last equation, one obtains

The first equation in (1) can be rewritten as

which by integration will result in another relationship between the initial and final fractions of susceptibles

One can also integrate the fourth equation in (1) on both sides to obtain

The last result provides a relationship between asymptomatic and infectious fractions in the course of the whole epidemic

By combining the relationships (3), (5) and (7), the following equation for the fraction of unaffected individuals is obtained

Using the fact that f (x) = ln x is steeper than g(x) = x for x < 1, one can verify that

This implies that the transcendental equation (8) has a solution only if

The last inequality is the condition of existence of endemic equilibrium S = S(∞); E = A = I = 0; R = S(0) − S(∞) of the system (1).

To study the linear stability of the equilibrium states, one should linearize the system (1). The respective Jacobian matrix J = [∂Ḃ/∂C]; B, C ∈ {S, E, A, I, R}, reads

At the disease-free state for which S = 1 and E = A = I = R = 0, the Jacobian has rather simple form

Because the first and last columns are zero, this Jacobian has two trivial eigenvalues zero, while the other three are the roots of the polynomial

We show in Appendix A that the three eigenvalues of the Jacobian have negative real part, if the following relationship holds

The obtained inequality is opposite of the condition of the existence of endemic equilibrium, as one can expect. If (13) holds, (10) does not, the disease-free state is stable and epidemic will not occur. In the opposite, if (13) is not satisfied, the equilibrium (1, 0, 0, 0, 0) is unstable, the epidemic will ensue, and the size of unaffected population can be obtained from (8) . Thus the threshold at which epidemic can emerge is the following relationship

The linear stability analysis of the endemic equilibrium can be applied by the same procedure as for the diseasefree one. In this case in the Jacobian (11) one should take S = S(∞) and E = A = I = 0. By applying the same procedure which is explained in Appendix A, it will be obtained that the only difference from the disease-free case is that instead of α and β it should be used S(∞)α and S(∞)β respectively. This would simply change the condition for stability of the endemic equilibrium to

From (8) one has the following relationship for the fraction containing the parameters

Plugging the last relationship in (15) , it will be obtained that the endemic equilibrium is stable once the following holds

The last inequality can be rearranged as

which holds always since S(0) > S(∞).

Subsequent investigation of a disease spreading model when one needs to account for the contacts between the individuals, is to use epidemic spreading on complex networks framework. We consider discrete-time evolution version of the proposed SEAIR model with finite population of N individuals. The network of contacts is conveniently modeled with fixed undirected graph, in which the vertices are the individuals, while the links exist between those persons which have contact to each other. This means that the disease can be transmitted only between neighbors in the graph. The dynamical variables of interest in this model are the set of probabilities for each node being in certain state, S, E, A, I, or R at given moment n. We will also denote the parameters with the same Greek letters as in the compartmental model. They correspond to the same transitions and have a meaning of probabilities instead of rates. To be more precise, α is the probability that asymptomatic person will infect a susceptible neighbor at one time step, while β is the respective probability in the case of contact between infectious and susceptible individual. Once becoming exposed, the person can proceed into asymptomatic phase with probability γ at one time step, or remain in the same state with probability 1 − γ. The respective probability to show symptoms by asymptomatic individual is σ. Again, as in the compartmental model, we assume identical probability µ to become recovered, for both the asymptomatic and infectious state.

The dynamics of the discrete-time version of the SEAIR model is build upon the model considered in [14] . Denote the probabilities that the individual i at the discrete moment n is in respective state with p S,i (n), p E,i (n), p A,i (n), p I,i (n) and p R,i (n). Our model assumes a reactive process of epidemic spreading, which means that in each time step every individual has a contact with every neighbor [29] . Then, certain susceptible person i at the moment n + 1 will not receive the infection from any of its neighbors with probability [14, 29] 

where N i denotes the set of neighbors of i. The individual will remain susceptible at moment n + 1, if he, or she, did not receive the contagion, which means that

Otherwise, an individual can become exposed if he, or she has been susceptible before and received the contagion, or continue to be exposed at the next moment if the incubation has not finished, with probability

The probability of being asymptomatic at the next moment is

where the last term accounts for the situation that neither the symptoms will appear, nor healing will happen in one time step. The node i will be in state I at moment n + 1 with probability

where the former term describes the probability to show symptoms, if in the previous moment the node was asymptomatic, while the last one corresponds to recovering. Finally, the probability of being recovered at some moment is

The set of equations (20) to (24) determines discrete-time dynamical system of equations for evolution of probabilities of the states for each node in the network. This can be solved numerically for arbitrary initial condition and determine the progress of the epidemic at each moment. In practice one can make such studies with networks with size depending on the computational capacities at hand.

A. Epidemic spreading on random regular graphs

We will pursue our analysis of spreading processes on random regular graphs where each node has the same degree k. For infinitely large random regular graphs, the probabilities of the states are equal for all nodes and one can drop the index of the node. The probability to avoid infection (19) will be simplified to

Then the system of equations for discrete-time epidemic spreading on random regular graph reads

Here we have also two equilibrium points: one where all individuals are susceptible p = (1, 0, 0, 0, 0) and the other when the fraction of susceptible individuals is such that it prevents further spread of the disease p = (p * S , 0, 0, 0, 1 − p * S ). In the Appendix B we show how application of similar reasoning as for the compartmental model can allow to find a closed form equation for the number of susceptible individuals at the end of epidemic. It can be applied when the contagion probabilities α and β are very small. As is shown in the Appendix B, the equation for determination of the probability of susceptible state at the end of epidemic is very similar to the respective one for the compartmental case

The last result extends the one for all-to-all coupling considered in compartmental models, where effectively each individual can get the disease from any one in the population. Here, it is obtained that appropriately modified relationship holds for restricted number of contacts. By repeating the same analysis as for the compartmental model, one can also obtain that the condition for existence of endemic equilibrium is

We can further make linear stability analysis of equilibria by linearizing the evolution equations. The respective Jacobian matrix at the fixed points for which

Solving the characteristic equation Q(λ) = det(J − λI) = 0 of the Jacobian (29), for p * S = 1 will result in two trivial eigenvalues equal to one and other three. In vicinity of the epidemic start, the nontrivial eigenvalues λ are the roots of the polynomial in λ,

One can introduce new variable ν = λ − 1 and obtain

The last characteristic function is very similar to that which corresponds to the stability of the compartmental model (13) . The only difference is that here one should take kα and kβ instead of α and β respectively, and that the independent variable is ν. Thus, one can apply the same reasoning as in that case, which is explained in Appendix A. The conclusion would be that all roots ν would have negative real part only if

We remind that because ν = λ − 1, when ν has negative real part, λ < 1, and thus the same condition determines when the disease-free state is linearly stable. The linear stability of the endemic equilibrium is established from the leading eigenvalue of the same Jacobian matrix (29) as the disease-free one, but using p * S = p S (∞) obtained from (27) . Again, the procedure is nearly the same as for the compartmental version and the only difference is that instead of k one should use kp S (∞) in all analysis. Then the endemic equilibrium will be stable, if the condition similar to (32) is satisfied

Without showing the details, we will just mention that once the endemic equilibrium in this discrete-time disease spreading model exists, it is linearly stable.

Let us now consider the general case when the contacts between individuals are described with complex network. In studies of interacting units coupled in a network it is typical to define the state by stacking the state vectors of each node one on top of another. In this case another ordering is more appropriate [27] . First, create vector of the probabilities of susceptible states of all nodes p S = [p S,1 , p S,2 , . . . , p S,N ] T , then those of the exposed states p E = [p E,1 , p E,2 , . . . , p E,N ] T , and likewise for the remaining three p A , p I and p R . Also, denote with A the adjacency matrix of the network of contacts between the individuals with elements A i,j = 1 only if nodes i and j are neighbors and A i,j = 0 otherwise. Under general circumstances, determination of the probabilities at the end of epidemic in case it happens, is very complicated, if not impossible. However, when the contagiousness is weak, which means that α ≪ 1 and β ≪ 1, one can obtain similar expressions which relate initial and final probabilities of susceptible state as for the former two models. As it is explained in details in the Appendix C, when α ≪ 1 and β ≪ 1, the susceptibility probability vector can be calculated from the following self-consistent system

The solution of the system of transcendental equations equations (34) consists of the susceptibility vector at the endemic equilibrium p S (∞) and the vector of sums of probabilities of asymptomatic states

Such transcendental system should be solved numerically, and for large networks might be impossible task. However, one can at least obtain how the solution will look like, when the epidemic is weak in a sense that only small fraction of the population is infected during its course. In such case the probability of susceptibility will not change significantly p S (0) ≈ p S (∞). This situation might be present for example when the contagiousness of the pathogens is slightly over the threshold. Then one can keep only the leading terms in the expansion of the logarithm and obtain

The last approximation means that effectively the left hand sides of relationships (34) are equal. Then, after some algebra, by using (35), from those relationships one can obtain

The last relationship is eigenvalue equation of a matrix which is the adjacency matrix multiplied by the scalar

, which corresponds to eigenvalue equal to one. Thus, the vector of sums of probabilities of asymptomatic state (35) represents eigenvector of the adjacency matrix that corresponds to the eigenvalue Λ such that

To determine which is the eigenvalue Λ, observe that we can apply similar inequality as (9), which means that for each node i one has

This implies that one has the following vector inequality

which is obtained from (34) with simple algebra. When the epidemic parameters are such that

where Λ max is the largest eigenvalue of A, then (38) can not be satisfied for no one eigenvalue. That is the condition when the endemic equilibrium does not exist. To determine when it will emerge, one should increase the value of the fraction in (41), by modifying the epidemic parameters. Then, the first eigenvalue that can satisfy inequality as (38) will be exactly the largest eigenvalue Λ max . Thus the condition for existence of endemic equilibrium is

The last result is generalization of the case of random regular graph for which Λ max = k.

Because the leading eigenvalue of the adjacency matrix A determines the equation for p A , that vector is determined from the respective eigenvector, or the principal eigenvector of A. However, the eigenvalue equation (37) just determines the relative magnitudes of the components of p A . If it is used in the system (34) it will be obtained that for each node the change in the probability of being susceptible is the same which counters the fact that it corresponds to the respective component of the principal eigenvector. We further note that from (34) the relative magnitudes of the changes of the probabilities of susceptible state p S,i (0) − p S,i (∞), and the resulting probabilities of recovered state p R,i (∞), as collinear to p A , are also proportional to the principal eigenvector of the adjacency matrix. This is in accordance to the reasoning that the individuals with highest risk of infection are those with many contacts, and particularly those which have many high-degree neighbors. Thus the eigenvector centrality of the node [30] , which is the respective component in the principal eigenvector, determines the risk of infection of that node. We note that, by applying this procedure one can also show that the same conclusions about the role of the leading eigenvalue and principal eigenvector in epidemic spreading on complex networks hold for the simpler SEIR and SIR models.

We now proceed with the study of the stability of the disease-free equilibrium and determine the epidemic threshold. The associated Jacobian matrix is obtained by taking the respective derivatives in the equations (20) to (24) . Also, we remind that after making differentiation, at the epidemic inception in the Jacobian it should be taken p E,i = 0, p A,i = 0, p I,i = 0 and p S,i = 1. It can be verified that the Jacobian will have the following matrix form

where I is identity matrix of the same size N as the adjacency matrix of the network A -the number of nodes in the network. One can note the similarity in the structure between the last matrix and that in (29) . The eigenvalues of the last Jacobian are obtained from the characteristic equation R(λ) = det(J − λI 5N ) = 0, where we emphasize that the involved identity matrix has size 5N × 5N . One could use the approach given in [27] to determine the dependence of the epidemic threshold on the largest eigenvalue of the adjacency matrix. We have chosen alternative approach here, based on Schur's determinant identity

which is more general. Clearly, when a matrix has many zero submatrices, its application provides simpler results. By repetitive use of it, which is elaborated in the Appendix D, it can be shown that the nontrivial eigenvalues can be obtained from the polynomial corresponding to the following determinant

We note that the same determinant can be obtained by the procedure given in the Appendix F which even delivers the eigenvectors of the Jacobian. Currently, we cannot state whether the approach in Appendix D is just an alternative, or it might have potential to provide results when the latter is not useful. To continue with the analysis, one can substitute the multiplier of the identity matrix in the last equation as

and will obtain the characteristic function for determination of the eigenvalues of the adjacency matrix det(ΛI − A). From the relationship (46) it can be seen that for each eigenvalue of the adjacency matrix Λ correspond three eigenvalues of the Jacobian λ, which are obtained from the polynomial

The last relationship is identical to the respective one for the random regular graph (30) if one substitutes k with Λ. To determine which Λ corresponds to the largest λ max , observe first that the roots of the polynomial (47) are obtained in fact from cubic equation of the form

where in the last equation all other parameters except Λ are absorbed in the constants a, b and c. To find how the largest λ can be determined, consider two different eigenvalues Λ 1 < Λ 2 , to which the largest solutions of (48), λ 1 and λ 2 respectively, are both positive. To compare λ 1 and λ 2 , observe first that the cubic polynomial in (48) has shape that looks like the inverse of the letter N, which means that it is negative and decreasing from the largest root up to infinity. Now, by taking λ = λ 1 in the cubic polynomial for Λ 2 it will be obtained

The last inequality means that for Λ = Λ 2 , the value of the polynomial at λ = λ 1 has positive value, which means that it has zero λ 2 > λ 1 . Thus, the largest positive λ max corresponds to the largest eigenvalue of the adjacency matrix Λ max . When the epidemic can spread, the largest by modulus eigenvalue of the Jacobian is also positive for the following reason. The perturbations of all susceptible states at the disease-free state are either negative or zero and if λ max < −1, then after one iteration the perturbations will be negative and would make the probabilities of the susceptible states greater than one, which is impossible. Thus, the largest by modulus eigenvalue of the Jacobian (43), λ max , which determines the stability of the disease-free state is related to the largest eigenvalue of the adjacency matrix Λ max through the equation (46). For certain adjacency matrix one can analyze the stability of the disease-free state from the relationship (46) by using the leading eigenvalue Λ max for Λ. One can apply the same procedure as for the random regular graph and conclude that the disease-free equilibrium is stable once the following holds

The last relationship can be seen as generalization of the respective one for the random regular graphs (32) . Let us proceed with determination of the linear stability of the endemic state. To determine the respective Jacobian matrix, first observe the following derivatives

where A i,j is the i, j-th element of the adjacency matrix. The remaining partial derivatives in the Jacobian matrix are the same as for the disease-free state and are conveniently captured in the respective submatrices in (43). The form of the partial derivatives (51) is such that for each p S,i they have identical form and likewise for the p E,i . To obtain a more compact form for expressing such relationship, one can introduce a diagonal matrix Σ which contains the endemic equilibrium probabilities p S,i (∞) along the diagonal Σ i,i = p S,i (∞). Then, one can obtain that the partial derivatives between the susceptibility and exposed vectors with respect to the asymptomatic and infectious vectors can be compactly written as

Now, the Jacobian of the endemic equilibrium differs from that for the disease-free one, only in that it contains the matrix product ΣA instead of A. Respectively, the stability of the endemic equilibrium will depend on the leading eigenvalue L max of the matrix product ΣA. Thus, the stability condition is similar to that for the disease-free state (50)

In the Appendix E it is shown that in case of small contagiousness α ≪ 1 and β ≪ 1 and when epidemic affects small population during its course, the endemic equilibrium is linearly stable. Let us finally examine the behavior of the disease spreading in the early phase of epidemic. The one-step evolution of the probabilities of different states in disease spreading on complex networks is given by equations (20) to (24) . One can combine all probabilities in single column vector as p = [p T S , p T E , p T A , p T I , p T R ] T and the right hand sides of the probability evolution equations in a vector F . Then, one-step evolution of the probabilities can be compactly written in vector notation as

Consider early stages of the epidemics, when the states p(n) are sufficiently close to the disease-free equilibrium p DF = F [p DF ]. Denote with δp(n) = p(n) − p DF the deviation from the disease-free state. Then from (54) one has

The linear approximation of the nonlinear function F in vicinity of the disease-free state is

It means that consecutive perturbations satisfy simple relationship δp(n + 1) ≈ J DF δp(n).

Thus, at the eqarly phase of an epidemic, the perturbation at given moment n is approximately given as

Denote with z i the eigenvector of the Jacobian that corresponds to the eigenvalue λ i . Also, let the vectors z i are properly chosen to make orthonormal basis. Then the perturbation δp(0) can be expressed in terms of the Jacobian basis vectors as

Then, after n time steps the perturbation will approximately evolve to

It is clear that as the number of steps n increases, the projection along the principal eigenvector will dominate the others. It means that one can use the approximation

where p max is the projection of the initial perturbation along the principal eigenvector z max . In the Appendix F it is explained that when the epidemic starts the leading eigenvalue of the Jacobian λ max depends on that of the adjacency matrix Λ max and that the principal eigenvector of the Jacobian in that case is determined with the principal eigenvector of the adjacency matrix. Thus, the latter determines the evolution of the epidemics at the early stages. We emphasize that this is an approximation since as n grows, the Jacobian which is used, less accurately represents the nonlinear evolution, because the state of the system goes away from the disease-free one.

Although an approximation, the last result provides an estimate of the risk of being infected of the nodes in a network, by the respective eigenvector centrality. This observation needs more theoretical work for establishing the conditions when the principal eigenvector is really useful in estimation of the risk of infection.

The focus of numerical experiments in this work is put on validation of the theoretical results for the discretetime SEAIR model. The theoretical analysis of the compartmental model was classical and did not bring any significant novelty, which is not known for the other compartmental models. Thus, the potential of the compartmental model should be tested on real data, which is left for future study.

We have made simulations of disease spreading on random regular graphs by numerical solution of the evolution equations (26) . The aim of these numerical experiments was to check the validity of the epidemic threshold relationships, as well as the equation for the fraction of the susceptible individuals at the end of the epidemic (27) . For computational reasons, for this and the other numerical experiments in this work, we have considered networks with 1000 nodes. All versions of the SEAIR epidemic spreading model considered here, have five parameters α, β, σ, µ and γ. However, theoretical analysis in previous sections has shown that only the first four of them are relevant for determination of the epidemic threshold and the susceptible fraction at the end of the epidemic. The contagiousness parameters' values, α = 0.0025 and β = 0.002 < α were chosen arbitrarily, by caring to be small to ensure that the approximations made in the theoretical analysis are justified and using the observation that for COVID-19 contagiousness is bigger before the onset of symptoms. We have taken γ = 0.5 which should correspond to two days mean period of incubation, while the value of σ = 0.2 was chosen arbitrarily. The critical parameter µ 0 was calculated from the following quadratic equation

which is obtained from the condition for emergence of endemic state for the random regular graphs (28) . The value of the parameter µ was varied in vicinity µ 0 . All simulations were repeated for ten different networks and for each network ten different initializations were made by putting random node in exposed state (patient zero), while leaving the remaining ones as susceptible. The pathogen was considered as extinct at the moment when the total fraction of exposed, asymptomatic and infectious individuals is smaller than 10 −8 of the population. In the figure 1 are shown the average number of the susceptible individuals at the end of the epidemic. The number of susceptible individuals for each particular simulation is simply sum of the probabilities of the susceptible state over all nodes. Averaging is performed for all networks from the same type and for all initial conditions. In the blue diamonds are given the results for the random regular graph with node degree k = 50, while with red circles is the curve for the random graph with constant degree distribution in the interval [30, 70] . We emphasize that in this figure the threshold value µ 0 is obtained for the random regular graph and the same value is used for the other one. Both considered graphs have same average degree, and thus show similar results, particularly when one is far enough from the threshold µ = µ 0 . In vicinity of µ 0 , as was theoretically shown for general complex networks, the epidemic threshold depends on 27) for infinite-size random regular graph with node degree 50; blue diamonds -random regular graph with the same degree and 1000 nodes; red circles -random graph with uniform degree distribution in [30, 70] and 1000 nodes.

the leading eigenvalue (42) which for the graph with distributed node degree is greater than the average degree Λ max > k as the Perron-Frobenius theorem claims [31] . Thus, for same µ one expects more infected individuals for the graphs with distributed degree. The results from the simulations are further compared with theoretical values obtained from (27) for random regular graph with k = 50 nodes, which holds for network with infinite number of nodes. It can be seen a noticeable difference between the theoretical curve and the simulations. One reason for such discrepancy could be attributed in the fact that the theoretical results hold for infinitely large networks. The other factor is the initialization of the epidemics, which even in case of stable disease-free state, µ > µ 0 , produces at least a small fraction of infected individuals -at least the neighbors of the patient zero. Next, we have considered disease spreading on Erdős-Rényi (ER) [32] and Barabási-Albert (BA) [12] models of complex networks. Within the ER model, we have considered probability of existence of link between each pair of nodes p ER ∈ {0.01, 0.03, 0.05} and generated ten different networks for each case. For the BA complex networks we have taken four different values of the seed m ∈ {5, 10, 15, 20}. As for the random regular graphs, for each ER and BA network ten different initial conditions were considered. The parameter values for α, β, γ and σ were taken identical as for the random regular graphs. In the figure 2 are shown the average number of susceptible individuals and the correlation coefficient between the recovered probability vector and the principal eigenvector of the adjacency matrix at the end of epidemic. Here, the threshold value µ 0 was calculated for each network separately from the equation

where Λ max is the leading eigenvalue of the network which had to be calculated previously.

In the figure 2 can be seen that as the parameter µ increases towards the critical value µ 0 , the number of susceptible at the endemic equilibrium approaches the total number of individuals, as it is expected. We remark that it is not equal to the total population even when the conditions of epidemic are not satisfied, because there is certain probability that the patient zero will infect some neighboring nodes. However, this is finite size effect, and in infinitely large network the fraction of infected individuals will be infinitesimally small. We also remark that the results about the ER network look that it is more prone to epidemic. This deception appears because the horizontal axis is in the units of the threshold value µ 0 and not the absolute terms. As it is known, BA networks are more prone to epidemics and for infinite size networks the epidemic threshold is vanishing [6] .

The rather high value of the correlation coefficient ρ, when disease is spreading suggests that indeed the principal eigenvector predicts the pattern of infection. When the epidemic is not possible, µ > µ 0 , the correlation does not drop sharply, due to the finite size effects. Near the epidemic threshold there is nonzero probability of infecting the neighborhood by the initially exposed node, and particularly those with higher eigenvector centrality.

We have finally studied the behaviour of the epidemics at the onset in order to verify which nodes bear the highest risk of contracting the disease. As common wisdom suggests, highly connected nodes, and particularly those with well connected neighbors are most risky ones -just as eigenvector centrality ranks the nodes. For that reason we have calculated the evolution of correlation coefficient between the principal eigenvector of the adjacency matrix and the probability vector of recovered state as the epidemic unfolds. In the figure 3 is shown the correlation coefficient as function of time. For the BA network shown at right the parameters have the same values, while for the ER network (at left panel) α = 0.05 and β = 0.04, while γ and σ are the same as in the other simulations. One can note that generally in the early stages of the disease outbreak very high correlation is achieved, which confirms that the eigenvector centrality predicts rather well the riskiness of contraction of the pathogen. As the epidemics fades out the correlation might drop because in certain parameter combinations majority of population has high chance of becoming infected and this infection pattern can differ significantly from the predictions by the principal eigenvector of the adjacency matrix. However, the lowest curve for the ER network model shows that this is not always happening. In such situation, when epidemic is barely possible, only small fraction of population can be affected, particularly those which are close to the patient zero. This observation needs further investigation of the pattern of risk in case of such small outbreaks.

We have studied SEAIR epidemic spreading model inspired by the contagiousness features of COVID-19. Theoretical analysis were made for the compartmental version as well as for discrete-time epidemic spreading on random regular and complex networks. For the compartmental model the epidemic threshold was found and it was shown that it also determines the emergence of endemic equilibrium as well. When the contagiousness is weak, we have shown that for random regular graph and complex network the epidemic threshold obtained from stability analysis of the disease-free state depends in a similar way on the model parameters. As is known for many other disease spreading models, the epidemic threshold was obtained to depend on the leading eigenvalue of the adjacency matrix. We have also demon-strated that the endemic equilibrium when exists is linearly stable in the three considered models. Theoretical analysis in this work has shown that the risk to be infected certain node is dependent on its eigenvector centrality. In early stages of epidemics, the eigenvector centrality points which nodes are most likely to be first to contract the disease, while in case of mild epidemic on complex network, it shows which nodes have more risk to contract the disease during the whole course of the epidemic.

The analysis of the linear stability was based on two approaches. By appropriately organizing the probabilities of various states as dynamical variables it was obtained Jacobian matrix of the equilibria which has form that allows analytical treatment. The first approach was based on applying Schur's determinant identity which lead to result that the nontrivial eigenvalues of the Jacobian are related to those of the adjacency matrix. In the second approach we have furthermore shown that eigenvectors corresponding to the nontrivial eigenvalues of the Jacobian are combinations of scaled eigenvectors of the adjacency matrix. We believe that these two techniques can be applied in a range of studies where multidimensional dynamical systems interact through complex topology of contacts.

Although the motivation of studying the SEAIR epidemic spreading model was the COVID-19 pandemic, we did not make any testing about its relevance on real data. Naturally, it is the first study which should follow this one. One of the key issues would be inference of the fraction of the population which has contracted the disease, but has not shown symptoms at all. This could help in estimating the likelihood of reappearance of the epidemic and its possible size, once it weakens. Since, in general, contagiousness parameters change during epidemic, testing the validity of the relationships for the fraction of susceptible individuals at the end of epidemic might not be easy task. However, at the early phase of an epidemic these parameters could be considered as constant. Then, by using real data, it could be verified how well the eigenvector centrality anticipates which individuals bear highest risk of infection. If it proves to be useful predictor, then a follow up is investigation of its relevance to planning of vaccination and other protective measures.

As is given in the main text, the eigenvalues of the Jacobian at the disease-free state of the compartmental model are obtained from the determinant det(J DF − λI), that is

Besides the two trivial eigenvalues λ = 0, the remaining three are the roots of the polynomial which is obtained by expanding the last determinant

The cubic polynomial in λ in the last equation can be written in the form

where the coefficients are

By the Routh-Hurwitz criterion [28] , the roots of the third order polynomial of the form (A3) will have negative real parts if and only if a 2 > 0 and a 2 a 1 > a 0 > 0.

Since a 2 > 0, the condition a 0 > 0 is equivalent to

We should also verify that also a 2 a 1 > a 0 is satisfied which after multiplication of the respective values in (A4) will result in γµσ + γ 2 σ + γµ 2 + 2γ 2 µ − αγ 2 +µσ 2 + γσ 2 + µ 2 σ + 2γµσ − αγσ

By algebraic manipulation and rearrangement of the last inequality becomes

Now, from the condition (A5) one has the inequality

which implies the following relationships for the negative sign terms in (A7)

By using the last three inequalities in (A7), one will obtain that only positive terms will remain at the left hand side, which means that it is satisfied. Thus by the Routh-Hurwitz criterion the nontrivial eigenvalues of the Jacobian have negative real parts if and only if (A5) holds.

Appendix B: Endemic equilibrium for random regular graph

For determination of the population remaining unaffected by the epidemic in disease spreading on random regular graph, we follow the same approach as in the compartmental model. To proceed in that spirit, first sum up the first four equations in the system (26) and obtain p S (n + 1) + p E (n + 1) + p A (n + 1) + p I (n + 1)

We can sum the last relationship over all moments n, from the onset to finish of the epidemic, and assume negligibly small initial probabilities of infected individuals. Then, due to cancellation of the respective terms it will be obtained

Appendix D: Characteristic polynomial for the eigenvalues of the Jacobian of the discrete-time model

To obtain more compact expression for determination of the eigenvalues of the Jacobian (43), we will extensively use the Schur's determinant identity

One should note that the identity does not need the matrices to be square and if at least one of the matrices R or S is zero, then one has simpler relationship det Q R S T = det(Q) · det(T).

(D2)

First we can assign the role of the bottom-right submatrix T in (D1) to the bottom-right identity matrix in (43). Then one can note that to the respective submatrix R corresponds zero matrix and use (D2) instead to obtain

By repeating the same procedure one more time with taking top-right submatrix (1 − λ)I as the submatrix Q in the Schur's determinant identity, and observing that now the submatrix S is zero, one can obtain that

To simplify notation one could first stop repetitive writing of the part which contains the trivial eigenvalue λ = 1 which has multiplicity 2N , and focus on the remaining. Take the submatrix T = (1 − µ − λ)I which determinant contains trivial eigenvalues λ = 1−µ and respectively the remaining submatrices Q, R and S. We note that 1 − µ are not eigenvalues of the Jacobian, since in expanding the determinants as polynomial, the terms corresponding to 1 − µ − λ that appear in T will cancel with the same terms which will appear in the denominator in the remaining determinant as will be seen below. From the last determinant let us first consider the submatrix that corresponds to the product RT −1 S in the Schur's identity (D1). By using the properties of the inverse matrix one can obtain first Now, the part of the characteristic polynomial which contains the nontrivial eigenvalues is

We can apply the Schur's identity again. First observe the matrix product that corresponds to the RT −1 S term

After simplification of the scalar at the right-hand side of the last relationship and subtract the respective matrices in the form Q − RT −1 S from (D1) one will obtain the following characteristic polynomial of the eigenvalues

Again, the first determinant has trivial eigenvalues λ = 1 − σ − µ with multiplicity N as well and the nontrivial ones are contained in the second determinant. By observing the second determinant in (D9) one can note that in the denominator multiplying the adjacency matrix appear terms 1 − µ − λ and 1 − σ − µ − λ. Expansion of the determinants as polynomials will result in cancellation of those terms in the denominators with the respective ones in the determinants det [(1 − µ − λ)I] and det [(1 − σ − µ − λ)I]. Finally, the characteristic polynomial resulting from the last nontrivial determinant will not change if one multiplies it with a constant. So, a more convenient form of the last determinant, and the respective characteristic polynomial is

(D10)

Appendix E: Stability of the endemic equilibrium in disease spreading on complex networks

Since the Jacobian of the endemic and disease-free equilibrium differ only in the presence of the matrix Σ, the characteristic equation will have the same form for both cases. However, it was previously obtained that the leading eigenvalue of the Jacobian of the diseasefree equilibrium λ max depends on the leading one of the adjacency matrix Λ max . Accordingly, for the endemic state the dependence will be on the leading eigenvalue L max of the matrix product ΣA. We will verify that this eigenvalue is related with that of the adjacency matrix as L max < p S,max Λ max , where p S,max = max p S,i (∞), is the maximum of the probabilities of susceptible states at the end of the epidemic. To prove that, denote with x the unit eigenvector of ΣA, corresponding to L max , or ΣAx = L max x. Let Λ i and u i are the eigenvalues and the respective orthogonal basis vectors corresponding to the adjacency matrix. The vector x in the basis u i is given as

where a 2 i = 1 because x is unit vector. Then, multiplying the matrix A with x will result in some vector

Due to the orthonormality of the basis u T i u j = δ i,j , the squared magnitude of y reads

which can be bounded as

This means that the vector y has length not bigger than Λ max . In connected network each node will be infected with nonzero probability p S,i (∞) < 1. Thus the matrix Σ is symmetric positive semi-definite, and all its eigenvalues are strictly less than one. Let us now express the vector y in the orthonormal basis v i of the matrix Σ

Then the vector Σy can be expressed as

since Σ is diagonal matrix with eigenvalues p S,i (∞). The squared magnitude of Σy is bounded as

Now, combining (E2), (E4) and (E7) will result in

Thus, we have just bounded the leading eigenvalue of the matrix ΣA as

Recall that in the stability analysis of the endemic equilibrium one has the matrix product ΣA instead of A which is used in the disease-free state. So, the stability of the endemic equilibrium depends on L max as the other case depends on Λ max . Correspondingly, the endemic equilibrium will be linearly stable, once the following inequality holds [refer to the respective condition (50)] µ(µ + σ) > L max (αµ + βσ).

(E10)

Now, consider the system of transcendental equations (34) and use the fact that p A is principal eigenvector of the adjacency matrix A, or Ap A = Λ max p A . By al-

Histoire de l'Acad

The prevention of malaria

22nd International Symposium on Reliable Distributed Systems

This research was partially supported by the Faculty of Computer Science and Engineering, at the SS Cyril and Methodius University in Skopje, Macedonia.

Summation of the infinite ladder of equations (B3) and using p I (0) ≈ 0 = p I (∞) leads to the same relationship between the probabilities of asymptomatic and infectious states as in the case of the compartmental model (7),The first equation in (26) can be rearranged as p S (n + 1) p S (n)If we take logarithm of the last equation, and use approximation αp A (n) ≪ 1 and βp I (n) ≪ 1, that holds for weak spreading α ≪ 1 and β ≪ 1, it will be obtainedSumming the last relationship for all moments, after cancellations results inUsing the relationships (B2), (B4) and (B7) one can obtain an estimate of the number of unaffected individuals in epidemic spreading on random regular graphs fromAppendix C: Endemic equilibrium for complex network For small contagiousness parameters α ≪ 1 and β ≪ 1, one can approximate the probability that susceptible individual will not receive the virus asThen the evolution of all probabilities can be compactly written asWe will follow the same technique as for the previous two scenarios. Summing up the first four equations in the last system will result in p S (n + 1) + p E (n + 1) + p A (n + 1) + p I (n + 1)(C3) Now, lets sum over all moments and use the fact that the probabilities of infected states at the beginning and ending of epidemic are vanishing. Then from the last relationship will be obtainedRearrangement of the fourth equation in (C2) and summing over all moments will lead to result that generalizes (B4)One can write the evolution equation of probability of susceptible state for each node i asFurther, take logarithm on both sides of the last equation and keep only leading terms in α and β in the expansion of the logarithm of the multipliers to obtainSumming (C7) over all moments will result in(C8) Denote with ln p S (n) the vector which components are the logarithms of probabilities of susceptible states ln p S,i (n). Then, the relationship (C8) for all nodes can be compactly written asFrom another side, applying (C5) in (C9) will lead toThe last two relationships are system of equations for determination of vectors of the probabilities of the susceptible state and the infinite sum of the vectors of asymptomatic states. gebraic manipulations, from the system (34) it can be shown that for each component of the susceptible probability vector holds relationship similar to (8) from which one hasCombining the endemic equilibrium stability condition (E10) with the last relationship (E12) will result inRearranging the terms in the last inequality will result in more convenient formThe last inequality is satisfied since one can use (E9) and the fraction at the left-hand side is always greater than one. Thus, when an epidemic occurs such that small fraction of the population is affected, the respective endemic equilibrium is linearly stable.Appendix F: Eigenvalues and eigenvectors of the Jacobian at the disease-free state for epidemic spreading on complex networksDenote the eigenvectors of the Jacobian matrix in the disease-free equilibrium with w = [w S , w E , w A , w I , w R ] T where w S , w E , w A , w I and w R are the parts which correspond to the states S, E, A, I and R respectively. Then, the eigenvalue equation for the Jacobian Jw = λw in more detailed form isFrom the fourth row in (F1), which corresponds to the infectious state, one can obtain thatwhich can be rearranged intoThe last equation relates the magnitudes of the vectors w A and w I and shows that they are collinear. In similar manner, from the last row in (F1), one can show that the vector w R is collinear with the previous ones and moreoverLikewise, from the third row in (F1) it follows that the exposed probability vector w E is also collinear to the previous ones, or more preciselyNow consider the second row in (F1), from which one haswhich after using (F3) and (F4) will result in(F7) The last relationship could be rearranged further as(F8) We have obtained eigenvalue equation for the adjacency matrix. Thus, every vector w A must be eigenvector of the adjacency matrix A. Since the eigenvalues of the adjacency matrix Λ are independent on any dynamical process evolving on the network, it means that the eigenvalues of the Jacobian λ must satisfy the relationshipThe last result relates the eigenvalues of the Jacobian with those of the adjacency matrix. By expanding the terms, one can see that it is cubic polynomial in λ, and thus for each eigenvalue Λ one has three possibly different eigenvalues λ. Thus, N eigenvalues of the adjacency matrix would generate 3N eigenvalues of the Jacobian. We remind that as is given in the Appendix D, there is one trivial eigenvalue λ = 1 with algebraic multiplicity 2N . The eigenvectors corresponding to this eigenvalue are those that span the subspace consisting of susceptible or recovered states only and zeros at the remaining states. It can be easily verified from (F1), that each vector of the formis eigenvector of the Jacobian.Finally, from the first row in (F1) it follows thatin which one can use (F8) to obtainThus, w S vector is collinear with the rest as well. This implies that besides the vectors (F10), the remaining eigenvectors of the Jacobian w consist of scaled copies of the eigenvectors of the adjacency matrix. More precisely, by using the relationships (F3), (F4), (F5) and (F12), the eigenvector w isSince w A is eigenvector of the adjacency matrix, there are 3N linearly independent eigenvectors of the form given in (F13). Together with 2N linearly independent vectors of the form (F10), they constitute a set of 5N linearly independent vectors. When properly normalized they can be made to form an orthonormal basis.