### Model development

The following model, which we call the SEIR-X model, is proposed for modelling the transmission of COVID-19. The model consists of the following mutually exclusive compartments: Susceptible \(\mathrmS(\mathrmt)\), uninfected individuals who are susceptible to the COVID-19 infection; Exposed \(\mathrmE(\mathrmt)\), those who are infected but have not yet entered the active COVID-19 stage; \(\mathrmM\left(\mathrmt\right)\), Mild individuals who are infected, infectious and have mild respiratory illness symptoms such as nasal congestion, runny nose and a sore throat; \(\mathrmC\left(\mathrmt\right)\), Critical individuals who are infected, infectious and have severe symptoms including shortness of breath, chest discomfort and bluish face. After seeking medical advice, critical individuals are classified as either \(\mathrmN_\mathrmH\)—non-hospitalised individuals but still infected, and \(\mathrmH(\mathrmt)\)—hospitalised individuals who are still infected. The last two compartments are the Recovered \(\mathrmR(\mathrmt)\), who were previously infected and were successfully treated, and Death \(\mathrmD(\mathrmt)\). At the same time t, an individual is classified into one and only one compartment. The total population size \(\mathrmN(\mathrmt)\) is assumed to be a constant at time t and well mixed:

$$\mathrmN\left(\mathrmt\right)=\mathrmS\left(\mathrmt\right)+\mathrmE\left(\mathrmt\right)+\mathrmM\left(\mathrmt\right)+\mathrmC\left(\mathrmt\right)+\mathrmN_\mathrmH\left(\mathrmt\right)+\mathrmH\left(\mathrmt\right)+\mathrmR\left(\mathrmt\right)+\mathrmD\left(\mathrmt\right).$$

(1)

To ensure the population size remains constant, we replace all deaths by newborns in the susceptible class. This includes deaths through natural causes, which occur in all states at a constant rate \(\upmu\), and COVID-19 related deaths which occur at a constant rate \(\upomega\). Susceptible individuals may be infected with a circulating strain of COVID-19 at the rate \(\lambda =\upbeta \left(\mathrmM\left(\mathrmt\right)+\mathrmC(\mathrmt)\right)\) and move to the corresponding exposed class \(\mathrmE(\mathrmt)\). Here, \(\upbeta\) is the probability of susceptible individual contracts infection after contact with Mild or Critical individuals with COVID-19. Those with latent infection progress to mild (the M(t) class) due to reactivation of the latent infection at an average period \(\mathrm\alpha \). However, some Mild individuals move to the recovery class \(\mathrmR (\mathrmt)\) at an average period \(\uprho\) due to natural recoveries and the rest of the Mild class individuals move to the critical compartment at an average period \(\upphi\) due to the progression and possibly comorbidities with other diseases, including hypertension, diabetes, cardiovascular disease, and respiratory system disease^{19}. A proportion (type of ratio) of the Critical individuals move to the Non-hospitalised and Hospital compartments at an average period \(\upgamma _1\) and \(\upgamma _2\), respectively. Some Non-hospitalised individuals progress to the recovered compartment \(\mathrmR\left(\mathrmt\right)\) at an average period \(\uptau _1\) through treatment and the rest progress to the death compartment (D) at an average period \(\updelta _1\). Similarly, some of the Hospitalised individuals progress to the recovered compartment \(\mathrmR\left(\mathrmt\right)\) at an average period \(\uptau _2\) through treatment and the rest progress to the death compartment (D) at an average period \(\updelta _2\). A flow diagram of our proposed model is presented in Fig. 3.

In this case, the model can be expressed by the following deterministic system of nonlinear ordinary differential equations:

$$\frac\mathrmdS\mathrmdt=\mathrm\mu N-\upbeta \left(\mathrmM+\mathrmC\right)\mathrmS-\mathrm\mu S$$

(2)

$$\frac\mathrmdE\mathrmdt=\upbeta \left(\mathrmM+\mathrmC\right)\mathrmS-(\mathrm\alpha +\upmu )\mathrmE$$

(3)

