SEAIR Epidemic spreading model of COVID-19.

We study Susceptible-Exposed-Asymptomatic-Infectious-Recovered (SEAIR) epidemic spreading model of COVID-19. It captures two important 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.

Introduction

Understanding epidemic spreading of contagious diseases and effectiveness of various countermeasures is of high interest for the public health and the society in general, with important contributions provided by epidemiologists, mathematicians and physicists as well. Although earliest theoretical work in the 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, some associations were found between the epidemic threshold and the network properties like the degree distribution or the leading 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 diseases 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 in this setting are systems of differential equations for the evolution of the size of those compartments. In such models, there is an 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 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 one should mention the studies of disease localization [16], [17], spreading on multilayer networks [17], [18], the assumption of non-exponential distribution of periods between consecutive events [19], [20] and the effects of the delay in the recovery [21] and in the infection [22].

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], [23]. 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 [23].

In this theoretical work we consider 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, 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 one based on statistical physics framework, in which is considered discrete-time epidemic model on complex networks. We study the latter model on random regular graphs and on complex networks separately. For the three versions we obtain the epidemic threshold, equation for the fraction of susceptible individuals at the end of epidemic and we study the linear stability of the disease-free and endemic equilibria. The results for the random regular graphs hold when the contagiousness is weak, while for the complex networks it is also needed that epidemic is small-scale one in a sense that it can affect a tiny fraction of the population. We furthermore study the roles of the leading eigenvalue and the principal eigenvector of the adjacency matrix in the spreading. It is already known that the leading eigenvalue determines the epidemic threshold [14], [24]. We show that the principal eigenvector and its associated eigenvector centrality have important role in estimation of the risk of infection. We finally note that, the techniques which are applied for analysis of the Jacobian matrix at the equilibria might be relevant in studies of linear stability of coupled multidimensional dynamical systems.

The paper is organized as follows. First, in Section 2 we present some of the relevant literature about epidemic spreading models of COVID-19 and their relation to the proposed model. Then, in Section 3 we introduce and analyze the SEAIR compartmental model. In the following Section 4 is studied the discrete-time version on random regular graphs and complex networks, while in Section 5 we present some results of the numerical experiments. The paper ends with the Conclusion.

Epidemic spreading models of COVID-19 and their relation to the SEAIR model

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 established that a key feature is the possibility that an infected person can spread the disease before onset of symptoms or without having them at all [25], [26]. Another important observation is that the median incubation period of the disease is approximately five days [27]. Thus, the asymptomatic spreaders and individuals in the incubation period are hidden disease carriers. Their presence poses a challenge in the control of the epidemic, in the planning of healthcare capacity buildup or the relaxation of lockdown measures, and even for estimation of the population affected by the pathogen. The mathematical models of the epidemic spreading are adapted accordingly to capture the key features of the spreading of the COVID-19, for which the most basic and popular ones, SIS and SIR are not satisfactory. The existence of incubation period is addressed by including compartment of exposed individuals (E), which in certain studies in the literature is denoted as latent. The absence or presence of symptoms is addressed in different ways in the literature. In one approach [28], it is considered presence of presymptomatic, or prodromic phase, where the individual is already infectious, which is followed by either asymptomatic or symptomatic phase which can have three different severity levels. In another study [29], there are separate latent compartments for the asymptomatic and symptomatic individuals, and it is considered different contagiousness by the infected individuals depending on the presence or absence of symptoms. The observation of existence of super spreaders of the disease was also included in a mathematical model [30]. In the analysis performed in [31], the infected persons which are detected, either asymptomatic or symptomatic, are accounted in separate class from those that are not. Since the quarantine is one of the key defence measures against the spreading, the quarantined compartment was considered as well in many works [32], [33], [34], [35]. Also, to account for the disease outcomes and the burden on the healthcare capacities were introduced compartments for the hospitalized, for those under intensive care, and for those which did not survive the disease [28], [29], [34], [35].

