Θ-SEIHRD mathematical model of Covid19-stability analysis using fast-slow decomposition

Introduction

The coronavirus belongs to the severe acute respiratory syndrome (SARS) family. It usually does not cause disease in humans, but infects animals—mammals and poultry. If humans are infected with the virus, the disease usually causes mild cooling, and the infection passes without any treatment. Conversely, if the infected person has a weak immune system for some reason, the coronavirus can be fatal. It is a common virus that attacks almost every person at least once in his/her lifetime—especially during childhood — and, as is usually the case, it is not dangerous. However, the current epidemic is a recently mutated virus that has become fatal (which is why the World Health Organization initially called it “the new coronavirus” 2019). It should be noted that, in 2002–2003, there was an outbreak of a similar virus from the SARS family, causing the death of more than 800 people (General Health Fund in Israel, https://www.clalit.co.il/he/yourhealth/family/Pages/coronavirus.aspx). The virus passes from person to person with a drop infection—as in the case of influenza: droplets from an infected person carry the virus into the air and penetrate the respiratory system of uninfected humans. The virus can also be passed through contact with an infected person. Hence, the importance of frequent hand washing with water and soap or alcohol is propagated. This is also why shaking hands and touching the face are to be avoided. It should be noted that the infected person can be either symptomatic (has the symptoms of the disease), pre-symptomatic (infected but has yet to develop symptoms), or asymptomatic (infected with the virus, but has no symptoms and will not develop the disease). By now, we already know that a significant proportion of those infected are asymptomatic and, therefore, without being aware of it, can endanger those who are not infected. Hence, the directive that everyone, without exception, should adhere to the guidelines on maintaining social distancing. In light of the current situation, and in the absence of a scientifically proven drug or vaccine, scientists from all over the world have been trying to find an effective vaccine for this virus. They include researchers from various fields of science such as biology, chemistry, physics and even mathematics. Researchers have developed various mathematical models that describe the evolution of the virus in the population (Liu et al., 2020; Ivorra & Ramos, 2020a, 2020b; Yan & Cao, 2019; Kucharski et al., 2020; Li et al., 2020; Ivorra et al., 2014; Johns Hopkins University (JHU), 2020). All countries whose population was infected with the virus have tried to prevent its spread by imposing a complete lockdown of the population. This effort is an attempt to contain the occurrence of new infections. The mathematical description for this situation in which there are no new infections, which means that the spread has stopped, is called equilibrium. Therefore, the optimal state is equilibrium stability. In this paper, we investigate the stability of the θ-SEIHRD model, a mathematical model proposed by Ivorra et al. (2020).

In general, given a large mathematical model that contains a system of linear or nonlinear differential equations, it is difficult, and sometimes impossible, to analyze the model analytically. Even numerical analysis is sometimes complicated and requires a very expensive computer running time. Therefore, the need to reduce the number of equations (i.e., reduced model) while maintaining the physical/biological dynamics of the system is required. In this paper, we introduce and apply a method called singular perturbed vector field (SPVF) (Bykov, Goldfarb & Gol’dshtein, 2006; Nave, 2017), which is a new version of the ILDM method. Because a large number of differential equations (a mathematical model) are presented in a hidden hierarchy form, there is no explicit time scale of the system. Hence, our aim is to first expose the time scale of a given system and then apply a reduction method. This is the main result of the SPVF method. The SPVF method transfers the original system to the form of the singularly perturbed system (SPS), that is, a system with an explicit hierarchy of the dynamic variables of the model. Once we transfer the system to an SPS one, the new system can be treated by the highly effective standard SPS theory for model reduction and decomposition, as we mentioned above, without losing the essential dynamics of the original system.

The rest of the paper is organized as follows. In “The θ-SEIHRD Model of COVID-19”, we describe the mathematical model of the coronavirus called θ-SEIHRD model. In “Results and Analysis”, we present the results of our comparative analysis of the cases reported by the Chinese authorities. In addition, we apply the SPVF model to the θ-SEIHRD model and find the stability of the equilibrium points. “Discussion” presents our discussion. Finally, “Conclusions” presents the conclusion of this paper.

The θ-SEIHRD Model of COVID-19

In this section, we introduce the mathematical model of the θ-SEIHRD model as presented in (Ivorra et al., 2020). This model is based on the Be−CoDiS model presented in (Ivorra, Ramos & Ngom, 2015). It is important at this point to note that this model is not the SIR-SEIR standard model. It considers the known special characteristics of the considered disease, such as the existence of infectious undetected cases and the sanitary conditions and infectiousness of hospitalized patients. The assumptions of the model can be found in (Ivorra et al., 2020). Here, we present the mathematical formulation of the COVID-19 spread. The system of equations includes nine nonlinear ordinary differential equations and has the following forms: (1) $$\begin{array}{ll}{\displaystyle \frac{dS}{dt}}& =-{\displaystyle \frac{S}{N}\left(m_{{}_E}{\mathrm{\beta }}_{E}E+m_{{}_I}{\mathrm{\beta }}_{I}I+m_{{}_{I_u}}{\mathrm{\beta }}_{I_u}{I}_u\right)}-{\displaystyle \frac{S}{N}\left(m_{{}_{H_R}}{\mathrm{\beta }}_{H_R}{H}_R+m_{{}_{H_D}}{\mathrm{\beta }}_{H_D}{H}_D\right)}\\ & -{\mathrm{\mu }}_{m}S+{\mathrm{\mu }}_n\left(S+E+I+I_u+R_d+R_u\right)\equiv F_S(\overrightarrow{V})\end{array}$$ (2) $$\begin{array}{ll}{\displaystyle \frac{dE}{dt}}& ={\displaystyle \frac{S}{N}\left(m_{{}_E}{\mathrm{\beta }}_{E}E+m_{{}_I}{\mathrm{\beta }}_{I}I+m_{{}_{I_u}}{\mathrm{\beta }}_{I_u}{I}_u\right)}+{\displaystyle \frac{S}{N}\left(m_{{}_{H_R}}{\mathrm{\beta }}_{H_R}{H}_R+m_{{}_{H_D}}{\mathrm{\beta }}_{H_D}{H}_D\right)}\\ & -{\mathrm{\mu }}_{m}E-{\mathrm{\gamma }}_{{}_E}E+{\mathrm{\tau }}_1-{\mathrm{\tau }}_2\equiv F_E(\overrightarrow{V})\end{array}$$ (3) $${\displaystyle \frac{dI}{dt}={\mathrm{\gamma }}_{{}_E}E-\left({\mathrm{\mu }}_m-{\mathrm{\gamma }}_{{}_I}\right)I\equiv F_I(\overrightarrow{V})}$$ (4) $${\displaystyle \frac{d{I}_u}{dt}=\left(1-\mathrm{\theta }\right){\mathrm{\gamma }}_{{}_I}I-\left({\mathrm{\mu }}_m+{\mathrm{\gamma }}_{{}_{I_u}}\right)I_u\equiv F_{I_u}(\overrightarrow{V})}$$ (5) $${\displaystyle \frac{d{H}_R}{dt}=\mathrm{\theta }\left(1-{\displaystyle \frac{\mathrm{\omega }(t)}{\mathrm{\theta }}}\right){\mathrm{\gamma }}_{{}_I}I-{\mathrm{\gamma }}_{{}_{H_R}}{H}_R\equiv F_{H_R}(\overrightarrow{V})}$$ (6) $${\displaystyle \frac{d{H}_D}{dt}=\mathrm{\omega }(t){\mathrm{\gamma }}_{{}_I}I-{\mathrm{\gamma }}_{{}_{H_D}}{H}_D\equiv F_{H_D}(\overrightarrow{V})}$$ (7) $${\displaystyle \frac{d{R}_d}{dt}={\mathrm{\gamma }}_{{}_{H_R}}{H}_R-{\mathrm{\mu }}_m{R}_d\equiv F_{R_d}(\overrightarrow{V})}$$ (8) $${\displaystyle \frac{d{R}_u}{dt}={\mathrm{\gamma }}_{{}_{I_u}}{I}_u-{\mathrm{\mu }}_m{R}_u\equiv F_{R_u}(\overrightarrow{V})}$$ (9) $${\displaystyle \frac{dD}{dt}={\mathrm{\gamma }}_{{}_{H_D}}{H}_d\equiv F_D(\overrightarrow{V})}$$where $\overrightarrow{V}=(S,E,I,I_u,H_R,H_D,R_d,R_u,D)$; ${\mathrm{\gamma }}_E,{\mathrm{\gamma }}_I,{\mathrm{\gamma }}_{H_R}$and ${\mathrm{\gamma }}_{I_u}$ are the transition rates of E, I, H_R, and I_u, respectively $[{\displaystyle \frac{1}{\mathrm{d}\mathrm{a}\mathrm{y}}]}$; β(·) is the disease contact rate $[{\displaystyle \frac{1}{\mathrm{d}\mathrm{a}\mathrm{y}}]}$; μ is the fatality rate $[{\displaystyle \frac{1}{\mathrm{d}\mathrm{a}\mathrm{y}}]}$; and τ indicates the infected people, who move from one territory to another per day. The initial conditions of the system at time t₀ (for each country, however, the initial time can be changed. For example, in Wuhan, China, t₀ = 10 A. CET (technical adjustment) because it was on 7 December 2019 when the first case of patients with symptoms was confirmed) are as follows: (10) $$\begin{array}{l}S(t_0)=S_0,E(t_0)=E_0,I(t_0)=I_0,I_u(t_0)=I_{u0},H_R(t_0)=H_{R0},\\ H_D(t_0)=H_{D0},R_d(t_0)=R_{d0},R_u(t_0)=R_{u0},D(t_0)=D_0\end{array}$$

As can be seen above, the presentation of the above mathematical model has no explicit hierarchy and, hence, it is impossible to apply different asymptotic methods. In the next section, we apply the SPVF methods to expose the hierarchy of the considered model.

Results and Analysis

In this section, we apply the SPVF method to the system of Eqs. (1)–(9), as presented in (Nave, 2017). While we apply the SPVF method, we expose the hierarchy of the model. This procedure enables us to rewrite the model under consideration in new coordinates and, hence, the “new” model can be decomposed into the so-called fast and slow subsystems.

The set of parameters and functions used in our calculations is as follows:

Functions:

(12) $$\begin{array}{l}\begin{array}{l}{\mathrm{\beta }}_{I_u}=0.375+{\displaystyle \frac{0.375}{(1-\mathrm{\omega }(t))}(1-\mathrm{\theta }),}\\ \mathrm{\omega }(t)=0.003\cdot m_I+0.997\cdot (1-m_I),\\ m_{{}_E}=m_{{}_{H_R}}=m_{{}_{H_D}}=m_{{}_{I_u}}=m_I=e^{-0.01(t-23)},\\ {\mathrm{\gamma }}_{{}_{H_D}}={\displaystyle \frac{1}{d_{I_u}+6(1-m_{{}_I})+{\mathrm{\delta }}_R}}\end{array}\end{array}$$

Here, we present only the absolute values of the eigenvalues of the SPVF method (and not the eigenvectors) to avoid overwhelming the readers with too much information: (13) $$\begin{array}{l}\begin{array}{l}{\mathrm{\lambda }}_1=2.1356\cdot {10}^{25}\\ {\mathrm{\lambda }}_2=7.8618\cdot {10}^{16}\\ {\mathrm{\lambda }}_3=5.3590\cdot {10}^{16}\\ {\mathrm{\lambda }}_4=3.4835\cdot {10}^{15}\\ {\mathrm{\lambda }}_5=1.0870\cdot {10}^{15}\\ {\mathrm{\lambda }}_6=1.0870\cdot {10}^{15}\\ {\mathrm{\lambda }}_7=5.4587\cdot {10}^{14}\\ {\mathrm{\lambda }}_8=5.4586\cdot {10}^{14}\\ {\mathrm{\lambda }}_9=3.8120\cdot {10}^{13}\end{array}\end{array}$$As we observed from the above results, and according to the SPVF method, the maximal gap of the eigenvalues is between λ₁ and λ₂. This gap implies that the first dynamic variable of the system, which is written in the new coordinates, is fast compared with the rest of the variables.

The results of the SPVF algorithm are presents in Figs. 1 and 2. We compute the value of |λ_{i + 1}|/|λ_i| (i = 1,…,9), for 365 days and plot only the maximum values of this quotient for every day. According to Fig. 1 we can see that for every day the gap indeed exists (for every day the graph is not zero). This means that the SPVF algorithm exposes the hierarchy of the system at the new coordinates. But it is not clear from these results what is the fast subsystem and what is the slow one. Hence, in order to know what is the fast subsystem and what is the slow one we match every eigenvalue to its corresponding equation at the model such that the eigenvalue λ_i belongs to the i equation at the new coordinates. Once we did this matching and with the results of the maximal gap (present at Fig. 1) we can know what is the fast subsystem and what is the slow one. For example, if the maximal gap is |λ₆|/|λ₇| then the fast subsystem is Eqs. (1)–(6) and the slow subsystem are Eqs. (7)–(9). The results are present in Fig. 2. We can see that at some days the fast subsystem includes only the first equation and for some other days the fast subsystem includes the first two equations of the model in the new coordinates. One can analyze the results of the SPVF algorithm from the other direction. That is, if we look at Fig. 2 on the 100th day, for example, then we see that only the first equation is the fast subsystem. But if we look at the 250th day for example, then the first two equations are the fast subsystem.

The SPVF algorithm aims to find a frame of coordinates for which the system has an SPS form. Hence, we have a degree of freedom to choose any frame of coordinates for this purpose. At our analysis, we have chosen the frame of coordinates, that is, eigenvectors, that belong to the 49th day where on this day we received the maximal gap of all the gaps that we obtained from the SPVF algorithm. On this day the fast subsystem is the first equation of the model.

We denote the new variables of the model in the new coordinates as (14) $$\begin{array}{l}\tilde{\overrightarrow{V}}=(\tilde{S},\tilde{E},\tilde{I},{\tilde{I}}_u,{\tilde{H}}_R,{\tilde{H}}_D,{\tilde{R}}_d,{\tilde{R}}_u,\tilde{D})\end{array}$$

The new model is rewritten in the new coordinates using the eigenvectors ${\overrightarrow{u}}_1,{\overrightarrow{u}}_2,...,{\overrightarrow{u}}_9$, which correspond to the eigenvalues λ₁, λ₂,…,λ₉. This means that the new dynamic variables of the model are linear combinations of its old variables when the coefficients of the linear combinations are taken from the eigenvectors, that is, (15) $$\left(\begin{array}{c}\tilde{S}\\ \tilde{E}\\ \tilde{I}\\ {\tilde{I}}_u\\ {\tilde{H}}_R\\ {\tilde{H}}_D\\ {\tilde{R}}_d\\ {\tilde{R}}_u\\ \tilde{D}\end{array}\right)=\left(\begin{array}{ccccccccc}⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ & & & & & & & & \\ & & & & & & & & \\ {\overrightarrow{u}}_1^t& {\overrightarrow{u}}_2^t& {\overrightarrow{u}}_3^t& {\overrightarrow{u}}_4^t& {\overrightarrow{u}}_5^t& {\overrightarrow{u}}_6^t& {\overrightarrow{u}}_7^t& {\overrightarrow{u}}_8^t& {\overrightarrow{u}}_9^t\\ & & & & & & & & \\ & & & & & & & & \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\end{array}\right)\cdot \left(\begin{array}{c}S\\ E\\ I\\ I_u\\ H_R\\ H_D\\ R_d\\ R_u\\ D\end{array}\right)$$where t refers to the transpose operation. In matrix form, system (15) can be written as (16) $$\tilde{\overrightarrow{V}}=\mathcal{B}\overrightarrow{V}$$

Here, we denote the matrix $\mathcal{B}$ as a matrix whose columns are the eigenvectors of the SPVF method. The right-hand side (RHD) of this system is a function of the variables $\overrightarrow{V}$, whereas its left-hand side is a function of the new variables $\tilde{\overrightarrow{V}}$. To write the system of Eq. (16), with the same variables $\tilde{\overrightarrow{V}}$, we should differentiate this system with respect to time, that is, (17) $$\dot{\tilde{\overrightarrow{V}}}=\mathcal{B}\dot{\overrightarrow{V}}=\mathcal{B}{\overrightarrow{F}}_{\overrightarrow{V}}(\overrightarrow{V})$$where (18) $${\overrightarrow{F}}_{\overrightarrow{V}}\left(\overrightarrow{V}\right)=(F_S,F_E,F_I,F_{I_u},F_{H_R},F_{H_D},F_{R_d},F_{R_u},F_D)$$

is the vector field of systems (1)–(9). To rewrite both sides of the model using the new variables, we should express the old variables as a function of the new variables using the inverse matrix of the eigenvectors, that is, (19) $${\mathcal{B}}^{-1}\tilde{\overrightarrow{V}}=\overrightarrow{V}$$

Substituting Eqs. (19) into (17), we obtain the mathematical model in the new coordinates with the initial conditions as (20) $$\begin{array}{l}\begin{array}{l}\dot{\tilde{\overrightarrow{V}}}=\mathcal{B}{\overrightarrow{F}}_{\overrightarrow{V}}({\mathcal{B}}^{-1}\overrightarrow{V})\equiv H_{\tilde{\overrightarrow{V}}}(\tilde{\overrightarrow{V}})\\ \tilde{\overrightarrow{V}}(0)=\mathcal{B}\overrightarrow{V}(0)\end{array}\end{array}$$where (21) $$H_{\tilde{\overrightarrow{V}}}\left(\tilde{\overrightarrow{V}}\right)=(H_{\tilde{S}},H_{\tilde{E}},H_{\tilde{I}},H_{{\tilde{I}}_u},H_{{\tilde{H}}_R},H_{{\tilde{H}}_D},H_{{\tilde{R}}_d},H_{{\tilde{R}}_u},H_{\tilde{D}})$$is the vector field of the new model (the RHD of the model is written in the new coordinates).

Stability analysis

After we transformed and presented the model in the new coordinates using the eigenvectors of the SPVF method, the model can be decomposed into the fast and slow subsystems based on the gap of the eigenvalues. The first eigenvalue, λ₁, indicate that the first variable of the system, say $\tilde{S}$, is the fast variable of the system and that $\tilde{E}$, $\tilde{I}$, ${\tilde{I}}_u$, ${\tilde{H}}_R$, ${\tilde{H}}_D$, ${\tilde{R}}_d$, ${\tilde{R}}_u$ and $\tilde{D}$ are the slow ones. Hence, the model presented in the system of Eq. (20) can be decomposed as follows: $$\begin{array}{l}\mathrm{f}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{s}\mathrm{y}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{m}\{\begin{array}{l}\frac{d\tilde{S}}{\frac{}{dt}}=H_{\tilde{S}}(\tilde{\overrightarrow{V}})\end{array}\end{array}$$ $$\begin{array}{l}\mathrm{s}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{s}\mathrm{y}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{m}\{\begin{array}{l}\frac{d\tilde{E}}{dt}=H_{\tilde{E}}\left(\tilde{\overrightarrow{V}}\right)\\ \frac{d\tilde{I}}{dt}=H_{\tilde{I}}\left(\tilde{\overrightarrow{V}}\right)\\ \frac{d{\tilde{I}}_u}{dt}=H_{{\tilde{I}}_u}\left(\tilde{\overrightarrow{V}}\right)\\ \frac{d{\tilde{H}}_R}{dt}=H_{{\tilde{H}}_R}\left(\tilde{\overrightarrow{V}}\right)\\ \frac{{\tilde{H}}_D}{dt}=H_{{\tilde{H}}_D}\left(\tilde{\overrightarrow{V}}\right)\\ \frac{{\tilde{R}}_d}{dt}=H_{{\tilde{R}}_d}\left(\tilde{\overrightarrow{V}}\right)\\ \frac{d{\tilde{R}}_u}{dt}=H_{{\tilde{R}}_u}\left(\tilde{\overrightarrow{V}}\right)\\ \frac{d\tilde{D}}{dt}=H_{\tilde{D}}\left(\tilde{\overrightarrow{V}}\right)\end{array}\end{array}$$

According to the SPVF method, the stability analysis can be investigated on the fast subsystem. To find the equilibrium points of the model, we should solve the following linear system: (22) $$H_{\tilde{S}}\left({\tilde{S}}^{\ast },{\tilde{E}}_0,{\tilde{I}}_0,{\tilde{I}}_{u0},{\tilde{H}}_{R0},{\tilde{H}}_{D0},{\tilde{R}}_{d0},{\tilde{R}}_{u0}{\tilde{D}}_0\right)=0$$

This implies that we look for the equilibrium points of the fast variable (denoted by an asterisk), that is, ${\tilde{S}}^{\ast }$, while the rest of the variables are frozen, and we take their values as the initial conditions of the system. System (22) is a system of one equation with one unknown (fast) variables: ${\tilde{S}}^{\ast }$. We solve this system and substitute the values ${\tilde{S}}^{\ast }$ in the slow subsystem and solve the system for the slow variables, that is, ${\tilde{E}}^{\ast }$, ${\tilde{I}}^{\ast }$, ${\tilde{I}}_u^{\ast }$, ${\tilde{H}}_R^{\ast }$, ${\tilde{H}}_D^{\ast }$, ${\tilde{R}}_d^{\ast }$, ${\tilde{R}}_u^{\ast }$ and ${\tilde{D}}^{\ast }$. Now, we have eight equations with eight unknown (slow) variables.

To check the stability of these equilibrium points, we compute the Jacobian matrix calculated at the equilibrium point. It is important at this point to note that the equilibrium points in the new coordinates have no biological meaning. Moreover, to provide a biological interpretation to the equilibrium points that we received, we have to transfer them to the old coordinates for which we need to calculate the inverse matrix of the matrix with the eigenvectors, that is, (23) $${\overrightarrow{V}}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}}^{\ast }={\mathcal{B}}^{-1}{\tilde{\overrightarrow{V}}}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}}^{\ast }$$