$$\frac\mathrmdM\mathrmdt=\mathrm\alpha E-(\upphi +\upmu +\uprho )\mathrmM$$

(4)

$$\frac\mathrmdC\mathrmdt=\mathrm\phi M-(\upomega +\upmu +\upgamma _1+\upgamma _2)\mathrmC$$

(5)

$$\frac\mathrmdN_\mathrmH\mathrmdt=\upgamma _1\mathrmC-\left(\upmu +\updelta _1+\uptau _1\right)\mathrmN_\mathrmH$$

(6)

$$\frac\mathrmdH\mathrmdt=\upgamma _2\mathrmC-(\upmu +\updelta _2+\uptau _2)\mathrmH$$

(7)

$$\frac\mathrmdR\mathrmdt=\mathrm\rho M+\uptau _1\mathrmN_\mathrmH+\uptau _2\mathrmH-\mathrm\mu R$$

(8)

$$\frac\mathrmdD\mathrmdt=\updelta _1\mathrmN_\mathrmH+\updelta _2\mathrmH$$

(9)

Given the non-negative initial conditions for the system above, it is straightforward to show that each of the state variables remains non-negative for all \(\mathrmt>0.\) Moreover, summing Eqs. (2)–(9) we find that the size of the total population, \(\mathrmN(\mathrmt)\) satisfies i.e. \(\mathrmN\left(\mathrmt\right)=\mathrmconstant\).

This shows that the total population size \(\mathrmN(\mathrmt)\) is a constant and it naturally follows that each of the compartment states \(\mathrmS,\mathrmE,\mathrm M,\mathrm C, \mathrmN_\mathrmH,\mathrm H,\mathrm R and D\) is bounded. Given the positivity and boundedness of the system solutions, the feasible region for Eqs. (2)–(9) is given by

$$\mathrmD_1=\left\ \left(\mathrmS,\mathrm E,\mathrm M,\mathrm C, \mathrmN_\mathrmH,\mathrm H,\mathrm R,\mathrm D\right)\in \mathbbR_+^8 :\mathrmS+\mathrmE+\mathrmM+\mathrmC+\mathrmN_\mathrmH+\mathrmH+\mathrmR+\mathrmD=\mathrmN\right\,$$

where \(\mathrmD_1\) is positively invariant.

### Basic reproduction number

The basic reproduction number (\(\mathrmR_0\)) is defined as the expected number of secondary cases created by a single infectious case introduced into a totally susceptible population. The disease can spread in a population only if the basic reproduction number is greater than one. An epidemic occurs when an infection spreads through and infects a significant proportion of a population. A disease-free population is possible when the basic reproduction number is less than one, which means that the disease naturally fades-out^{20,21}.

There are eight states in the modelling system in which five belong to the infected states, i.e. \(\mathrmE,\mathrm M,\mathrm C, \mathrmN_\mathrmH\) and \(\mathrmH\), and three are uninfected states, i.e. \(\mathrmS,\mathrm R\) and \(\mathrmD\). At the infection-free steady-state \(\mathrmE=\mathrmM=\mathrmC=\mathrmN_\mathrmH=\mathrmH=\mathrmR=\mathrmD=0,\) hence \(\mathrmS_0=\mathrmN,\) where \(S_0\) is the initial susceptible population. The Eqs. (3)–(7) are closed, in that they do not involve the derivation of \(\mathrmS\) from steady state value. Also, \(\mathrmR\) and \(\mathrmD\) do not appear in Eqs. (3)–(7)\(,\) and for \((\mathrmE,\mathrmM,\mathrmC,\mathrmN_\mathrmH,\mathrmH)\) we have the following equations:

$$\frac\mathrmdE\mathrmdt=\upbeta \left(\mathrmM+\mathrmC\right)\mathrmS-(\mathrm\alpha +\upmu )\mathrmE$$

(10)

$$\frac\mathrmdM\mathrmdt=\mathrm\alpha E-(\upphi +\upmu +\uprho )\mathrmM$$

(11)