From one side, the simpler models are easier to study and in general provide better estimation of the model parameters which are needed for making predictions for future development of an epidemic. From another side, the more complex models obviously allow for better description of certain features of particular epidemic like the COVID-19. Thus, one has to make a choice about what kind of model to apply in order to obtain particular insight. In this work, the aim is to work with a model that will address the key features of COVID-19 epidemic, the latency and the existence of asymptomatic transmission, and study it in frames of complex networks theory which accounts for nontrivial pattern of contacts between the individuals. Accordingly, we separate the disease transmitters in two compartments. In the asymptomatic one is included an individual which do not have symptoms at the moment of transmission of the virus, regardless whether he or she will obtain them later or not. In the infectious compartment are accounted the persons which can infect the others and have symptoms already. Thus, in the asymptomatic and infectious compartments are included the three kinds of disease spreaders: always asymptomatic, currently asymptomatic which will develop symptoms later and symptomatic. More detailed description of the SEAIR model are given in the following two sections. We finally note that the proposed model can be considered as similar to that in [36], where the individuals which are not detected as carriers of the virus are considered as unreported to the authorities and can transmit the disease, while the reported, or detected ones, do not contribute to the spreading. Another model closely related to SEAIR is that in [37], because it also has incorporated pre-symptomatic transmission and two different compartments for the asymptomatic and symptomatic individuals.

SEAIR compartmental model

Let the variables and 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 a 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 be while of infectious ones be . 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 the 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 modelWe 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, .

One trivial solution of the system (1) is the disease-free state 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 which will be applied here [38]. For the SEAIR model, if we sum the top four equations in (1), the following relationship will holdIn 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 aswhich by integration will result in another relationship between the initial and the final fractions of susceptibles

One can also integrate the fourth equation in (1) on both sides to obtainThe last result provides a relationship between asymptomatic and infectious fractions in the course of the whole epidemicBy combining the relationships (3), (5) and (7), the following equation for the fraction of unaffected individuals is obtainedUsing the fact that is steeper than for one can verify thatThis implies that the transcendental Eq. (8) has a solution only ifThe last inequality is the condition of existence of endemic equilibrium of the system (1).

To study the linear stability of the equilibrium states, one should linearize the system (1). The respective Jacobian matrix readsAt the disease-free state for which and the Jacobian has rather simple formBecause the first and last columns are zero, this Jacobian has two trivial zero eigenvalues, while the other three are the roots of the polynomialWe show in Appendix A that the three nontrivial eigenvalues of the Jacobian have negative real part, which implies linear stability of the disease-free state, if the following relationship holdsThe obtained inequality is opposite of the condition for the existence of endemic equilibrium, as one can expect. If (14) holds, (10) does not, the disease-free state is stable and epidemic will not occur. In the opposite, if (14) is not satisfied, the equilibrium 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 disease-free one. In this case, in the Jacobian (11) one should take and . 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 and respectively. This would simply change the condition for stability of the endemic equilibrium toFrom (8) one has the following relationship for the fraction containing the parametersPlugging the last relationship in (16), it will be obtained that the endemic equilibrium is stable once the following holdsThe last inequality can be rearranged aswhich holds always since .

Discrete-time SEAIR model

Subsequent investigation of a disease spreading model when one needs to account for the contacts between the individuals, is to study epidemic spreading on complex networks framework. Let us consider discrete-time evolution version of the proposed SEAIR model with finite population of 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. An exact approach for analysis of contact-based spreading on complex networks is based on using indicator random variable for each state of every node in the network and working with a system of size . However, as it is elaborated in details in [39], under the assumption that a node can be in certain state independently of the states of the other nodes, one can instead use the probability of being at that state as a more convenient variable of interest. We proceed in that spirit and build our model on the set of probabilities for each node being in certain state, S, E, A, I, or R, at given moment . We will also denote the parameters with the same Greek letters as in the compartmental model. They correspond to the same transitions and have 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 . 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 the infectious state.