Here, we inverse transform only the stable equilibrium points. We received only one stable point in the new coordinates, and, therefore, we should transform only one stable equilibrium point.

According to our numerical simulations, if the initial condition is on 12 January 2020, then the equilibrium point is reached after 49 days. According to official reports (from the Wuhan County authorities), the closure of Wuhan ended after 62 days. This means that the situation stabilized after 62 days. As we can see from the results of the mathematical model and what actually happened, there is a difference of 13 days. Medically, another 13 days of lockdown would not have affected the civilians. However, apparently, the economic considerations outweighed the medical considerations and the authorities decided to lift the lockdown earlier than expected. In Fig. 3, we present the evolution of the number of cases in China for different values of the parameters. As we can see from these results, the equilibrium point is attained at approximately 83K. The results are summarized in Table 1. In this table, we present only the important variables that are stable.

Table 1:

Equilibrium point.

Equilibrium point of the model vs. Reporting from the authorities.

	Equilibrium point of the model	Reporting from the authorities
Coronavirus Cases Deaths	82,921 4,597	82,992 4,634

DOI: 10.7717/peerj.10019/table-1

Discussion

As we have shown in the previous section, we obtain the stable equilibrium points of the mathematical model owing to the application of the SPVF method. The application of this method enables us to decompose the system into a fast and a slow subsystems. After the decomposition, we explored only the fast subsystem. In general, instead of decomposing a given system into fast and slow subsystems, one can solve the original system and determine the stable equilibrium points numerically. However, the big problem with the numerical method is that the equilibrium points are represented by their values and not as points, and they depend on the original system’s parameters as can be obtained by the considered decomposition. That is, if we want to change the system parameters and determine new equilibrium points, we have to resolve the mathematical model numerically each time, which is time consuming. Moreover, the SPVF method allows us to first find all the equilibrium points of the original model analytically. In addition, the equilibrium points depend on the parameters of the original system. Therefore, if we want to change the model parameters, we do not have to resolve the mathematical model again but only need to change the parameters in the stable equilibrium points that have been determined. As can be seen from the results obtained from the mathematical model and from the results reported by the authorities, the relative error of the model can be calculated for which, we define the relative error of the equilibrium points for each dynamic variable of the system as (24) $$E_{(\cdot )}={\displaystyle \frac{|U_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}-U_{\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{l}}|}{U_{\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{l}}}\cdot 100\mathrm{\%}}$$where U_model is the dynamic variable of the system, that is, the stable equilibrium point of the original model (after the inverse transform of the stable equilibrium points) and U_real denotes the values reported by the authorities. The results are as follows: (25) $$\begin{array}{l}{E}_S=0.134\mathrm{\%},E_E=2.349\mathrm{\%},E_I=1.169\mathrm{\%},E_{I_u}=3.276\mathrm{\%},E_{H_R}=1.155\mathrm{\%},\\ E_{H_D}=4.633\mathrm{\%},E_{R_D}=0.122\mathrm{\%},E_{R_u}=4.091\mathrm{\%},E_D=1.445\mathrm{\%}\end{array}$$