$$\frac\mathrmdC\mathrmdt=\mathrm\phi M-(\upomega +\upmu +\upgamma _1+\upgamma _2)\mathrmC$$

(12)

$$\frac\mathrmdN_\mathrmH\mathrmdt=\upgamma _1\mathrmC-\left(\upmu +\updelta _1+\uptau _1\right)\mathrmN_\mathrmH$$

(13)

$$\frac\mathrmdH\mathrmdt=\upgamma _2\mathrmC-(\upmu +\updelta _2+\uptau _2)\mathrmH$$

(14)

Here, these Ordinary Differential Equations (ODEs) in (10)–(14) are referred to as the infection subsystem, as they only describe the production of newly infected individuals and changes in the states of already infected individuals.

By setting \(\mathbfx=(\mathrmE,\mathrmM,\mathrmC,\mathrmN_\mathrmH,\mathrmH)^\mathrm^\prime\), where the prime denotes transpose, the infection subsystem can be written in the following form:

$$\dot\mathbfx=\left(\mathrmT+\Sigma \right)\mathbfx.$$

(15)

The matrix \(\mathrmT\) corresponds to the transmission, and the matrix \(\Sigma\) to transitions. All epidemiological events that lead to new infections are incorporated in the model via \(\mathrmT\) and other events via \(\Sigma\). If the infected states are indicated with \(\mathrmi\) and \(\mathrmj\) with \(\mathrmi,\mathrmj\in \1, 2, 3, 4, 5\\), then the entry \(\mathrmT_\mathrmij\) is the rate at which individuals in infected state \(\mathrmj\) give rise to individuals in infected state \(\mathrmi\). The matrices \(\mathrmT\) and \(\Sigma\) admit the form

$$\mathrmT=\left(\beginarrayccc0&\upbeta \mathrmS_0& \beginarrayccc\upbeta \mathrmS_0& 0& 0\endarray\\ 0& 0& \beginarrayccc0& 0& 0\endarray\\ \beginarrayc0\\ 0\\ 0\endarray& \beginarrayc0\\ 0\\ 0\endarray& \beginarrayc\beginarrayccc0& 0& 0\endarray\\ \beginarrayccc0& 0& 0\endarray\\ \beginarrayccc0& 0& 0\endarray\endarray\endarray\right)\mathrmand$$

$$\Sigma = \left( \beginarray*20l – \left( \upalpha + \upmu \right) & 0 & 0 & 0 & 0 \\ \alpha & – \left( \upphi + \uprho + \upmu \right) & 0 & 0 & 0 \\ 0 & \upphi & – \left( \upomega + \upgamma _1 + \upgamma _2 + \upmu \right) & 0 & 0 \\ 0 & 0 & \upgamma _1 & – \left( \updelta _1 + \uptau _1 + \upmu \right) & 0 \\ 0 & 0 & \upgamma _2 & 0 & – \left( \updelta _2 + \uptau _2 + \upmu \right) \\ \endarray \right)$$

The next-generation matrix, \(\mathrmK\), is given by^{22} (note the essential minus sign)

$$\textK = – \textT\Sigma ^ – 1 = \textT\left( – \Sigma ^ – 1 \right) = \left( \beginarray*20c \textA_1 & \textA_2 & \textA_3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ \endarray \right)$$

where

$$\mathrmA_1=\frac\mathrmS_0\mathrm\alpha \beta (\mathrm\alpha +\upmu )(\upphi +\uprho +\upmu )+\frac\mathrmS_0\mathrm\alpha \beta \phi (\mathrm\alpha +\upmu )(\upphi +\uprho +\upmu )(\upgamma _1+\upgamma _2+\upomega +\upmu )$$

$$\mathrmA_2=\frac\mathrmS_0\upbeta (\mathrm\alpha +\upmu )(\upphi +\uprho +\upmu )+\frac\mathrmS_0\mathrm\beta \phi (\mathrm\alpha +\upmu )(\upphi +\uprho +\upmu )(\upgamma _1+\upgamma _2+\upomega +\upmu )$$