The dynamics of the discrete-time version of the SEAIR model is built similarly to the model considered in [14]. Denote the probabilities that the individual at the discrete moment is in respective state with and . Our model assumes reactive process of epidemic spreading, which means that in each time step every individual has contact with every neighbor [40]. Then, certain susceptible person at the moment will not receive the infection from any of its neighbors with probability [14], [40] where denotes the set of neighbors of . The individual will remain susceptible at the moment if he, or she, did not receive the contagion, which means thatOtherwise, 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 probabilityThe probability of being asymptomatic at the next moment iswhere the last term accounts for the situation that neither the symptoms will appear, nor healing will happen in one time step, which imposes a restriction . The node will be in state I at the moment with probabilitywhere 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 isThe set of Eqs. (21) to (25) determines a discrete-time dynamical system of equations for evolution of probabilities of the states for each node in the network. It can be solved numerically for arbitrary initial condition and one can thus observe 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.

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 . 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 (20)will be simplified toThen 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 and the other when the fraction of susceptible individuals is such that it prevents further spread of the disease . In the Appendix Bwe 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 the epidemic. It can be applied when the contagion probabilities and are very small. As it is shown in the Appendix B, the equation for determination of the probability of the susceptible state at the end of epidemic is very similar to the respective one for the compartmental caseThe last result extends the one for all-to-all coupling considered in compartmental models, where effectively each individual can get the disease from anyone 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 isOne can note that the last relationship is similar to the respective one for the compartmental model (10), and the only difference is the presence of the node degree in the discrete-time case.

We can further make linear stability analysis of equilibria by linearizing the evolution equations. The respective Jacobian matrix at the fixed points for which and isSolving the characteristic equation of the Jacobian (30), for will result in two trivial eigenvalues equal to one and other three. We note that in the discrete-time case the trivial eigenvalues have value equal to one that also corresponds to the marginal stability, which are zero in the continuous-time scenario. The nontrivial eigenvalues are the roots of the polynomialOne can notice the similarity to the respective characteristic polynomial for the compartmental case (13). The only difference is presence of the degree which in the latter expression multiplies and and one has in the discrete-time case instead of . In the Appendix Cit is shown that the disease-free state is stable once the following relationship is satisfiedwhich is similar to the one for the compartmental case. Again, the only difference is the presence of the node degree .

The linear stability of the endemic equilibrium is established from the leading eigenvalue of the same Jacobian matrix (30) as the disease-free one, but using obtained from (28). The procedure is nearly the same as for the disease-free equilibrium and the only difference is that instead of one should use in all analysis. Then the endemic equilibrium will be stable, if the condition similar to (32) is satisfiedWithout showing the details, we will just mention that once the endemic equilibrium in this discrete-time disease spreading model exists, it is linearly stable.

Epidemic spreading on complex networks

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 of the whole system by stacking the state vectors of each node one on top of another. In this case another ordering is more appropriate [24]. First, create vector of the probabilities of susceptible states of all nodes then those of the exposed states and likewise for the remaining three and . Also, denote with the adjacency matrix of the network of contacts between the individuals with elements only if nodes and are neighbors and 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 and 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 D, when and the susceptibility probability vector can be calculated from the following self-consistent systemThe solution of the system of transcendental Eq. (34) consists of the susceptibility vector at the endemic equilibrium and the vector of sums of probabilities of asymptomatic statesSuch 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 . 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 obtainThe 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 obtainThe 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 the probabilities of the asymptomatic state (35) represents eigenvector of the adjacency matrix that corresponds to the eigenvalue such thatTo determine which is the eigenvalue observe that we can apply similar inequality as (9), which means that for each node one hasThis implies that one has the following vector inequalitywhich is obtained from (34) with simple algebra. When the epidemic parameters are such thatwhere is the largest eigenvalue of then (38) can not be satisfied for no one eigenvalue. That is the condition when 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 the equation as (38) will be exactly the largest eigenvalue . Thus the condition for existence of endemic equilibrium isThe last result is generalization of the case of random regular graph for which .

Because the leading eigenvalue of the adjacency matrix determines the equation for that vector is determined from the respective eigenvector, or the principal eigenvector of . However, the eigenvalue Eq. (37)just determines the relative magnitudes of the components of . 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 and the resulting probabilities of recovered state as collinear to 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 [41], 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, SIR and SIS models.

