Novel Method to Analytically Obtain the Asymptotic Stable Equilibria States of Extended SIR-type Epidemiological Models

: We present a new analytical method to ﬁnd the asymptotic stable equilibria states based on the Markov chain technique. We reveal this method on the SIR-type epidemiological model that we developed for viral diseases with long-term immunity memory pandemic. This is a large-scale model containing 15 nonlinear ODE equations, and classical methods have failed to analytically obtain its equilibria. The proposed method is used to conduct a comprehensive analysis by a stochastic representation of the dynamics of the model, followed by ﬁnding all asymptotic stable equilibrium states of the model for any values of parameters and initial conditions.


Introduction and Related Work
A large group of epidemiological models is extensions of the Susceptible-Infected-Recovered (SIR) model [1]. These models were used for both prediction of the pandemic spread and to find optimal intervention policies for multiple types of diseases such as Polio [2], COVID-19 [3], Ebola [4], and Influenza [5]. These models use a wide range of analyses and extensions for the SIR model to properly represent the epidemiological and biological properties unique for each disease. One can use these models to obtain multiple intervention policies and properties of the pandemic dynamics [3]. For example, to predict the required number of intensive care units (ICU) to treat all severely infected individuals [6], to estimate the influence of the pandemic on the economy [3,7].
However, as the models become large, it becomes complex to numerically solve them and to obtain the model analytical properties. A specific property of interest is asymptotic stable equilibria states [8]. The way of obtaining these steps is significantly dependent on the system and its representation. For example, one can obtain the asymptotic stable equilibria states of ordinary differential equation (ODE) system numerically using a proportional-derivative controller [8]. On the other hand, one may use analytical methods which are obtaining first the equilibria states by setting the gradient of the system's state to zero and solving to find the system's state. In order to find the stability properties of such equilibria, one may use several methods such as adding error to the equilibria and show if it either decreasing or increasing under what condition [9]. Another option is to show that under given conditions the eigenvalues of the Jacobin matrix in the equilibria states are negative [2]. An additional common option is using the Lyapunov stability theorem [10][11][12][13]. While these methods are useful, they require some level of expertise to understand how to configure a specific dynamic system for each method. Furthermore, these methods require to first find the equilibria states of the model which may be a time and resource-consuming task by itself. Therefore, the analytical analysis of large-scale systems a complex task remains challenging.
One method demonstrated to be useful to model dynamics appearing in epidemics is Markov chain [14,15]. In the context of epidemics, the Markov chain method represents the dynamics, that is, the rapid spread of the disease in the population using a transmission matrix [14].
Since the pandemic spread is subject to multiple complex factors whose nature is uncertain [15], and that these factors change based on the current state of the disease spread, Markov chain is a natural approach to use to model such dynamics [16][17][18]. It was shown that Markov chain approximates the deterministic SIR model very well [19,20]. In addition, using the Markov chain method it is possible to obtain multiple analytical properties such as equilibria states and asymptotic states [21].
We propose a novel method to obtain all asymptotic stable equilibria states of an extended SIR for three age groups and for five epidemiological-states, which we develop to describe long-term immunity memory, airborne infection pandemic model [3]. Our method approximates the continuous extended SIR model using a discrete, stochastic, Markov chain representation. The paper is organized as follows. First, we introduce an extended SIR model, consists of 15 ODEs. Second, we find the asymptotic equilibrium state of the proposed model. Third, we compare the proposed method with a classical method. Finally, we discuss the main advantages and limitations of the proposed new method.