$$\mathrmA_3=\frac\mathrmS_0\upbeta (\upgamma _1+\upgamma _2+\upomega +\upmu )$$

The dominant eigenvalue of \(\mathrmK\) is the basic reproduction number for COVID-19, which represents the average number of new infections produced by one infected individual. Here, the basic reproduction number can be expressed as:\(\mathrmR_0=\frac\mathrmS_0\mathrm\alpha \beta (\mathrm\alpha +\upmu )(\upphi +\uprho +\upmu )+\frac\mathrmS_0\mathrm\alpha \beta \phi (\mathrm\alpha +\upmu )(\upphi +\uprho +\upmu )(\upgamma _1+\upgamma _2+\upomega +\upmu ).\)

The first term represents the probability of becoming infectious once infected with mean infectious period \(\frac1(\upphi +\uprho +\upmu )\). The second term indicates the expected number of infected individuals due to the comorbidities with other diseases with mean infectious period \(\frac1(\upphi +\uprho +\upmu )(\upgamma _1+\upgamma _2+\upomega +\upmu )\).

We provide detailed analysis, including the existence and stability of the equilibrium points, for the proposed COVID-19 model (2)–(9) in the supplementary materials, see the sections existence of equilibria section and global stability of disease-free equilibrium.

### Parameter estimation and model fitting

We estimated the COVID-19 model parameters from fitting different combinations of parameters in Eqs. (2)–(9) to the actual reported cases in metropolitan and rural health districts in NSW^{23}. In order to parameterize the model, we obtained some of the initial parameter values from literature (see Table 1), and others were estimated from data fitting. The estimation of parameters was carried out using the least squares method which minimises summation of the square errors given by \(\sum \left(\mathrmY\left(\mathrmt,\mathrmq\right)-\mathrmX_\mathrmreal\right)^2\) subject to the COVID-19 model (2)–(9), where \(\mathrmX_\mathrmreal\) is the real reported data, and \(\mathrmY(\mathrmt,\mathrmq)\) denotes the solution of the model corresponding to the number of cases over time \(\mathrmt\) with the set of estimated parameters, denoted by \(\mathrmq\). We assume the initial condition for the state variables in the following way: in metropolitan health district, \(\mathrmN\left(0\right)=\mathrm226,5170,\mathrm E\left(0\right)=\mathrm7,286,\mathrm M\left(0\right)=910,\mathrm C\left(0\right)=300, \mathrmN_\mathrmH\left(0\right)=100,\mathrm H\left(0\right)=40,\mathrm R\left(0\right)=10,\mathrm D\left(0\right)=0,\mathrm S\left(0\right)=\mathrmN\left(0\right)-\mathrmE\left(0\right)-\mathrmM\left(0\right)-\mathrmC\left(0\right)-\mathrmN_\mathrmH\left(0\right)-\mathrmH\left(0\right)-\mathrmR\left(0\right)-\mathrmD\left(0\right)=\mathrm2,256,524\); and in rural health district, \(\mathrmN\left(0\right)=\mathrm103,4370,\mathrm E\left(0\right)=5310,\mathrm M\left(0\right)=620,\mathrm C\left(0\right)=200, \mathrmN_\mathrmH\left(0\right)=80,\mathrm H\left(0\right)=25,\mathrm R\left(0\right)=6,\mathrm D\left(0\right)=0,\mathrm S\left(0\right)=\mathrmN\left(0\right)-\mathrmE\left(0\right)-\mathrmM\left(0\right)-\mathrmC\left(0\right)-\mathrmN_\mathrmH(0)-\mathrmH\left(0\right)-\mathrmR\left(0\right)-\mathrmD\left(0\right)=\mathrm1,028,129\). Figure 4 shows the incidence data of COVID-19 (red dash) and the model fitted curve (blue solid curve).

### Sensitivity of the model to parameters