It should be emphasized that, there are works in the literature which point that the principal eigenvector of the adjacency matrix can have important role in disease spreading. The probability of infectious state of a node in SIS spreading model on complex networks was found to be proportional to the eigenvector centrality, in vicinity of the epidemic threshold [16], similarly as the analysis above claims. This finding was further extended to multilayer networks [17]. It was reported in the same contributions, that if the structure of the network is such that the principal eigenvector is localized, the spreading will be limited to small number of nodes, even for large networks. Thus, the association between the eigenvector centrality and probability of becoming infected might not be observed in certain scenarios. Such example could be observed, for example, in structured network with communities, when the virus starts spreading in a node in one community, while principal eigenvector has significantly large components in the other communities. Thus, more research on this issue is needed for better understanding of the conditions when the principal eigenvector is really useful in estimation of the risk of infection. Finally, it is worth noting that it was found that the principal eigenvector of another matrix – the submatrix corresponding to the infectious states – also determines the disease spreading pathways. This observation has appeared in the studies of disease spreading between spatial regions in a waterborne disease [42] and COVID-19 [34].

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 Eqs. (21) to (25). Also, we remind that after making differentiation, at the epidemic inception, in the Jacobian it should be taken and . It can be verified that the Jacobian will have the following matrix formwhere is identity matrix of the same size as the adjacency matrix of the network – the number of nodes in the network. One can note the similarity in the structure between the last matrix and that in (30). The eigenvalues of the last Jacobian are obtained from the characteristic equation where we emphasize that the involved identity matrix has size . One could use the approach given in [24]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 identitywhich 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 E, it can be shown that the nontrivial eigenvalues can be obtained from the polynomial corresponding to the following determinantWe note that the same determinant can be obtained by the procedure given in the Appendix G which even delivers the eigenvectors of the Jacobian. Currently, we cannot state whether the approach in Appendix E 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 asand will obtain the characteristic function for determination of the eigenvalues of the adjacency matrix, . From the relationship (46)it can be seen that to each eigenvalue of the adjacency matrix correspond three eigenvalues of the Jacobian which are obtained from the polynomialThe last relationship will be identical to the respective one for the random regular graph (31), if one substitutes with . Because all eigenvalues of the adjacency matrix are real [43], the coefficients in the last polynomial in are also real for each . Instead of checking whether the roots of the last polynomial are within the unit circle, one can determine the eigenvalues of the related polynomial that corresponds to the Jacobian of the compartmental model (13), and then use the relationship between the roots of the characteristic polynomials of the compartmental and discrete-time models (C.5). The disease-free state of the model on complex network is unstable, if there is at least one eigenvalue of the adjacency matrix, for which there is a real positive root of the polynomial (47) in which one should use instead of . This situation happens when parameters are such that the following inequality holdsIn determination of the epidemic threshold in SIS and SIR models, usually the fraction between the contagiousness and recovery parameters is varied. Since the last inequality is a bit more complex than those in the SIS and SIR models, for the SEAIR model one can use the parameters and as bifurcation parameters. Then from the last inequality the smallest and for which the disease-free state is unstable are those obtained for the largest that is the leading eigenvalue of the adjacency matrix . Thus, the same condition (42)implies existence of endemic equilibrium and instability of the disease-free state. If the opposite is true,then the disease-free state will be stable. The last inequality 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 derivativeswhere is the 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 (50)is such that for each they have identical form and likewise for the . To write a more compact form for expressing such relationship, one can introduce a diagonal matrix which contains the endemic equilibrium probabilities along the diagonal . 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 asNow, the Jacobian of the endemic equilibrium differs from that for the disease-free one, only in that it contains the matrix product instead of . Respectively, the stability of the endemic equilibrium will depend on the leading eigenvalue of the matrix product . Thus, the stability condition is similar to that for the disease-free state (49) In the Appendix F it is shown that in case of small contagiousness and 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 Eqs. (21) to (25). One can combine all probabilities in single column vector as and the right hand sides of the probability evolution equations in a vector . Then, one-step evolution of the probabilities can be compactly written in vector notation asConsider early stages of the epidemics, when the states are sufficiently close to the disease-free equilibrium . Denote with the deviation from the disease-free state. Then from (53) one hasThe linear approximation of the nonlinear function in vicinity of the disease-free state iswhere is the Jacobian at the disease-free state. It means that consecutive perturbations satisfy simple relationshipThus, at the early phase of an epidemic, the perturbation at given moment is approximately given as

Denote with the eigenvector of the Jacobian that corresponds to the eigenvalue . Consider situation when the eigenvectors constitute an orthonormal basis, in which the perturbation can be expressed in terms of the Jacobian basis vectors asThen, after time steps the perturbation will approximately evolve toIt is clear that as the number of steps increases, the projection along the principal eigenvector will dominate the others. It means that one can use the approximationwhere is the projection of the initial perturbation along the principal eigenvector . In the Appendix Git is explained that when the epidemic starts the leading eigenvalue of the Jacobian depends on that of the adjacency matrix 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 grows, the Jacobian which is used, represents the nonlinear evolution less accurately, because the state of the system goes away from the disease-free one. Although being 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.