Model Definition
We describe an extended SIR epidemiological model proposed by us in [22], with the addition of the new age group (elderly) and dividing the infection-state to asymptomatic and symptomatic sub-populations. A full description of the proposed model is as follows. The model considers a constant population with a fixed number of individuals N. Each individual belongs to one of the five sub-populations: susceptible (S), infected asymptomatic (I a ), infected symptomatic (I s ), recovered (R) and dead (D) such that N = S + I s + I a + R + D such that each sub-population is non-negative. When an individual in the susceptible sub-population (S) is exposed to the infection, it is transferred to the either asymptomatic or symptomatic infected sub-population (I a , I s ) in rates β a , β s . Individuals in the symptomatic infected sub-population (I s ) stays in this sub-population on average d s I→R days, after which it is transferred to either the recovered (R) or dead (D) subpopulation. Therefore, in each time unit, some portion of infected individuals recover while others die or remain seriously ill. Individuals in the asymptomatic infected sub-population (I a ) stays in this sub-population on average d a I→R days, after which it is transferred to the recovered sub-population (R). Thus, our extended SIR model that consists of Susceptible, Infected-Asymptomatic, Infected-Symptomatic, Recovered and Deceased sub-populations is called SIIRD. A schematic view of the transition of an individual in the population between the model's states is shown in Fig. 1.
The population is divided into three classes based on their age: children, adults, and elderly because these sub-populations experience diseases in varying degrees of severity and have different infection probabilities. Individuals below age A 1 are associated with the "children" age-class while individuals below age A 2 are associated with the "adult" age-class and the complementary sub-population are associated with the "elderly" age-class. The specific threshold ages (A 1 , A 2 ) may differ in different locations but the main goal is to divide the population into three representative age classes. Since it takes A 1 years from birth to move from a child to an adult age sub-population and A 2 − A 1 from an adult to elderly, the conversion rate is set as α 1 := 1 In addition, children are born and elderly are die unrelated to the pandemic at a rate λ. We assume that different age-group spend most of their time in separation from each other which results in a relatively small rate of infected individuals infecting a susceptible individual from different group-age. Therefore, we neglecting these dynamics by setting these settings to zero. By expanding the designation to three age classes, we let S c , I a c , I s c , R c , D c , S a , I a a , I s a , R a , D a , and S e , I a e , I s e , R e , D e represent susceptible, asymptomatic infected, symptomatic infected, recovered, and dead sub-populations for children, adults, and elderly, respectively, such that {x ∈ {c, a, e} | N x := S x + I a x + I s x + R x + D x }, and N = Σ x∈{c,a,e} N x . In addition, we mark n = 15 to be the number of sub-population in the model. Afterward, in order to obtain the distribution of the sub-populations sizes in the whole population, we divide each sub-population in the overall size of the population to obtain: {x ∈ {c, a, e}, p ∈ {S, I a , I s , R, D} | p x := p x /N}.
Eqs (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15) describe the epidemic's dynamics. In Eq (1), dt is the dynamical amount of susceptible individual children over time. It is affected by the following four terms. First, with a rate β cs each symptomatic infected child infects susceptible children. Second, with a rate β ca each asymptomatic infected child infects the susceptible children. Third, children grow and pass from the children's age-class to the adult's age-class with transition rate α 1 , reduced from the children's age-class. Finally, with a rate λ children burn prospective to the death of the elderly which is not related to the pandemic.
In Eq (2), dt is the dynamical amount of symptomatic infected individual children over time. It is affected by the following four terms. First, with a rate β cs each symptomatic infected child infects the susceptible children. Second, individuals recover from the disease with a rate γ cr . Third, individuals die from the disease with a rate γ cd . Finally, children grow and pass from the children's age-class to the adult's age-class with transition rate α 1 , reduced from the adult's age-class.
In Eq (3), dt is the dynamical amount of asymptomatic infected individual children over time. It is affected by the following three terms. First, with a rate β ca each asymptomatic infected child infects the susceptible children. Second, individuals recover from the disease with a rate γ cr . Finally, children grow and pass from the children's age-class to the adult's age-class with transition rate α 1 , reduced from the adult's age-class.
In Eq (4), dt is the dynamical amount of recovered individual children over time. It is affected by the following two terms. First, in each point, a portion of the symptomatic and asymptomatic infected children recover with a rate γ cr . Second, children grow from birth and pass from the children's age-class to the adult age-class with transition rate α 1 , reduced from the children's age-class.
In Eq (5), is the dynamical amount of dead individual children over time. It is affected by the symptomatic infected children that die with a rate γ cd .
In Eq (6), dt is the dynamical amount of susceptible adult individuals over time. It is affected by the following four terms. First, children grow and pass from the children's ageclass to the adult's age-class with transition rate α 1 , added to the adult age-class. Second, adult grow and pass from the adults' age-class to the elderly age-class with transition rate α 2 , reduced from the adult's age-class. Third, with a rate β as , each symptomatic infected adult infects susceptible adults. Finally, with a rate β aa , each asymptomatic infected adult infects susceptible adults.
In Eq (7), dt is the dynamical amount of symptomatic infected individual adults over time. It is affected by the following five terms. First, children grow and pass from the children's age-class to the adult's age-class with transition rate α 1 , added to the adult age-class. Second, adult grow and pass from the adults' age-class to the elderly age-class with transition rate α 2 , reduced from the adult's age-class. Third, with a rate β as each symptomatic infected adult infects susceptible adults. Forth, individuals recover from the disease with a rate γ cr . Finally, individuals die from the disease with a rate γ ad .
In Eq (8), dt is the dynamical amount of asymptomatic infected individual adults over time. It is affected by the following four terms. First, with a rate β aa , each asymptomatic infected adult infects susceptible adults. Second, individuals recover from the disease with a rate γ ar . Third, children grow and pass from the children's age-class to the adult's age-class with transition rate α 1 , reduced from the adult's age-class. Finally, adult grow and pass from the adults' age-class to the elderly age-class with transition rate α 2 , reduced from the adult's age-class.
In Eq (9), is the dynamical amount of recovered individual adults over time. It is affected by the following three terms. First, in each point, a portion of the symptomatic and asymptomatic infected adult recover with a rate γ ar . Second, children grow from birth and pass from the children's age-class to the adult age-class with transition rate α 1 , reduced from the children's age-class. Finally, adult grow and pass from the adults' age-class to the elderly age-class with transition rate α 2 , reduced from the adult's age-class.
In Eq (10), is the dynamical amount of dead individual adults over time. It is affected by the symptomatic infected adult that dies with a rate γ ad .
In Eq (11), dt is the dynamical amount of susceptible elderly individuals over time. It is affected by the following four terms. First, adults grow and pass from the adults' age-class to the elderly age-class with transition rate α 2 , added to the elderly age-class. Second, the elderly naturally die with transition rate λ, reduced from the elderly age-class. Third, with a rate β es , each symptomatic infected elderly infects susceptible elderly. Finally, with a rate β ea , each asymptomatic infected elderly infects susceptible elderly.
In Eq (12), dt is the dynamical amount of symptomatic infected individual elderly over time. It is affected by the following five terms. First, adults grow and pass from the adult's age-class to the elderly age-class with transition rate α 2 , added to the elderly age-class. Second, the elderly naturally die with transition rate λ, reduced from the elderly age-class. Third, with a rate β es , each symptomatic infected elderly infects susceptible elderly. Forth, individuals recover from the disease with a rate γ er . Finally, individuals die from the disease with a rate γ ed .
In Eq (13), dt is the dynamical amount of asymptomatic infected individual elderly over time. It is affected by the following four terms. First, with a rate β ea each asymptomatic infected elderly infects susceptible elderly. Second, individuals recover from the disease with a rate γ er . Third, adults grow and pass from the adult's age-class to the elderly age-class with transition rate α 2 , reduced from the adult's age-class. Finally, elderly die naturally with transition rate λ, reduced from the elderly age-class.
In Eq (14), dt is the dynamical amount of recovered individual elderly over time. It is affected by the following three terms. First, in each point, a portion of the symptomatic and asymptomatic infected elderly recover with a rate γ er . Second, adult grow from birth and pass from the adult's age-class to the elderly age-class with transition rate α 2 , reduced from the children's age-class. Finally, the elderly naturally die with transition rate λ, reduced from the elderly age-class.
In Eq (15), is the dynamical amount of dead individual elderly over time. It is affected by the symptomatic infected elderly that die due to the pandemic with a rate γ ed . Therefore, the system takes the form: dI a e (t) dt = α 2 I a a (t) + β ea I a e (t)S e (t) − (λ + γ er )I a e (t), dR e (t) dt = α 2 R a (t) + γ er (I s e (t) + I a e (t)) − λR e (t), dD e (t) dt = γ ed I s e (t).
In this notation, the parameters are rates and defined the changes in the population entirely and not in the individual level. One may model the pandemic dynamics using a stochastic process due to the unstable nature of the parameters of the pandemic used in the model such as the infection rates (β cs , β ca , β as , β aa , β es , β ea ) and recovery rates (γ cr , γ cd , γ ar , γ ad , γ er , γ ed ) which differ over time as these models are affected by multiple parameters that unnecessarily taken into consideration or even immeasurable in real-world settings. Therefore, it is possible to treat these parameters as an average probability that an event would happen. Following these assumptions, one can represent the epidemiological dynamics as a transition matrix between two consecutive states of the model which represented by an n-dimensional vector, corresponding to the number of sub-populations, as follows: where h ∈ R is an arbitrary small step in time and T ∈ R n×n is the transformation matrix. The model's state at time t is defined by which therefore Eq. (18) takes the form: The transformation matrix is defined as T : where the model's parameters P (see Eq. (17)) are probabilities rather than rates as they represent the probabilities for state transfer in the individual level. The transformation matrix (T) obtained by solving dM(t) dt = ΦM(t), where dM(t) dt is taken from Eq. (16), after performing linearization on the β kl I l k S k terms to be: k ∈ {c, a, e}, l ∈ {a, s} : β kl I l k S k → ξ kl I l k , (22) such that ξ kl = β kl S k (t). Therefore, the parameter ξ kl is the probability that an infected individual would infect other individuals in the population while β kl is the probability that a suspicious would be infected by an infected individual. The parameter ξ kl changes over time as S k (t) change over time but can be treated as a constant as ξ kl is a random variable in nature, so incorporating sufficient variability to capture the dynamics of S k (t) over time.
The motivation of using this linearization is that the alternative where β kl I l k S k → β kl S k is a worse approximation. For example, consider the case where k ∈ {c, a, e}, l ∈ {a, s} : I l k = 0 ∧ S k > 0. Following the approximation where β kl I l k S k → β kl S k means that some portion of the suspicious population become infected which is impossible from an epidemiological perspective. On the other hand, following the linearization in Eq. (22), will result in k ∈ {c, a, e}, l ∈ {a, s} : I l k = 0 in this scenario. Therefore, matrix T is the stochastic, linear, approximation of the transformation between two states of the ODE-based model (see Eq. (16)). Nevertheless, the models described in Eq. (16) (ODE, deterministic model) and Eqs. (20-21) (linear transformation matrix, stochastic model) are analytically differ since in Eq. (16) the parameters P can be assigned any real value. While it may no longer describe epidemiological dynamics, the mathematical model is well-defined in such scenario. On the other hand, Eqs. (20)(21) required the parameters to be ∀p ∈ P : p ∈ (0, 1] this is according to Eq. (21) which is stochastic matrix and therefore satisfies that each row sums to 1. This condition not met if ∀p ∈ P : p ∈ (0, 1] not hold. Therefore, the model represented by Eq. (16) inclusive the model represented by Eqs. (20)(21).
However, for the sub-space where both models are define, they are numerically equal for any finite time. Indeed, for a given norm function || · || : R n → R, start condition M(0), and time interval [0, t max ]. The state of the stochastic SIIRD (Eqs. (20-21)) M s (t) and the state of the deterministic SIIRD (Eq. (16)) M d (t) satisfies: where the parameters {c, a, e}, l ∈ {a, s} : ξ kl (t) = β kl S k (t) are updated at each point in time t.
By approximating the deterministic representation (Eq. (16)) system using the forward Euler method [23] in each one of the states, the approximation introduce O(h 2 ) error for each step in time. Now, one needs to take t max h steps in time to cover [0, t max ] which introduce an overall t max h · O(h 2 ) = O(t max h) error. Therefore: ||M d − M s || < t max h. As a result, for h < /t max , the condition ∀t ∈ [0, t max ] : ||(M s (t) − M d (t)|| < satisfied. Afterward, we define a stochastic process of the dynamics in Eqs. (20)(21) at each point at time (t), for each sub-population M i (t) ∈ M(t), there are three possible options for each individual in the population in respect to this sub-population. First, an individual can be transform from M i (t) to M j (t + 1) (i = j ∈ [1, . . . , n]) in a probability α which results in Φ i,j = α. Second, an individual can transform from M j (t) to M i (t + 1) in a probability ξ which results in Φ j,i = ξ in a symmetric way to the first case. Third, an individual in M i (t) can stay in M i (t + 1). This is a default case and happens in probability 1 − Σ n j=1 (Φ i,j ) as the complementary probability to all the probabilities of an individual to transform from M i (t). These are the only options possible for an individual in each sub-population M i (t) as ∀t : S c (t + h) = (1 − α 1 )S c (t) − ξ ca I s c (t) − ξ cs I a c (t) + λ(S e (t) + I s e (t) + I a e (t) + R e (t)), R e (t + h) = (1 − λ)R e (t) + α 2 R a (t) + γ ar (I s e (t) + I a e (t)), D e (t + h) = D e (t) + γ ad I s e (t).
The representation in Eq. (23) is isomorphic to the one in Eqs. (20)(21). However, Eq. (23) treats the dynamic as a stochastic process in nature rather than approximating the deterministic ODE-based dynamics while restoring the underline behavior of the epidemiological system.

Asymptotic Stable Equilibria States
In epidemiology, there are two types of cases that broadly interesting decision-makers. First, the state of the population in a long-range after the end of the pandemic. Second, the equilibria points and their stable or unstable nature.
The state of the population in a long range after the end of the pandemic can be mapped to the asymptotic state of the pandemic in time as after long enough time (e.g. t → ∞), the population is either survive the pandemic and restore regular dynamics without the pandemic or extinct. While the second scenario is trivial as the population is distributed between the different death states of the model, the first scenario holds a larger space of options. Specifically, the pandemic can die out (i.e. the size of the infected population is zero) and as a result after a few generations only susceptible individuals would remain. On the other hand, in some settings the pandemic may not die out but kept under control such that the pandemic converges to a steady-state.
The equilibria states are important for decision-makers as these hold a promise of a scenario that remains the same unless some action is taken or a major event takes place. However, the equilibria states should be divided into two groups. On the one hand, the unstable equilibria states which provide some level of stability are still problematic due to their unstable nature where even relatively small change results in drastic outcomes. On the other hand, stable equilibria do not have this issue. Therefore, in this section, we analyze the model's asymptotic equilibria states. One may try to obtain the asymptotic equilibria states and their stability properties from the ODE-based representation (e.g. Eq. (16)). Nevertheless, this approach would require solving a system of n-dimensional, non-linear, heterogeneous, ODE system which is both numerically and analytically complex and time-consuming task. On the other hand, by defining a non-homogeneous discrete-time Markov chain represented by the transformation function (Eq. (23)) with state space M(t) (Eq. (18)), in respect to the states in Eq. (16) one can find the asymptotic equilibrium as follows. First, in order to model the dynamics as a Markov chain one needs to show that where {M} t is a stochastic process with values in the state space for all t ≥ 0 and all states i 0 , ...i t , j and h is an arbitrary small step in time [24]. T satisfies Eq (24) if any value in M(t + h) depended only on the previous state M(t). Indeed, as T does not depend on t or any value of M(t) or previous state, the condition is satisfied. Therefore, we show Eq. (23) describes a Markovian process. As a result, given the model's initial condition (M(0)), the model's state at some time t is defined by M(t) = T t M(0) [21]. Now, assume any initial condition M(0). From Eq. (23) and Fig. (1), it is possible to notice a few sub-processes in the dynamics. First, an individual that at some time t reaches a death state (corresponding to lines 5, 10, and 15 in Eq. (23)), stays there as the coefficients of D c , D a , and D e is 1 for any parameter's values. Second, if the pandemic ended -namely, I s c + I a c + I s a + I a a + I s e + I a e = 0, there are two possible cases -the population is extended or some portion of the population (or even the whole population) is survived. In the case the population is extended, the obtained state is a distribution over the {D c , D a , D e } states while the other sub-populations are 0. While in the second option, the population is distributed over the {S c , S a , S e , R c , R a , R e , D c , D a , D e } as the infection states are 0. However, after t > 1/α 1 all the individuals at R c transform to R a . Similarly, after time t > 1/α 2 all the individuals at R a transform to R e and finally after t > λ all the individuals at R e transform to S c . As a result, after time t > 1/α 1 + 1/α 2 + 1/λ the sub-populations R c = R a = R e = 0. While in the same time the remain population at S c , S a , and S e circulate between these states. Finally, in the case where the pandemic is not over, after some time t the pandemic will do end as ∀t : T t which means these sub-populations are eventually decreased over time. As a result, for any start condition M(0) and parameters ∀p ∈ P : p ∈ (0, 1], at t → ∞ the model's asymptotic state (lim t→∞ M(t)) takes the form: where N is the total size of the population as defined in Section 2.
The values {ν} 6 i=1 are dependent on the initial conditions and the model's parameters and therefore complex and time-consuming to find. Therefore, Eq. point on. Therefore, one can take advantage of this property in order to obtain the asymptotic stable equilibria as they necessarily follow Eq. (25). In order to obtain the asymptotic equilibrium state, we set Eq. (25) in Eq. (18) and obtain: It is possible to divide the types of equilibria into two sub-groups, where ν 1 = ν 3 = ν 5 = 0 and otherwise. The first option corresponding to the scenario in which the population is extended due to the pandemic. By setting ν 1 = ν 3 = ν 5 = 0 in Eq. (27), one obtain that all combinations of {ν 2 , ν 4 , ν 6 } such that ν 2 + ν 4 + ν 6 = N are equilibrium. Resulting in (N + 1) 2 for the population in size N as its combinatory equivalent to dividing N items into three (allowing empty) groups. On the other hand, assuming ν 1 = 0, ν 3 = 0, ν 5 = 0 the state is equilibrium if and only if: which means the asymptotic state is also equilibrium if the following condition is fulfilled: As a result, for any initial condition and the model's parameters values, there is a time t 0 such that for all t > t 0 , Eq. (25) is fulfilled. In the case where either ν 1 = ν 3 = ν 5 = 0 or α 1 ν 1 = α 2 ν 3 = λν 5 the model is in asymptotically stable equilibrium state.

Comparison With Classical Methods
We compare the proposed method shown in Section 3 with a classical method used to obtain equilibria states and their stability properties for dynamical systems [25,26].

Equilibrium
The equilibria states of the model (Eq. (16)) are defined as the states of the model in which the gradient is zero. Therefore, Eq. −(β cs I s c + β ca I a c + α 1 )S c + λ(S e + I s e + I a e + R e ) = 0, β cs I s c S c − (α 1 + γ cr + γ cd )I s c = 0, β ca I a c S c − (α 1 + γ cr )I a c = 0, γ cr (I s c + I a c ) − α 1 R c = 0, γ cd I s c = 0, α 1 S c − (α 2 + β as I s a + β aa I a a )S a = 0, α 1 I s c + β as I s a S a − (α 2 + γ ar + γ ad )I s a = 0, α 1 I a c + β aa I a a S a − (α 2 + γ ar )I a a = 0, α 1 R c + γ ar (I s a + I a a ) − α 2 R a = 0, γ ad I s a = 0, α 2 S a − (λ + β es I s e (t) + β ea I a e )S e = 0, α 2 I s a + β es I s e S e − (λ + γ er + γ ed )I s e = 0, α 2 I a a + β ea I a e S e − (λ + γ er )I a e = 0, α 2 R a + γ er (I s e + I a e ) − λR e = 0, γ ed I s e = 0. (29) One can notice that from Eqs. (16,29) it follows that where ν 2 , ν 4 , ν 6 are arbitrary constants such that {ν 2i } 3 i=1 ≥ 0 and ∑ 3 i=1 ν 2i ≤ N, because if either I s c , I s a , or I s e is not equal zero then the gradient of D c , D a , or D e is not zero and therefore the state is not equilibrium by definition. Therefore, I s c = 0, I s a = 0, I s e = 0 that leads to D c = ν 2 , D a = ν 4 , D e = ν 6 due to Eqs. (5, 10, and 15). As a result, one left with −(β ca I a c + α 1 )S c + λ(S e + I a e + R e ) = 0, β ca I a c S c − (α 1 + γ cr )I a c = 0, γ cr I a c − α 1 R c = 0, α 1 S c − (α 2 + β aa I a a )S a = 0, α 1 I a c + β aa I a a S a − (α 2 + γ ar )I a a = 0, α 1 R c + γ ar I a a − α 2 R a = 0, α 2 S a − (λ + β ea I a e )S e = 0, α 2 I a a + β ea I a e S e − (λ + γ er )I a e = 0, α 2 R a + γ er I a e − λR e = 0.
It is easy to check that for the considered system (1) the Jacobian matrix J(x * ) coincides with the matrix Φ given in Eq. (21). A state is an asymptotic stable if and only if the matrix is negative by defined for any values of the parameters (Eq. (17)). Unfortunately, the matrix Φ is neither diagonal nor triangle and therefore, one is able to determinate if it is negative definite or not by analytically obtaining the determine and investigate its properties. However, this will result in a 15-ordered polynomial which is much more time and resource consuming compared to the presented Section 3.

Discussion
We propose a novel method to analytically obtain all asymptotic stable equilibria states. We present this method for an extended SIR model, for the three age groups SIIRD model. This method is based on the Markov chain model whose parameters are deterministic (Eq. (23)). Using this representation, one is able to obtain all asymptotic stable equilibrium states of the model for any given start condition and properties using the stationary state (Eq. (31)).
When comparing the proposed method with classical methods of obtaining equilibria and its stability properties, it is clear that for large-scale SIR models that the proposed method is superior due to several benefits. First, the classic method requires some level of algebraic expertly to unstitch the equations describing the dynamics while the proposed method treats them as a single process and therefore redundant this process at all. Second, the classic method by itself is not able to identify all asymptotic stable equilibria by itself as one requires to manually find all equilibria states and investigate each one independently which is time and resource consuming. On the other hand, the proposed model obtains analytically all asymptotic stable equilibria as it finds the stationary state of the stochastic process that represents the dynamics.
Naturally, converting the deterministic biological related rate coefficients such as the recovery rates γ or infection rates ξ into transformation probabilities may create cases that are not corresponding to the biological dynamics on the individual level and therefore should be treated with care. For example, a susceptible individual (p ∈ S a ) can be infected and transformed into the asymptomatic infected sub-population (I a a ) in a given time t. Immediately afterward, in time t + 1, there is a chance γ ar that the same individual is recovered and transformed to the recovered sub-population (R a ).
The stochastic representation of pandemic dynamics allows more flexibility and credibility than treating model parameters as deterministic values due to the fact that data often involves uncertainty [15]. This approach allows simulating pandemic dynamics that based on an extended SIR model using distributed systems models which allows adding upon the epidemic dynamics additional social [3], non-pharmaceutical and pharmaceutical intervention (NPI / PI) policies [27], and economical policies [7].