As can be seen from the results of the relative error, the reported cases are indeed small in number, as was expected from the decomposition method.

Conclusions

In this paper, we investigated the stability of the θ-SEIHRD mathematical model of COVID-19. The mathematical model of coronavirus is presented with a hidden hierarchy. This implies we cannot know which of the dynamic variables of the model is progressing fast and which one is progressing slowly. The hidden hierarchy of the model neither allows for an asymptotic analysis in general, nor an analytical investigation in particular, but only the running of numerical simulations. In particular, the equilibrium points of the system and their stability cannot be investigated and obtained from the present model. Therefore, in this study, we implemented the SPVF method, which explicitly exposes the system hierarchy. This method transforms the model into new coordinates using the eigenvectors of the given vector field. In the new coordinates, the model is presented in the form of the SPS, that is, with an explicit hierarchy. The hierarchy of the new mathematical model enables us to decompose the model and divide it into subsystems. In our case, we only divided it into two subsystems: a fast and a slow subsystem. According to the SPVF theory, the fast subsystem can be investigated while the slow system is frozen. We found the equilibrium points the fast system in the new coordinates analytically. We analyzed the stability of the equilibrium points and obtained a single equilibrium point for the new model. To obtain a biological interpretation of the equilibrium point, we inverse transformed the stable equilibrium point into the old coordinates. According to our analytical results and the official reports of the Chinese authorities, we found that the stable equilibrium points obtained from the mathematical model are very close to those of the official reports.