It is essential to discover how sensitive the COVID-19 model (2)–(9) is to variations in each of its parameters in order to advise intervention strategies that will support in bringing down the infection trajectory. Further, sensitivity analysis will help understand what should be prepared or avoided to mitigate the outbreak of the COVID-19^{28,29}. For this purpose, we calculate the partial rank correlation coefficient (PRCCs)^{30,31} between each of the model parameters and several output variables (Mild and Critical cases) using a Latin Hypercube Sampling. Specially, a uniform distribution is allocated from half to fourfold baseline value (see Table 1) for each model parameters and assigned 100,000 simulations for each. Figures 5 and 6 display the correlations between Mild and Critical cases of metropolitan and rural health areas and the corresponding parameters \(\upbeta ,\mathrm \alpha ,\upphi ,\uprho , \upgamma _1, \upgamma _2\) and \(\upomega\). From Figs. 5, 6 and 7, it is observed that Mild and Critical cases have a strong positive correlation with parameters \(\upbeta\) (transmission rate) and \(\mathrm\alpha \) (progression rate from E to M) in both metropolitan and rural health areas, implying that increasing \(\upbeta\) and \(\mathrm\alpha \) will rise Mild and Critical cases. On the other hand, parameters \(\upphi ,\uprho , \upgamma _1, \upgamma _2\) and \(\upomega\) have a negative correlation with Mild and Critical cases, implying that increasing those parameter values will reduce Mild and Critical cases.

As discussed in earlier sections, the scale and severity of COVID-19 transmission are directly related to the basic reproduction number \(\mathrmR_0.\) Here, we assessed the sensitivity indices of the reproduction number \(\mathrmR_0\). The indices identify how significant each parameter is to \(\mathrmR_0\) and thus the COVID-19 transmission dynamics, and recognize which area should be focused in terms of intervention strategies.

Thus, from the explicit formula for \(\mathrmR_0\), the analytical expression for the sensitivity indices can be derived as a comparative variation in \(\mathrmR_0\) when each parameter changes using the following equation ^{32}:

$$\Upsilon _\mathrmq^\mathrmR_0= \frac\partial \mathrmR_0\partial \mathrmq\times \frac\mathrmq\mathrmR_0,$$

where \(\Upsilon _\mathrmq^\mathrmR_0\) is the sensitivity index of a differentiable \(\mathrmR_0\) for any parameter, \(\mathrmq\).

The sensitivity indices for the model (2)–(9) are graphically presented in Fig. 7. It can be perceived that of all the positives indices, the effective contact rate,\(\upbeta\), is the highest in both metropolitan and rural areas, and therefore the most sensitive parameter. The value of the sensitive index suggests that an increase (or a decrease) in the value of \(\upbeta\) will increases (or decrease) \(\mathrmR_0\) by 100%. However, of all the negative indices displayed in Fig. 7, the recovery rate for the Mild class, denoted by \(\uprho\), is the most sensitive parameter in both metropolitan and rural areas. An increase (or a decrease) of the value of \(\uprho\) will decrease (or increase) \(\mathrmR_0\) by 85%. Moreover, we observed that the progression rate \(\mathrm\alpha \) is comparatively more sensitive in the metropolitan area than in the rural area.

Additionally, contour plots of \(\mathrmR_0\) as a function of other parameters are displayed in Fig. 8 to determine how variations in these parameters affect the basic reproduction number, \(\mathrmR_0\). Figure 8A shows a decrease in \(\mathrmR_0\) with increasing recovery rates for both progression rates from Critical (C) class to Non-hospitalised and Hospital classes. Further, Fig. 8B shows a decreasing trend of the basic reproduction number with both transmission rate, \(\beta\), and recovery rate from M to R, \(\rho\). It is conjectured that \(R_0\) can be brought lower than the threshold of one if efforts are geared towards dropping the contact rate while concurrently improving control and treatment of COVID-19 cases.

### Ethical approval

This study is based on aggregated COVID-19 surveillance data from the New South Wales (NSW) in Australia provided by the NSW government. No confidential information was included because mathematical analyses were performed at the aggregate level. We compiled data from the publicly available website https://data.nsw.gov.au/search/dataset/ds-nsw-ckan-aefcde60-3b0c-4bc0-9af1-6fe652944ec2/details?q = .

link