Numerical experiments and discussion

The focus of numerical experiments in this work is put on validation of the theoretical results for the discrete-time 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 Eq. (27). 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 (28). 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, and 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 which should correspond to two days mean period of incubation, while the value of was chosen arbitrarily. The critical parameter was calculated from the following quadratic equationwhich is obtained from the condition for emergence of endemic state for the random regular graphs (29). The value of the parameter was varied in the vicinity of . All simulations were repeated for ten different networks and for each network ten different initializations were made by putting a randomly chosen 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 of the population. In the Fig. 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. The averaging was 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 while with red circles are those for random graph with constant degree distribution in the interval . We emphasize that in this figure the threshold value is obtained for the random regular graph and the same value is used for the others. Both considered kinds of graphs have the same average degree, and thus show similar results, particularly when one is far enough from the threshold . In vicinity of as was theoretically shown for general complex networks, the epidemic threshold depends on the leading eigenvalue as is given in (42), which for the graph with distributed node degree is greater than the average degree as the Perron-Frobenius theorem claims [43]. Thus, for the 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 (28) for random regular graph with 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 way of initialization of the epidemics, which even in case of stable disease-free state, produces a small fraction of potentially infected individuals – at least the neighbors of the patient zero.

Fig. 1

Disease spreading on random regular graph and random graph with constant degree distribution.The curves represent the dependence of the number of susceptible individuals at the end of the epidemic on the parameter . The meaning of the symbols is the following: orange stars – theoretical values from eq. for (28) 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. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Next, we have considered disease spreading on Erdős-Rényi (ER) [44] 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 and generated ten different networks for each case. For the BA complex networks we have taken four different values of the seed . 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 Fig. 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 was calculated for each network separately from the equationwhere is the leading eigenvalue of the given network.

Fig. 2

Discrete-time epidemic model on complex networks at the end of epidemic. In all panels ten different complex networks with 1000 nodes are considered. The seed of generating the BA networks, and the link probability for ER networks is given in the inset. In the top panels are shown the average number of susceptible individuals while at bottom are average correlation coefficients between the number of recovered individuals and the principal eigenvector of the respective adjacency matrix. The horizontal axis is given in units of the critical value of the parameter at the epidemic threshold which is calculated for each network separately.

In the Fig. 2 can be seen that as the parameter increases towards the critical value 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 is expected to be infinitesimally small in general. We remind that, in the related case, when the epidemic threshold is barely passed, for networks in which the principal eigenvector of the adjacency matrix is localized, only infinitesimal fraction of the population will be affected [16], although for general networks it will be finite. One can also notice 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 and not the absolute terms. We emphasize that, as it is well known, for infinite size BA networks the respective leading eigenvalue is infinite, and thus the threshold value of the contagiousness parameter is vanishing [6]. In our analysis of such networks, where the parameter is chosen to be varied, its critical value diverges for infinite networks.

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, 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 the eigenvector centrality ranks the nodes. For that reason we have calculated the evolution of the correlation coefficient between the principal eigenvector of the adjacency matrix and the probability vector of the recovered state as the epidemic unfolds. In the Fig. 3 is shown the correlation coefficient as function of time. For the BA network shown at right the parameters have the same values as previously, while for the ER network (at left panel) and 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 for 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 suggests further investigation of the pattern of risk in case of such small outbreaks.

Fig. 3

Evolution of correlation coefficient between the principal eigenvector of the adjacency matrix and the vector of probability of recovered state in ER (left panel) and BA (right panel) complex networks. The considered networks have 1000 nodes. Each curve is obtained by averaging ten networks with the same parameters and ten randomly chosen initially infected nodes for each network.

Conclusions

We have studied SEAIR epidemic spreading model aimed to capture 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 graphs 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 graphs and complex networks 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 demonstrated that when endemic equilibrium exists, it is linearly stable in the three considered models. The theoretical analysis in this work has shown that the risk for infection of 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 the 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 for 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.

CRediT authorship contribution statement

Lasko Basnarkov: Conceptualization, Methodology, Software, Writing - original draft, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.