James A. Njongwe^1,2, Miranda I. Teboh-Ewungkem³, Gideon A. Ngwa^1,4*
1Applied Mathematical and Computer Assisted Modelling Unit, Department of Mathematics, University of Buea, Buea, Cameroon.
²Department of Mathematics and Computer Science, The University of Bamenda, Bamenda, Cameroon.
³US Department of Defense, Government, Forte Meade, MD, USA.
⁴Department of Mathematics and Applied Mathematics, University of Pretoria, Pretoria, South Africa.
DOI: 10.4236/jamp.2026.149181   PDF    HTML   XML   0 Downloads   29 Views   Citations

Abstract

A non-linear ordinary differential equation mathematical model is developed and used to study long-term dynamics for viral disease transmission. The model explicitly incorporates human behaviour and environmental transmissions from viral particles shed into the environment by infected humans. The proposed compartmental model is sufficiently general to capture key epidemiological features common to a broad class of viral infections, while accounting for behavioural changes as individuals transition between disease states and adapt their contact patterns within each stage of infection. The model allows a contaminated environment to act as a reservoir of infection, thereby capturing environmentally mediated transmission arising from viral particles shed by infected humans. The population is stratified into quarantined and non-quarantined classes, allowing the model to assess the impact of non-pharmaceutical interventions on disease progression. Analytical results are used to establish the existence and stability of disease-free and endemic equilibria. The basic reproduction number ${ℛ}_0$ is calculated taking into consideration different combinations of transition and new infection pathways, and its threshold-like nature is well established. Numerical simulations are performed using COVID-19 informed parameter values to illustrate the model dynamics under different scenarios. The theoretical and numerical results from the model’s output highlight the significant role of human behavioural dynamics and environmental transmission in shaping long-term disease outcomes.

Keywords

Behavioural Dynamics, Environmental Contributors, Viral Disease Dynamics

Share and Cite:

1. Introduction and Background

Viral diseases remain a major global public health challenge, causing substantial mortality and persistent social and economic disruption [1]-[3]. Their rapid transmission, continuous evolution, and ability to spread through human, environmental, and behavioural pathways have contributed to some of the most devastating epidemics and pandemics in history, including influenza, Ebola, SARS, HIV/AIDS, and COVID-19 [4]-[8]. Increasing global mobility, urbanization, and pathogen evolution further increase the likelihood of large-scale outbreaks, while limited therapeutic options and viral mutation can complicate their control [9]. Consequently, effective disease control requires not only medical interventions but also surveillance, quarantine, isolation, public health education, and sustained behavioural compliance [10] [11].

Mathematical modelling provides a powerful framework for understanding disease transmission, identifying key epidemiological thresholds, and evaluating intervention strategies. Since the foundational compartmental models of Ross, Hamer, McKendrick, and Kermack [12] [13], numerous models have incorporated quarantine, behavioural responses, environmental transmission, and other epidemiological mechanisms [14]-[18]. However, these features are often considered separately. For example, models incorporating quarantine may neglect reinfection or environmental contamination, while environmentally mediated models may omit behavioural heterogeneity, treatment-centre capacity, or disease transmission from quarantined individuals [16]-[18]. This motivates the development of a more integrated framework capable of capturing the interacting mechanisms that shape long-term viral disease dynamics.

In this study, we develop a general deterministic compartmental framework that unifies key epidemiological, environmental, and behavioural mechanisms within a system of nonlinear ordinary differential equations. The model incorporates non-constant recruitment, limited treatment-centre capacity, reinfection, environmental viral particles, disease-induced mortality and cadaver-mediated transmission, contact tracing, quarantine and isolation, and behavioural heterogeneity among susceptible, suspected, and infectious freely roaming individuals. In particular, non-quarantined populations are stratified according to their adoption of barrier measures, allowing behavioural responses to influence transmission dynamically. To the best of our knowledge, this provides a unified framework for simultaneously capturing these interacting mechanisms while remaining sufficiently flexible for application to a broad range of viral infections through appropriate parameterization. The resulting model provides a mathematical platform for analysing disease thresholds, transmission pathways, intervention strategies, and behavioural feedback in long-term viral disease dynamics.

The rest of the paper is organised as follows. In Section 2, we outline the derivation of the model showing the state variables and parameters, their possible ranges and baseline values used and how they relate together in a conceptual framework; we present a mathematical analysis of the derived model to ascertain that the results are physically realizable and then re-parametrized the model. In Section 3, we present detail analysis of the general disease free model and on the epidemiological model with no barrier measures. Section 4 provides sensitivity analysis and numerical simulations. Finally, discussion of results and conclusion are provided in Section 5.

2. The Mathematical Model

In this section, we give a detailed description of the mathematical model studied. We start with the assumption that there exist $N$ barrier/control measures $b_1,b_2,\cdots ,b_N$ , put in place by a supervisory authority, and that the barrier measures have been combined into sets $C_k$ , $k=0,1,2,\cdots ,N$ , already ordered so that $\varnothing =C_0\subset C_1\subset C_2\subset \cdots \subset C_N$ . Each $C_k$ , is a set that contains up to $k$ barrier measures so that any individual $X$ , indexed by the subscript $k$ , $X_k$ , is considered to adhere to all barrier measures $b_j$ , $j\le k$ . We briefly write $X=\left(X_0,X_1,\cdots ,X_N\right)$ to represent a vector of length $N+1$ containing individuals with different adherent indices. Thus each individual can be adhering to any number of these barrier measures. Our proposed model shall therefore focus on dynamics involving non-pharmaceutical interventions such as quarantining and behavioural changes during a viral outbreak, captured in the form of barrier measures.

2.1. Model Formulation and Description

At any time $t$ , we divide the total human population into compartments representing disease status as well as adherent state to barrier measure(s). The compartments include: $S\left(t\right)$ a vector of susceptible individuals, $S_F\left(t\right)$ , a vector of suspected asymptomatic freely-roaming individuals, $S_Q\left(t\right)$ suspected asymptomatic quarantined individuals, $I_F\left(t\right)$ , a vector of symptomatic non-isolated infected/infectious freely-roaming individuals, $I_H\left(t\right)$ confirmed isolated symptomatic infectious at homes, $I_T\left(t\right)$ confirmed isolated symptomatic infectious at treatment centres, $R_R\left(t\right)$ recovered, $R_D\left(t\right)$ cadavers, and to these, we add a compartment for environmental transmission contamination by $E_v\left(t\right)$ . So, from a base set of $N$ barrier measures, the vectors $S,S_F,I_F$ , are each of length $N+1$ as explained above, and we construct a system comprising $3N+9$ compartments linked by the viral disease pathogen in a dynamical system. The conceptual framework for the model is shown in Figure 1. The following model assumptions may be highlighted:

A1. Non-constant living human population.

A2. Infected individuals can shed viral particles onto the environment.

A3. Individuals can become infected when they come in contact with infected/infectious individuals, the cadavers of the viral disease victims or with a contaminated environment.

A4. Non-Pharmaceutical control is introduced in our model via barrier measures.

A5. Infected, suspected and infectious individuals in the population can be either quarantined or freely roaming.

A6. Quarantined/isolated individuals are assumed to adhere to all the barrier measures.

A7. All freely roaming humans can choose to adhere to none, one or more forms of barrier measures.

A8. On recovery, recovered individuals gradually lose their immunity to the viral disease and return into the susceptible compartment.

A9. The number of available hospital beds is limited by a fixed maximum capacity.

Figure 1. Conceptual framework of the compartmental viral outbreak model. Compartmental flows are described in the text, with notation defined in Tables 1-4.

The state variables, their descriptions, and their corresponding quasi-dimensions are summarised in Table 1.

Table 1. Description of state variables and their quasi-dimensions.

Variable	Description	Quasi-dimension
$H_L$	Density of living humans at time $t$	$\mathscr{H}$
$H$	Density of living humans together with cadavers at time $t$	$\mathscr{H}$
$S$	Vector of dimension $N+1$ representing susceptible humans at time $t$ .	$\mathscr{H}$
$S_F$	Vector of dimension $N+1$ representing suspected humans at time $t$ .	$\mathscr{H}$
$S_Q$	Density of suspected asymptomatic quarantined humans at time $t$	$\mathscr{H}$
$I_F$	Vector of dimension $N+1$ representing infected and infectious individuals at time $t$	$\mathscr{H}$
$I_T$ , $I_H$	Densities of confirmed symptomatic infectious isolated humans at the treatment centres and at homes at any time $t$ respectively	$\mathscr{H}$
$R_R$	Density of recovered humans at time $t$	$\mathscr{H}$
$R_D$	Density of viral disease related deaths at time $t$	$\mathscr{H}$
$R_M$	Density of humans due to natural deaths at time $t$	$\mathscr{H}$
$E_v$	Density of viral particles at time $t$	$\mathcal{V}$

2.2. The Model’s Equations

In this section, we describe the equations that govern the time rate of change of each of the entities in Table 1.

(i) The force of infection: The exposure rate of susceptible humans to infection. We denote by ${\lambda }_{k_N}$ , $k=0,1,2,\cdots ,N$ , the exposure rate of susceptible humans of class $k$ to the infection. We take into consideration susceptible individual’s contacts with infected/infectious humans, cadavers and the contaminated environments. We use a standard incidence model where we append the size of cadavers to $H_L$ , the total living human population to get $H$ . So that:

$\begin{array}{l}{H}_L\left(t\right)=S_Q\left(t\right)+I_T\left(t\right)+I_H\left(t\right)+R_R\left(t\right)+{\displaystyle \sum }_{k=0}^N\left(S_k+S_{Fk}+I_{Fk}\right)\left(t\right);\\ H\left(t\right)=H_L\left(t\right)+R_D\left(t\right).\end{array}$ (1)

For all $k=0,1,2,\cdots ,N$ , the equation for the force of infection takes the following form:

${\lambda }_{k_N}=\frac{1}{H}\left[{\displaystyle \sum }_{j=0}^N\left(\underset{\left(1\right)}{\underbrace{{\tau }_{k,j}{S}_{Fj}}}+\underset{\left(2\right)}{\underbrace{{\phi }_{k,j}{I}_{Fj}}}\right)+\underset{\left(3\right)}{\underbrace{{\epsilon }_k{S}_Q}}+\underset{\left(4\right)}{\underbrace{{\varphi }_k{I}_T}}+\underset{\left(5\right)}{\underbrace{{\zeta }_k{I}_H}}+\underset{\left(6\right)}{\underbrace{d_k{R}_D}}\right]+\underset{\left(7\right)}{\underbrace{{\chi }_{k}Z\left(E_v\right)}}.$ (2)

where $H\left(t\right)>0$ is defined in Equation (1). Each of the terms in Equation (2) are explained below by listing the numbers in the underbraces.

(1) Contacts between susceptible humans with adherence index $k$ and suspected asymptomatic freely-roaming humans with adherence index $j$ . Here, we have simply considered a contact pattern through the formula $S_i+S_{Fj}\stackrel{{\tau }_{i,j}}{\to }{S}_i{S}_{Fj}$ according to standard mass action law with rate constant ${\tau }_{i,j}$ , ${\tau }_{k,j}\ne {\tau }_{j,k}$ . So the mixing pattern $S_k{S}_{Fj}$ will meet with contact rate ${\tau }_{k,j}$ while $S_j{S}_{Fk}$ will meet with contact rate ${\tau }_{j,k}$ . We shall demand for simplicity that the contact weight, ${\overline{\tau }}_{k,j}$ , shall be the same, and write,

${\tau }_{k,j}={\phi }_{S_k}{\overline{\tau }}_{k,j},{\overline{\tau }}_{k,j}=\frac{1}{k+j+1},k,j=0,1,2,\cdots ,N,{\phi }_{S_k}>0.$(3)

Here, ${\phi }_{S_k}>0$ is interpreted for each $k$ as the effective contact rate. Two major requirements for the choice of ${\overline{\tau }}_{k,j}$ are: (i) monotonicity and (ii) positivity. That is, ${\overline{\tau }}_{0,0}>{\overline{\tau }}_{N,N}>0$ and ${\overline{\tau }}_{k,j}>{\overline{\tau }}_{k,j+1}>0$ . One form is shown in Equation (3) which we can identify with the entries of a finite Hilbert matrix [19]. Other forms for ${\overline{\tau }}_{k,j}$ , such as ${\overline{\tau }}_{k,j}={\text{e}}^{-\left(k+j\right)}$ are possible.

(2) Infection from contacts between susceptible humans with adherence index $k$ and freely-roaming infected individual with adherence index $j$ . Using similar argument as in (1) above, we set,

${\phi }_{k,j}={\phi }_{I_k}{\overline{\tau }}_{k,j},k,j=0,1,2,\cdots ,N,{\phi }_{I_k}>0.$(4)

Here, ${\phi }_{I_k}>0$ , $k=0,1,2,\cdots ,N$ is the effective contact rate. From our assumptions, one will expect that ${\phi }_{S_k}\le {\phi }_{I_k}$ .

(3) Infections from contacts between susceptible humans with adherence index $k$ and suspected asymptomatic quarantined individuals.

(4) Infections from contacts between susceptible humans with adherence index $k$ and confirmed isolated symptomatic individuals in the treatment centres.

(5) Infections arising from contacts between susceptible humans with adherence index $k$ and confirmed self isolated symptomatic individuals at homes.

(6) Infections arising from contacts between susceptible humans with adherence index $k$ and the cadavers of the viral infected victims.

(7) Infections from contacts between susceptible humans with adherence index $k$ with the contaminated environments. While some authors, see for example [20], have considered an unbounded linear form for the function $Z\left(E_v\right)$ , we shall assume that $Z:\left[0,\infty \right)\to {ℝ}^{+}$ is a bounded and monotone increasing function. In our

formulation, $Z\left(E_v\right)$ is indexed with two parameters: $B_v$ that measures the supremum of $Z$ and $A_v$ to measure the concentration level where $\text{sup}\left(Z\right)=\frac{1}{2}{B}_v$ .

We demand that, $Z\left(E_v\right)\to 0$ as $E_v\to 0$ and $Z\left(E_v\right)\to B_v$ as $E_v\to \infty$ . One such function is:

$Z\left(E_v\right)=\frac{B_v{E}_v}{A_v+E_v},$ (5)

a Holling’s Type II functional response.

(ii) The susceptible individuals: This class of individuals represented by the variable $S\in {ℝ}^{N+1}$ , comprises humans that have not in anyway been infected by the viral disease. The individuals in the $S$ class are categorized into $N+1$ groups; $S_0,S_1,\cdots ,S_N$ determined by their adherence to barrier measures. For example, $S_k$ , $k=0,1,2,\cdots ,N$ are individuals that adhere up to $k$ forms of barrier measures. The rate of change of the number of susceptible individuals with time decreases when members of this class become exposed and are either moved to the suspected asymptomatic quarantined compartment or to the suspected asymptomatic non-quarantined freely-roaming compartment. While the density increases when some confirmed quarantined individuals who returned negative tests, false suspected non-quarantined freely-roaming individuals who are not sick of the virus and those individuals who recovered from the disease and have lost their immunity all return to the susceptible pool. Figure 2(a) shows the various flows in and out of an arbitrary class $S_k$ , of the susceptible compartment. So, the net rate of change of individuals in the $S_k$ compartment, for $k=0,1,2,\cdots ,N$ at any time $t$ can be written in the form:

(a)

(b)

Figure 2. Conceptual diagram showing the flows for an arbitrary group.

$\begin{array}{l}\frac{\text{d}{S}_k}{\text{d}t}={\left(\text{New}\text{Births}\right)}_k+\left({\displaystyle \sum }_{j=0,j\ne k}^N{a}_{j,k}{S}_j\right)+\left({\displaystyle \sum }_{j=0}^N{\alpha }_{j,k}{S}_{Fj}\right)+{\alpha }_{Qk}{S}_Q\\ +{\alpha }_{Rk}{R}_R-{\lambda }_{k_N}{S}_k-\left(A_k+\mu \right)S_k,\text{where}{A}_k=\left({\displaystyle \sum }_{j=0,j\ne k}^N{a}_{k,j}\right),\end{array}$ (6)

where the parameters in ${\lambda }_{k_N}$ are as given in Equation (2). We suppose that the total out flow rate for individuals from the $S_k$ class due to exposure to the viral infection is ${\lambda }_{k_N}$ , then we expect the conservation equation ${\displaystyle {\sum }_{j=0}^N{\lambda }_{k,j}{S}_k}+{\lambda }_{kQ}{S}_k={\lambda }_{k_N}{S}_k$ , where ${\lambda }_{k_N}$ is as defined in Equation (2) to hold. Now, let ${\theta }_k$ be the fraction of ${\lambda }_{k_N}{S}_k$ individuals that were identified and placed in quarantine and the remaining ${\theta }_{k,j}$ be the fraction that escaped detection and so enter the suspected freely roaming class. Our conservation equation will now take the form; ${\theta }_k+{\displaystyle {\sum }_{j=0}^N{\theta }_{k,j}}=1$ ; for ${\lambda }_{kQ}={\theta }_k{\lambda }_{k_N}$ , and ${\lambda }_{k,j}={\theta }_{k,j}{\lambda }_{k_N}$ .

(iii) The suspected asymptomatic freely-roaming individuals: This class comprises individuals who came in contact with or where near a known viral infected patient. They are described by the variable $S_F\in {ℝ}^{N+1}$ and are categorized into $N+1$ freely-roaming suspected groups represented by the variables $S_{F0},S_{F1},\cdots ,S_{FN}$ . Outflow from this class is driven by the development of symptoms or that they were not infected of the viral disease at all. The rate of change of the density of members in this group decreases when some of these individuals return to the susceptible pool because they were not sick of the viral disease, some developed symptoms and progressed to the infected/infectious freely roaming compartment and others who developed symptoms, were identified and isolated either in the treatment centres or at homes. The density of these individuals increases when the susceptible individuals are exposed by having come in contact with or being associated with any of the possible infectious cases. For a particular class, $S_{Fk}$ , $k=0,1,2,\cdots ,N$ of the $S_F$ compartment, the flows in and out of this subclass is as shown on Figure 2(b) and the equation that governs the rate of change of the $S_{Fk}$ individuals at any time $t$ , is given by:

$\frac{\text{d}{S}_{Fk}}{\text{d}t}=\left({\displaystyle \sum }_{j=0}^N{\lambda }_{j,k}{S}_j\right)+\left({\displaystyle \sum }_{j=0,j\ne k}^N{b}_{j,k}{S}_{Fj}\right)-\left(f_{k_N}+G_k+\mu \right)S_{Fk}$ (7)

where $f_{k_N}={\xi }_{Fk}+{\omega }_{Fk}+{\displaystyle \sum }_{j=0}^N\left({\alpha }_{k,j}+{\beta }_{k,j}\right)$ and $G_k=\left({\displaystyle \sum }_{j=0,j\ne k}^N{b}_{k,j}\right)$ , $k=0,1,2,\cdots ,N$ .

(iv) The suspected asymptomatic quarantined individuals: This class comprises the other group of suspected individuals who were in contact with or in the vicinity of a known infected patient but are traced, via contact tracing or some other means and quarantined. Members in this class, described by the variable $S_Q$ , are either self-quarantined or isolated in approved facilities. Individuals in this class move to the susceptible class (following a negative test result) or they move to the confirmed infected/infectious symptomatic isolated classes (following a positive test result) and depending on the availability of space at the treatment centres. The rate of change of the density of individuals increases when a fraction of the susceptible individuals who got exposed get traced and quarantined. The density is reduced when some of these individuals returned to the susceptible pool after returning negative viral test results, while those who are confirmed positive, either move to the treatment centres depending on whether there is space or are forced to isolate at homes. Thus:

$\frac{\text{d}{S}_Q}{\text{d}t}=\left({\displaystyle \sum }_{j=0}^N{\lambda }_{jQ}{S}_j\right)-\left(f_{Q_N}+\mu \right)S_Q;f_{Q_N}={\displaystyle \sum }_{j=0}^N{\alpha }_{Qj}+\gamma +\beta .$ (8)

(v) The symptomatic non-isolated infected/infectious freely-roaming population: This class is made up of those asymptomatic non-quarantined freely-roaming individuals who are infectious but do not exhibit symptoms or may exhibit symptoms but for some reason, cannot get a test. These are the super spreaders. They are described by the variable $I_F\in {ℝ}^{N+1}$ and are categorized into $N+1$ subgroups represented by the indexed variables $I_{F0},I_{F1},\cdots ,I_{FN}$ . We believe these are key disease drivers in the community and their level of contribution to disease spread via human-to-human contacts and rate of viral shedding into the environment would be determined by the respective level of compliance to the barrier measures. We assume that the rate of viral shedding in to the environment, in terms of the subgroups is in the order $I_{F0}\ge I_{F1}\ge \cdots \ge I_{FN}$ whereby individuals in the subgroup $I_{F0}$ could be termed super disease spreaders. The rate of change of the density of symptomatic infected/infectious freely roaming individuals in the population increases when a fraction of the suspected asymptomatic freely-roaming individuals developed symptoms, progressed to this class. Outflow from this class occurs when some of these individuals are removed by recovery from the disease at rate $r_{Fk}$ , by death caused by the virus at rate ${\delta }_{Fk}$ or to isolated units at homes or in treatment centres at rates ${\rho }_{Fk}$ and ${\sigma }_{Fk}$ respectively while others die due to natural death. For an arbitrary class, $I_{Fk}$ , of the $I_F$ compartment, the equation that governs the rate of change of the $I_{Fk}$ individuals at any time $t$ , is given by:

$\frac{\text{d}{I}_{Fk}}{\text{d}t}=\left({\displaystyle \sum }_{j=0}^N{\beta }_{j,k}{S}_{Fj}\right)+\left({\displaystyle \sum }_{j=0,j\ne k}^N{c}_{j,k}{I}_{Fj}\right)-\left(g_k+{\mathcal{C}}_k+\mu \right)I_{Fk},k=0,1,\cdots ,N;$ (9)

where ${\mathcal{C}}_k=\left({\displaystyle \sum }_{j=0,j\ne k}^N{c}_{k,j}\right)$ and $g_k={\sigma }_{Fk}+{\rho }_{Fk}+r_{Fk}+{\delta }_{Fk}$ .

(vi) The confirmed isolated symptomatic infectious individuals at homes: This class of individuals comprises all those individuals who are now presenting with symptoms. They are described by the variable $I_H$ . Recruitment into this class depends on whether the treatment centre has attained it carrying capacity, $B_{\text{max}}$ or not. The rate of change of the number of $I_H$ increases as a result of flow from the $S_Q$ and were unable to be admitted at the treatment centres because there is no space in the treatment centres or they chose to isolate themselves at homes, from $S_F$ class who developed symptoms, and then forced to isolate themselves at homes. Further more, increment in this class is also caused by cases from the $I_F$ class. The rate of change of the density of $I_H$ class is reduced when some of the individuals are removed, or moved to the treatment centres. Removal could be as a result of recovery at rate $r_H$ or as a result of disease induced death at rate ${\delta }_H$ , or of natural cause at rate $\mu$ . At time $t$ , we have:

$\frac{\text{d}{I}_H}{\text{d}t}=F_{H_N}\left(S_Q,S_F,I_F\right)-\left({\beta }_H+\mu \right)I_H,$ (10)

where the function $F_{H_N}$ measures the flow rate into the $I_H$ compartment as a function of the size of the treatment centre and it is given by:

$F_{H_N}\left(S_Q,S_F,I_F\right)=\{\begin{array}{ll}\left(\beta +\gamma \right)S_Q+{\displaystyle \sum }_{j=0}^N{\mathcal{K}}_j,\hfill & \text{if}〚I_T〛>B_{\text{max}}\hfill \\ \beta S_Q+{\displaystyle \sum }_{j=0}^N\left({\xi }_{Fj}{S}_{Fj}+{\rho }_{Fj}{I}_{Fj}\right),\hfill & \text{if}〚I_T〛\le B_{\text{max}}\hfill \end{array},$ (11)

for ${\mathcal{K}}_j={\omega }_{Fj}{S}_{Fj}+{\xi }_{Fj}{S}_{Fj}+{\sigma }_{Fj}{I}_{Fj}+{\rho }_{Fj}{I}_{Fj}$ , $j=0,1,2,\cdots ,N$ , and ${\beta }_H=r_H+{\eta }_H+{\delta }_H$ is the total out flow rate from this compartment with:

${\eta }_H=\{\begin{array}{ll}0,\hfill & \text{if}〚I_T〛>B_{\text{max}}\hfill \\ \eta ,\hfill & \text{if}〚I_T〛\le B_{\text{max}},\hfill \end{array}$ (12)

where $〚I_T〛$ is is the floor of $I_T$ , which is also the greatest integer smaller or equal to $I_T$ . Mathematically, $〚I_T〛=I\iff n\le I<n+1$ , $n\in \mathbb{N}\cup \left\{0\right\}$ . Thus $〚I_T〛$ is the effective number of individuals at the treatment centres, capped by the maximum treatment centre carrying capacity $B_{\text{max}}$ .

(vii) The confirmed isolated symptomatic infectious individuals in treatment centres: This is the class of individuals who at some point had tested positive and were asymptomatic, not exhibiting any viral disease-like symptoms but are now presenting with symptoms. They are described by the variable $I_T$ , flow into this compartment is greatly conditioned on the size of the treatment centre. The rate of change of the $I_T$ individuals increases as a result of flow from the $S_Q$ class by individuals who returned positive viral tests, from $S_F$ class who developed symptoms, were traced and then taken to the treatment centres. Increment is also caused by cases from the $I_H$ class, when it is noticed that there is space at the treatment centre and, from the $S_F$ class. The rate of change of the density of $I_T$ individuals is reduced when some of the individuals are removed as a result of recovery at rate $r_T$ or as a result of disease related death at rate ${\delta }_T$ and lastly, as a result of natural death at rate $\mu$ . At time $t$ , we have:

$\frac{\text{d}{I}_T}{\text{d}t}=F_{T_N}\left(S_Q,I_H,S_F,I_F\right)-\left({\beta }_T+\mu \right)I_T;{\beta }_T=r_T+{\delta }_T,$ (13)

where the function $F_{T_N}$ measures the flow rate into the $I_T$ compartment as a function of the size of the treatment centre and it is modelled by:

$F_{T_N}\left(S_Q,I_H,S_F,I_F\right)=\{\begin{array}{ll}0,\hfill & \text{if}〚I_T〛>B_{\text{max}}\hfill \\ \gamma S_Q+\eta I_H+{\displaystyle \sum }_{j=0}^N\left({\omega }_{Fj}{S}_{Fj}+{\sigma }_{Fj}{I}_{Fj}\right),\hfill & \text{if}〚I_T〛\le B_{\text{max}}\hfill \end{array}.$ (14)

(viii) The recovered individuals: The class of recovered individuals denoted by $R_R$ , contains all individuals who once got infected by the virus and did not die of the disease. Since recovery from the disease may not confer permanent immunity against further infections, these recovered individuals eventually return into the susceptible pool and can be reinfected if proper precautions are not taken. The equation governing the rate of change of recovered individuals at time $t$ , takes the form:

$\frac{\text{d}{R}_R}{\text{d}t}=\left({\displaystyle \sum }_{j=0}^N{r}_{Fj}{I}_{Fj}\right)+r_T{I}_T+r_H{I}_H-\left({\alpha }_{R_N}+\mu \right)R_R;$ (15)

where ${\alpha }_{R_N}={\displaystyle {\sum }_{j=0}^N{\alpha }_{Rj}}$ .

(ix) The cadavers: Here, we have two kinds of deaths; (i) disease induced deaths, collected into the compartment $R_D$ at rates ${\delta }_{Fk}$ , ${\delta }_T$ and ${\delta }_H$ form the $I_F$ , $I_T$ and $I_H$ compartments respectively. (ii) natural deaths collected into the compartment $R_M$ at a constant rate $\mu$ from all living population classes. It is known that the dead bodies of some viral disease victims are still very infectious and can still infect susceptible individuals upon effective contact. Thus, in this work, we consider $r_D$ to be the rate at which cadavers are properly disposed of through burial or cremation. The rate of change of the population of disease induced deceased individuals at time $t$ is:

$\frac{\text{d}{R}_D}{\text{d}t}=\left({\displaystyle \sum }_{j=0}^N{\delta }_{Fj}{I}_{Fj}\right)+{\delta }_T{I}_T+{\delta }_H{I}_H-r_D{R}_D.$ (16)

In Equation (16), $r_D{R}_D$ are disease related disposed cadavers. It is a collection class modelled by; $\frac{\text{d}{D}_D}{\text{d}t}=r_D{R}_D$ . Next, we define $R_M$ , another collection class

that tracks all natural deaths, occurring at rate $\mu$ , from all the living population classes. It is modelled by:

$\frac{\text{d}{R}_M}{\text{d}t}=\mu S_Q+\mu I_H+\mu I_T+\mu R_R+{\displaystyle \sum }_{j=0}^N\left(\mu S_j+\mu S_{Fj}+\mu I_{Fj}\right)=\mu H_L.$ (17)

We use the classes $R_M$ and $D_D$ , only as place holders.

(x) The environmental transmission contributors: This compartment comprises all viral particles shed by infected humans into the environment and it is described by the variable $E_v$ . The shedding of these viral particles by infected/infectious humans occur at different rates depending on the adherence index of the infected human. Thus, the rate of change of the density of viral particles shedded into the environment at time $t$ , is modelled by the equation:

$\frac{\text{d}{E}_v}{\text{d}t}={\displaystyle \sum }_{j=0}^N\left({\psi }_{Fj}{S}_{Fj}+{\nu }_{Fj}{I}_{Fj}\right)+{\nu }_Q{S}_Q+{\nu }_H{I}_H+{\nu }_T{I}_T-d{E}_v.$ (18)

In the absence of viral disease, $E_v$ decays steadily to zero at rate $d$ .

(xi) Dynamics of the total human populations: The total living human population is not constant since the model accounts for births and deaths. We assume that all new births, enter the class of susceptible individuals who adhere to all $N$

control measures. Let $b\left(H_L\right)=r_h\left(1-\frac{H_L}{K_h}\right)$ , where $r_h$ and $K_h$ are positive

constants. So that the net rate of new births that enters the susceptible compartment, $S$ , at any time $t$ is given by:

${\left(\text{New Births}\right)}_k=\{\begin{array}{ll}b\left(H_L\right)H_L,\hfill & k=N\hfill \\ 0,\hfill & k\ne N\hfill \end{array}.$ (19)

At any given time $t$ , the rate of change of the total living population is:

$\frac{\text{d}{H}_L}{\text{d}t}=\left(b\left(H_L\right)-\mu \right)H_L-{\delta }_T{I}_T-{\delta }_H{I}_H-\left({\displaystyle \sum }_{j=0}^N{\delta }_{Fj}{I}_{Fj}\right).$ (20)

Equation (20) shows the dependence of the size of the living population on disease related deaths. If ${\delta }_T={\delta }_H={\delta }_{Fj}=0$ for $j=0,1,2,\cdots ,N$ , then all deaths in the system are only due to natural causes and the total living human population

satisfies the equation: $\frac{\text{d}{H}_L}{\text{d}t}=r_h\left(1-\frac{H_L}{K_h}\right)H_L-\mu H_L$ ; which can be solved directly to have $H_L\left(t\right)=\frac{K}{1+B{\text{e}}^{-rt}}$ , where $K=\frac{K_h\left(r_h-\mu \right)}{r_h}$ , $B=\frac{K_h}{r_h}{\text{e}}^{-\left(r_h-\mu \right)c}$ , $r=r_h-\mu$

and $c\in ℝ$ . Clearly, if (i) $r_h>\mu$ , $H_L\left(t\right)\to K$ as $t\to \infty$ and if (ii) $r_h\le \mu$ , $H_L\left(t\right)\to 0$ as $t\to \infty$ . In either case, $H_L\left(t\right)$ remains bounded as $t\to \infty$ . We shall assume that $r_h>\mu$ . Also, applying similar arguments, the equation for the rate of change of the total living human population and the cadavers is given as:

$\frac{\text{d}H}{\text{d}t}=r_h\left(1-\frac{H_L}{K_h}\right)H_L-\mu H_L-r_D{R}_D.$ (21)

Putting all the equations together, we get:

$\begin{array}{l}\frac{\text{d}{S}_k}{\text{d}t}={\left(\text{New Births}\right)}_k+\left({\displaystyle \sum }_{j=0,j\ne k}^N{a}_{j,k}{S}_j\right)+\left({\displaystyle \sum }_{j=0}^N{\alpha }_{j,k}{S}_{Fj}\right)+{\alpha }_{Qk}{S}_Q\\ +{\alpha }_{Rk}{R}_R-{\lambda }_{k_N}{S}_k-\left(A_k+\mu \right)S_k,k=0,1,2,\cdots ,N;\\ \frac{\text{d}{S}_{Fk}}{\text{d}t}=\left({\displaystyle \sum }_{j=0}^N{\lambda }_{j,k}{S}_j\right)+\left({\displaystyle \sum }_{j=0,j\ne k}^N{b}_{j,k}{S}_{Fj}\right)-\left(f_{k_N}+G_k+\mu \right)S_{Fk},k=0,1,2,\cdots ,N;\\ \frac{\text{d}{S}_Q}{\text{d}t}=\left({\displaystyle \sum }_{j=0}^N{\lambda }_{jQ}{S}_j\right)-\left(f_{Q_N}+\mu \right)S_Q;\\ \frac{\text{d}{I}_{Fk}}{\text{d}t}=\left({\displaystyle \sum }_{j=0}^N{\beta }_{j,k}{S}_{Fj}\right)+\left({\displaystyle \sum }_{j=0,j\ne k}^N{c}_{j,k}{I}_{Fj}\right)-\left(g_k+{\mathcal{C}}_k+\mu \right)I_{Fk},k=0,1,2,\cdots ,N;\\ \frac{\text{d}{I}_H}{\text{d}t}=F_{H_N}\left(S_Q,S_F,I_F\right)-\left({\beta }_H+\mu \right)I_H;\\ \frac{\text{d}{I}_T}{\text{d}t}=F_{T_N}\left(S_Q,I_H,S_F,I_F\right)-\left({\beta }_T+\mu \right)I_T;\\ \frac{\text{d}{R}_R}{\text{d}t}=\left({\displaystyle \sum }_{j=0}^N{r}_{Fj}{I}_{Fj}\right)+r_T{I}_T+r_H{I}_H-\left({\alpha }_{R_N}+\mu \right)R_R;\\ \frac{\text{d}{R}_D}{\text{d}t}=\left({\displaystyle \sum }_{j=0}^N{\delta }_{Fj}{I}_{Fj}\right)+{\delta }_T{I}_T+{\delta }_H{I}_H-r_D{R}_D;\\ \frac{\text{d}{E}_v}{\text{d}t}={\displaystyle \sum }_{j=0}^N\left({\psi }_{Fj}{S}_{Fj}+{\nu }_{Fj}{I}_{Fj}\right)+{\nu }_Q{S}_Q+{\nu }_H{I}_H+{\nu }_T{I}_T-d{E}_v.\end{array}\}$ (22)

The equations for the total living human populations are given by:

$\begin{array}{l}\frac{\text{d}{H}_L}{\text{d}t}=r_h\left(1-\frac{H_L}{K_h}\right)H_L-\mu H_L-{\delta }_T{I}_T-{\delta }_H{I}_H-\left({\displaystyle \sum }_{j=0}^N{\delta }_{Fj}{I}_{Fj}\right);\\ \frac{\text{d}H}{\text{d}t}=r_h\left(1-\frac{H_L}{K_h}\right)H_L-\mu H-\left(r_D-\mu \right)R_D.\end{array}\}$ (23)

where ${\lambda }_{k_N}$ is as given in Equation (2); $A_k$ , Equation (6); $G_k$ and $f_{k_N}$ , Equation (7); $f_{Q_N}$ , Equation (8); ${\mathcal{C}}_k$ and $g_k$ , Equation (9); $F_{H_N}\left(S_Q,S_F,I_F\right)$ , Equation (11) and $F_{T_N}\left(S_Q,I_H,S_F,I_F\right)$ Equation (14).

One form of initial conditions could be those that start off the process with some initial density of the form:

$S_k\left(0\right)=S_k^0,S_{Fk}\left(0\right)=S_{Fk}^0,I_{Fk}\left(0\right)=I_{Fk}^0,S_Q\left(0\right)=S_Q^0,I_H\left(0\right)=I_H^0,$ (24)

$I_T\left(0\right)=I_T^0,R_R\left(0\right)=R_R^0,R_D\left(0\right)=R_D^0,E_v\left(0\right)=E_v^0,$ (25)

where the variables with superscript $0$ are typical variables at time $t=0$ whose values will be provided as initial conditions.

Definition 2.1 (Realistic solution). A solution of system (22) is called realistic if it is non-negative and bounded.

For all intents and purposes, since the rates are polynomials and rational functions, we can easily establish using standard techniques as shown in [21] that system (22) is mathematically and physically well-posed. Also, the solutions lie in the bounded region,

$\begin{array}{l}\Omega =\{\left(S_i\left(t\right),S_{Fi}\left(t\right),S_Q\left(t\right),I_{Fi}\left(t\right),I_T\left(t\right),I_H\left(t\right),R_R\left(t\right),R_D\left(t\right),E_v\left(t\right)\right)\stackrel{}{}\\ \in {ℝ}_{+}^{3N+9};H_L\left(t\right)\le \frac{r_h{K}_h}{4\mu }\text{and}{E}_v\left(t\right)\le \frac{Q_{\text{max}}}{d}\}\subset {ℝ}_{+}^{3N+9},\\ Q_{\text{max}}=\left[{\displaystyle \sum }_{j=0}^N\left({\psi }_{Fj}+{\nu }_{Fj}\right)+{\nu }_Q+{\nu }_H+{\nu }_T\right]\left(\frac{r_h{K}_h}{4\mu }\right),i\in ℐ.\end{array}\}$(26)

2.3. Non-Dimensionalisation and Re-Parametrisation

To scale the system, we make the following change of variables:

$\begin{array}{l}{\tilde{S}}_k=\frac{S_k}{S_k^0},{\tilde{S}}_{Fk}=\frac{S_{Fk}}{S_{Fk}^0},{\tilde{S}}_Q=\frac{S_Q}{S_Q^0},{\tilde{I}}_{Fk}=\frac{I_{Fk}}{I_{Fk}^0},{\tilde{I}}_T=\frac{I_T}{I_T^0},{\tilde{I}}_H=\frac{I_H}{I_H^0},{\tilde{R}}_R=\frac{R_R}{R_R^0},\\ {\tilde{R}}_D=\frac{R_D}{R_D^0},{\tilde{R}}_M=\frac{R_M}{R_M^0},{\tilde{D}}_D=\frac{D_D}{D_D^0},\tilde{H}=\frac{H}{H^0},{\tilde{H}}_L=\frac{H_L}{H_L^0},{\tilde{E}}_v=\frac{E_v}{E_v^0},\tau =\frac{t}{T^0},\end{array}$

where $S_k^0$ , $S_{Fk}^0$ , $S_Q^0$ , $I_{Fk}^0$ , $I_T^0$ , $I_H^0$ , $R_R^0$ , $R_D^0$ , $R_M^0$ , $D_D^0$ , $E_v^0$ , $H_L^0$ and $H^0$ , are reference quantities associated with the different human and viral compartments, and $T^0$ is a characteristic time frame for the system. We set;

$\begin{array}{l}K=\frac{K_h\left(r_h-\mu \right)}{r_h};r=r_h-\mu ;\forall k\in ℐ,S_k^0=S_{Fk}^0=S_Q^0=I_{Fk}^0=I_H^0=K;\\ I_T^0=R_R^0=R_D^0=R_M^0=D_D^0=H_L^0=H^0=K;E_v^0=A_v;T^0=\frac{1}{\mu };\end{array}$ (27)

and then define the following dimensionless parameter groupings:

(28)

With $\ast 1\in \left\{Fk,T,H\right\}$ . This leads to the scaled system:

$\begin{array}{l}\frac{\text{d}{\tilde{S}}_k}{\text{d}\tau }={\gamma }_{S_k}[r_k\left(1-h{\tilde{H}}_L\right){\tilde{H}}_L+{\displaystyle \sum }_{j=0,j\ne k}^N{\tilde{a}}_{j,k}{\tilde{S}}_j+{\displaystyle \sum }_{j=0}^N{\tilde{\alpha }}_{j,k}{\tilde{S}}_{Fj}\\ \stackrel{}{}+{\tilde{\alpha }}_{Qk}{\tilde{S}}_Q+{\tilde{\alpha }}_{Rk}{\tilde{R}}_R-a_k{\lambda }_k\left(\tilde{y}\right){\tilde{S}}_k-{\tilde{S}}_k],\forall k=0,1,2,\cdots ,N;\\ \frac{\text{d}{\tilde{S}}_{Fk}}{\text{d}\tau }={\gamma }_{S_{Fk}}\left[{\displaystyle \sum }_{j=0}^N{\tilde{\theta }}_{j,k}{\lambda }_j\left(\tilde{y}\right){\tilde{S}}_j+{\displaystyle \sum }_{j=0,j\ne k}^N{\tilde{b}}_{j,k}{\tilde{S}}_{Fj}-{\tilde{S}}_{Fk}\right],\forall k=0,1,2,\cdots ,N;\\ \frac{\text{d}{\tilde{S}}_Q}{\text{d}\tau }={\gamma }_{S_Q}\left[{\displaystyle \sum }_{j=0}^N{\widehat{\theta }}_j{\lambda }_j\left(\tilde{y}\right){\tilde{S}}_j-{\tilde{S}}_Q\right];\\ \frac{\text{d}{\tilde{I}}_{Fk}}{\text{d}\tau }={\gamma }_{I_{Fk}}\left[{\displaystyle \sum }_{j=0}^N{\tilde{\beta }}_{j,k}{\tilde{S}}_{Fj}+{\displaystyle \sum }_{j=0,j\ne k}^N{\tilde{c}}_{j,k}{\tilde{I}}_{Fj}-{\tilde{I}}_{Fk}\right];\\ \frac{\text{d}{\tilde{I}}_H}{\text{d}\tau }={\gamma }_{I_H}\left[F_h\left({\tilde{S}}_Q,{\tilde{S}}_F,{\tilde{I}}_F\right)-{\tilde{I}}_H\right];\\ \frac{\text{d}{\tilde{I}}_T}{\text{d}\tau }={\gamma }_{I_T}\left[F_t\left({\tilde{S}}_Q,{\tilde{I}}_H,{\tilde{S}}_F,{\tilde{I}}_F\right)-{\tilde{I}}_T\right];\\ \frac{\text{d}{\tilde{R}}_R}{\text{d}\tau }={\gamma }_{R_R}\left[{\displaystyle \sum }_{j=0}^N{\tilde{r}}_{Fj}{\tilde{I}}_{Fj}+{\tilde{r}}_T{\tilde{I}}_T+{\tilde{r}}_H{\tilde{I}}_H-{\tilde{R}}_R\right];\\ \frac{\text{d}{\tilde{R}}_D}{\text{d}\tau }={\gamma }_{R_D}\left[{\displaystyle \sum }_{j=0}^N{\widehat{\delta }}_{Fj}{\tilde{I}}_{Fj}+{\widehat{\delta }}_T{\tilde{I}}_T+{\widehat{\delta }}_H{\tilde{I}}_H-{\tilde{R}}_D\right];\\ \frac{\text{d}{\tilde{E}}_v}{\text{d}\tau }={\gamma }_{E_v}\left[{\displaystyle \sum }_{j=0}^N\left({\tilde{\psi }}_{Fj}{\tilde{S}}_{Fj}+{\tilde{\nu }}_{Fj}{\tilde{I}}_{Fj}\right)+{\tilde{\nu }}_Q{\tilde{S}}_Q+{\tilde{\nu }}_H{\tilde{I}}_H+{\tilde{\nu }}_T{\tilde{I}}_T-{\tilde{E}}_v\right];\end{array}$(29)

and the total populations satisfy the scaled equations:

$\begin{array}{l}\frac{\text{d}{\tilde{H}}_L}{\text{d}\tau }=\widehat{r}\left(1-{\tilde{H}}_L\right){\tilde{H}}_L-{\tilde{\delta }}_T{\tilde{I}}_T-{\tilde{\delta }}_H{\tilde{I}}_H-{\displaystyle \sum }_{j=0}^N{\tilde{\delta }}_{Fj}{\tilde{I}}_{Fj};\\ \frac{\text{d}\tilde{H}}{\text{d}\tau }=\tilde{r}\left(1-h{\tilde{H}}_L\right){\tilde{H}}_L-\tilde{H}-{\tilde{r}}_D{\tilde{R}}_D.\end{array}$(30)

where:

$\begin{array}{c}{\lambda }_k\left(\tilde{y}\right)={\displaystyle \sum }_{j=0}^N\left({\widehat{\tau }}_{k,j}\frac{{\tilde{S}}_{Fj}}{\tilde{H}}+{\widehat{\phi }}_{k,j}\frac{{\tilde{I}}_{Fj}}{\tilde{H}}\right)+{\widehat{\epsilon }}_k\frac{{\tilde{S}}_Q}{\tilde{H}}+{\widehat{\varphi }}_k\frac{{\tilde{I}}_T}{\tilde{H}}+{\widehat{\zeta }}_k\frac{{\tilde{I}}_H}{\tilde{H}}\\ +{\widehat{d}}_k\frac{{\tilde{R}}_D}{\tilde{H}}+{\widehat{\chi }}_k\left(\frac{{\tilde{E}}_v}{1+{\tilde{E}}_v}\right);\end{array}$(31)

$\begin{array}{l}{F}_h\left({\tilde{S}}_Q,{\tilde{S}}_F,{\tilde{I}}_F\right)=\{\begin{array}{ll}\tilde{\gamma }{\tilde{S}}_Q+{\displaystyle \sum }_{j=0}^N\left({\tilde{\xi }}_{Fj}{\tilde{S}}_{Fj}+{\tilde{\rho }}_{Fj}{\tilde{I}}_{Fj}\right),\hfill & \text{if}〚{\tilde{I}}_T〛>b_{\text{max}}\hfill \\ \tilde{\beta }{\tilde{S}}_Q+{\displaystyle \sum }_{j=0}^N\left({\tilde{\omega }}_{Fj}{\tilde{S}}_{Fj}+{\tilde{\sigma }}_{Fj}{\tilde{I}}_{Fj}\right),\hfill & \text{if}〚{\tilde{I}}_T〛\le b_{\text{max}}\hfill \end{array};\hfill \end{array}$(32)

and

$\begin{array}{l}{F}_t\left({\tilde{S}}_Q,{\tilde{I}}_H,{\tilde{S}}_F,{\tilde{I}}_F\right)=\{\begin{array}{ll}0,\hfill & \text{if}〚{\tilde{I}}_T〛>b_{\text{max}}\hfill \\ \widehat{\gamma }{\tilde{S}}_Q+\tilde{\eta }{\tilde{I}}_H+{\displaystyle \sum }_{j=0}^N{\widehat{\mathcal{Q}}}_j,\hfill & \text{if}〚{\tilde{I}}_T〛\le b_{\text{max}}\hfill \end{array}.\hfill \end{array}$(33)

For, ${\widehat{\mathcal{Q}}}_j=\left({\widehat{\omega }}_{Fj}{\tilde{S}}_{Fj}+{\widehat{\sigma }}_{Fj}{\tilde{I}}_{Fj}\right)$ .

Table 2. Model’s parameter values and their Quasi-dimensions.

Parameter	Description	Units	Range of possible values	Baseline values	References
$r_h$	Recruitment rate of humans	$T^{-1}$	$\left[\frac{6}{1000\times 365},\frac{46.6}{1000\times 365}\right]$	$9.58904\times {10}^{-5}$	[22]
$\mu$	Natural mortality rate of humans.	$T^{-1}$	$\left[\frac{1}{89.8\times 365},\frac{1}{54.4\times 365}\right]$	$4.56621\times {10}^{-5}$	[22]
$K_h$	Human environmental carrying capacity	$\mathscr{H}$	Varies	106	Assumed
${\theta }_{k,j}\left({\theta }_k\right)$	Proportions of susceptible individuals with adherence index $k$ who become suspected asymptomatic freely-roaming individuals with adherence index $j$ , (respectively, suspected asymptomatic quarantined individuals).	1	$\left[0,1\right]$	0.75(0.25)	Assumed
$f_{k_N}$ ($f_{Q_N}$ )	Total outflow rate for suspected asymptomatic freely-roaming (and suspected quarantined) individuals with adherence index $k$ .	$T^{-1}$	$\left[\frac{1}{14},\frac{1}{2}\right]$	$\frac{1}{8}$	[23]
${\alpha }_{k,j}$	Flow rate from the $S_{Fk}$ class in the suspected freely-roaming compartment to the $S_j$ class in the susceptible compartment.	$T^{-1}$	$\left[\frac{1}{500},\frac{1}{2}\right]$	0.04375	Assumed, weighted with $f_{k_N}$ .
${\beta }_{k,j}$	Rate at which individuals from the $S_{Fk}$ class in the suspected freely-roaming compartment progressed to the $I_{Fj}$ class in the symptomatic infectious freely-roaming compartment.	$T^{-1}$	$\left[\frac{1}{500},\frac{1}{2}\right]$	0.025	Assumed, weighted with $f_{k_N}$ .
${\xi }_{Fk}$	Rate at which suspected freely roaming individuals are removed and isolated at homes.	$T^{-1}$	$\left[\frac{1}{500},\frac{1}{2}\right]$	0.04375	Assumed, weighted with $f_{k_N}$ .
${\omega }_{Fk}$	Rate at which suspected freely roaming individuals are removed and isolated into the treatment centres.	$T^{-1}$	$\left[\frac{1}{500},\frac{1}{2}\right]$	0.0125	Assumed, weighted with $f_{k_N}$ .
${\alpha }_{Qk}$	Rate at which suspected quarantined individuals moved to the susceptible class $k$ .	$T^{-1}$	$\left[\frac{1}{200},\frac{1}{2}\right]$	0.04375	Assumed, weighted with $f_{Q_N}$ .
$\gamma$	Rate at which suspected quarantined individuals who are confirmed to be infected with the virus progressed to the isolated infectious individuals at the treatment centre.	$T^{-1}$	$\left[\frac{1}{200},\frac{1}{2}\right]$	0.05	Assumed, weighted with $f_{Q_N}$ .
$\beta$	Rate at which suspected quarantined individuals who are confirmed to be infected with the virus progressed to the isolated infectious individuals at homes.	$T^{-1}$	$\left[\frac{1}{200},\frac{1}{2}\right]$	0.03125	Assumed, weighted with $f_{Q_N}$ .
${\alpha }_{Rk}$	Rate at which recovered individuals lose immunity.	$T^{-1}$	$\left[2.73\times {10}^{-6},\frac{1}{120}\right]$	1/120	[24].
$a_{k,j}$ , $b_{k,j}$ , $c_{k,j}$	Rates at which susceptible, suspected freely-roaming and infectious freely-roaming individuals change their adherence status, from adherence status index $k$ to adherence status index $j$ within the respective compartments.	$T^{-1}$	Varies	Varies	Estimated
$g_k$ , ${\beta }_H$ , ${\beta }_T$	Total outflow rates for individuals who are infected/infectious symptomatic freely-roaming, symptomatic isolated at homes and at treatment centres with adherence index $k$ , $k=0,1,2,\cdots ,N$	$T^{-1}$	$\left[\frac{3}{500},\frac{1}{5}\right]$	$\frac{3}{20}$	[25] [26]
${\sigma }_{Fk}$	Rate at which infectious freely roaming individuals are removed and isolated into the treatment centres.	$T^{-1}$	$\left[0,\frac{1}{5}\right]$	0.0375	Assumed, weighted with $g_k$ .

Table 3. Model’s parameter values and their Quasi-dimensions.

Parameter	Description	Units	Range of possible values	Baseline values	References
${\rho }_{Fk}$	Rate at which infectious freely roaming individuals are removed and isolated at homes.	$T^{-1}$	$\left[0,\frac{1}{5}\right]$	0.0375	Assumed, weighted with $g_k$ .
$r_{Fk}$	Recovery rate of confirmed infected/infectious symptomatic freely roaming individuals with adherence index $k$ .	$T^{-1}$	$\left[0,\frac{1}{5}\right]$	0.0375	Assumed, weighted with $g_k$ .
${\delta }_{Fk}$	Disease related death rate for infected/infectious symptomatic freely roaming individuals with adherence index $k$ .	$T^{-1}$	$\left[0,\frac{1}{5}\right]$	0.0375	Assumed, weighted with $g_k$ .
$r_H$	Rate at which confirmed isolated symptomatic infectious individuals at homes recover from the infection.	$T^{-1}$	$\left[0,\frac{1}{5}\right]$	0.075	Assumed, weighted with ${\beta }_H$ .
$\eta$	Rate at which confirmed isolated symptomatic infectious individuals at homes become confirmed isolated symptomatic infectious individuals in treatment centres.	$T^{-1}$	$\left[0,\frac{1}{5}\right]$	0.03	Assumed, weighted with ${\beta }_H$ .
${\delta }_H$	Disease related death rates for confirmed symptomatic self-isolated individuals at homes.	$T^{-1}$	$\left[0,\frac{1}{5}\right]$	0.045	Assumed, weighted with ${\beta }_H$ .
$r_T$	Recovery rate of confirmed isolated symptomatic infectious individuals in treatment centres.	$T^{-1}$	$\left[0,\frac{1}{5}\right]$	0.09	Assumed, weighted with ${\beta }_T$ .
${\delta }_T$	Disease related death rate for confirmed symptomatic isolated individuals in treatment centres.	$T^{-1}$	$\left[0,\frac{1}{5}\right]$	0.06	Assumed, weighted with ${\beta }_T$ .
${\overline{\tau }}_{k,j}$	Entries of Hilbert-Like matrix for contacts between susceptible humans with adherence index $k$ and suspected freely-roaming humans or infected/infectious freely roaming humans with adherence index $j$ , $k,j=0,1,2,\cdots ,N$ .	1	[0, 1]	Varies	Assumed
${\phi }_{I_F}$	Effective contact rate between susceptible humans with adherence index $k$ and infected/infectious freely-roaming individuals with adherence index $j$ .	$T^{-1}$	$\left[0.002,1.12\right]$	0.18	[20] [25] [27] [28]
$C_{\ast }$	Proportions of $I_F$ , $S_F$ , $I_H$ , $I_T$ , $S_Q$ and $R_D$ humans who come in contact with $S$ humans. Here, $\ast \in \left\{I_F,S_F,I_H,I_T,S_Q,R_D\right\}$	1	[0, 1]	Varies	Assumed
${\tau }_{k,j}$	Contact rate between susceptible humans with adherence index $k$ and suspected asymptomatic freely-roaming humans with adherence index $j$ .	$T^{-1}$	$\left[0.0002,1.12\right]$	0.15	Weighted with ${\phi }_{I_F}$ .
${\phi }_{k,j}$	Contact rate between susceptible humans with adherence index $k$ and confirmed infected/infectious symptomatic freely-roaming humans with adherence index $j$ .	$T^{-1}$	$\left[0.0002,1.12\right]$	0.18	Weighted with ${\phi }_{I_F}$ .
${\epsilon }_k$	Effective contact rate between susceptible humans with adherence index $k$ and suspected quarantined asymptomatic humans.	$T^{-1}$	$\left[0.0002,1.12\right]$	0.06	Weighted with ${\phi }_{I_F}$ .
${\varphi }_k$	Effective contact rate between susceptible humans with adherence index $k$ and confirmed isolated symptomatic infected/infectious humans in the treatment centre.	$T^{-1}$	$\left[0.0002,1.12\right]$	0.09	Weighted with ${\phi }_{I_F}$ .

Table 4. Model’s parameter values and their Quasi-dimensions.

Parameter	Description	Units	Range of possible values	Baseline values	References
${\zeta }_k$	Effective contact rate between susceptible humans with adherence index $k$ and confirmed isolated symptomatic infected/infectious humans at home.	$T^{-1}$	$\left[0.0002,1.12\right]$	0.12	Weighted with ${\phi }_{I_F}$ .
$d_k$	Effective contact rate between cadavers of confirmed isolated/non-isolated symptomatic individuals and susceptible individuals with adherence index $k$ .	$T^{-1}$	$\left[0.0002,1.12\right]$	0.03	Weighted with ${\phi }_{I_F}$ .
${\eta }_3$	The environmental contamination factor, standing for the per capita rate of people who interact with the environment daily	${\mathcal{V}}^{-1}$	$\left[0,1\right]$	0.000007	[20]
${\chi }_k$	Effective contact rate between susceptible humans with adherence index $k$ and viral contaminated environment.	$T^{-1}$	$\left[0.0002,1.12\right]$	$1.26\times {10}^{-6}$	Weighted with ${\phi }_{I_F}$ .
$d$	Decay rate of viral particles from the environment.	$T^{-1}$	$\left[\frac{1}{16},1.275\right]$	0.85	[20]
$r_D$	Rate at which cadavers are removed and buried.	$T^{-1}$	$\left[\frac{1}{112},\frac{1}{17}\right]$	$\frac{129}{3808}$	[29]-[31]
$B_v$	Measures the maximum possible concentration of viral particles in the environment.	$\mathcal{V}$	$\left[10,{10}^7\right]$	30,000	Assumed
$A_v$	Measures the concentration level when the maximum concentration of viral particles in the environment is exactly half $\left(\frac{1}{2}{B}_v\right)$	$\mathcal{V}$	$\left[10,{10}^7\right]$	15,000	Assumed
${\nu }_{Fk}$	Rate of shedding viral particles into the environment by the infected/infectious symptomatic freely-roaming individuals with adherence index $k$ .	$\mathcal{V}{\mathscr{H}}^{-1}{T}^{-1}$	$\left[0,0.5\right]$	$\frac{16}{129}$	[20] [32]
$d_{\ast \ast }$	Proportions with $\ast \ast \in \left\{I_F,S_F,I_H,I_T,S_Q\right\}$	1	[0, 1]	Varies	Assumed
${\psi }_{Fk}$	Rate of shedding viral particles into the environment by the suspected quarantined individuals with adherence index $k$ .	$\mathcal{V}{\mathscr{H}}^{-1}{T}^{-1}$	$\left[0,0.5\right]$	$\frac{64}{645}$	Weighted with ${\nu }_{Fk}$ .
${\nu }_H$	Rate of shedding viral particles into the environment by confirmed symptomatic infected/infectious self isolated individuals at homes.	$\mathcal{V}{\mathscr{H}}^{-1}{T}^{-1}$	$\left[0,0.5\right]$	$\frac{48}{645}$	Weighted with ${\nu }_{Fk}$ .
${\nu }_T$	Rate of shedding viral particles into the environment by confirmed isolated symptomatic infected/infectious individuals in the treatment centres.	$\mathcal{V}{\mathscr{H}}^{-1}{T}^{-1}$	$\left[0,0.5\right]$	$\frac{32}{645}$	Weighted with ${\nu }_{Fk}$ .
${\nu }_Q$	Rate of shedding viral particles into the environment by the suspected quarantined individuals.	$\mathcal{V}{\mathscr{H}}^{-1}{T}^{-1}$	$\left[0,0.5\right]$	$\frac{16}{645}$	Weighted with ${\nu }_{Fk}$ .
$B_k$ , $k\in ℐ$	Per capita birth rate of humans per unit time	${\mathscr{H}}^{-1}{T}^{-1}$	Varies	Varies
${\lambda }_{k,j}$	A real-valued function representing the force of infection. The index $k$ indicates the flow rate from the $S_k$ class in the $S$ compartment to the $S_{Fj}$ compartment, where $k,j=0,1,2,\cdots ,N$ .	$T^{-1}$	Varies	Varies	Estimated
${\lambda }_{kQ}$	A real-valued function representing the force of infection. The index $k$ indicates the flow rate from the $S_k$ class in the $S$ compartment to the $S_Q$ compartment, where $k=0,1,2,\cdots ,N$ .	$T^{-1}$	Varies	Varies	Estimated

3. Mathematical Analyses of the Scaled Model

In this section, we perform mathematical analyses on the scaled system (29).

3.1. The Scaled Disease Free Model

In the absence of the disease, system (29)-(30) reduces to:

$\frac{\text{d}{\tilde{S}}_k}{\text{d}\tau }={\gamma }_{S_k}\left[r_k\left(1-h{\tilde{H}}_L\right){\tilde{H}}_L+\left({\displaystyle \sum }_{j=0,j\ne k}^N{\tilde{a}}_{j,k}{\tilde{S}}_j\right)-{\tilde{S}}_k\right];$(34)

$\frac{\text{d}{\tilde{H}}_L}{\text{d}\tau }=\widehat{r}\left(1-{\tilde{H}}_L\right){\tilde{H}}_L.$(35)

Here, the dimensionless parameters; ${\gamma }_{S_k},r_k,h,{\tilde{a}}_{j,k},\widehat{r}$ and $\tilde{r}$ are as defined in Equation (28). When there is no barrier measure, that is when $N=0$ , ${\tilde{H}}_L={\tilde{S}}_0$ and Equation (34) reduces to

$\frac{\text{d}{\tilde{S}}_0}{\text{d}\tau }={\gamma }_{S_0}\left[r_0\left(1-h{\tilde{S}}_0\right){\tilde{S}}_0-{\tilde{S}}_0\right]\equiv \widehat{r}\left(1-{\tilde{H}}_L\right){\tilde{H}}_L,$(36)

a single equation which can be rearranged to give Equation (35). Equation (35) can be solved directly to have, ${\tilde{S}}_0\left(\tau \right)={\tilde{H}}_L\left(\tau \right)=\frac{\varsigma }{{\text{e}}^{-\widehat{r}\tau }+\varsigma }$, $\varsigma ={\text{e}}^{c_0}$ with $c_0\in ℝ$ .

Now, as $\tau \to \infty$ , ${\tilde{H}}_L\left(\tau \right)\to 1$ . Thus, based on the scaling adopted, the equation for the scaled total living human population is always bounded, i.e., $0\le {\tilde{H}}_L\left(\tau \right)\le 1$ . Let $x_{DF}={\left({\tilde{S}}_0,{\tilde{S}}_1,\cdots ,{\tilde{S}}_N\right)}^{\text{T}}$, $f_{DF}={\left(0,0,\cdots ,{\gamma }_{S_N}{r}_N\left(1-h{\tilde{H}}_L\right){\tilde{H}}_L\right)}^{\text{T}}$ and

$A_{DF}=\left(\begin{array}{ccccc}-{\gamma }_{S_0}& {\gamma }_{S_0}{\tilde{a}}_{1,0}& {\gamma }_{S_0}{\tilde{a}}_{2,0}& \dots & {\gamma }_{S_0}{\tilde{a}}_{N,0}\\ {\gamma }_{S_1}{\tilde{a}}_{0,1}& -{\gamma }_{S_1}& {\gamma }_{S_1}{\tilde{a}}_{2,1}& \dots & {\gamma }_{S_1}{\tilde{a}}_{N,1}\\ {\gamma }_{S_2}{\tilde{a}}_{0,2}& {\gamma }_{S_2}{\tilde{a}}_{1,2}& -{\gamma }_{S_2}& \dots & {\gamma }_{S_2}{\tilde{a}}_{N,2}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\gamma }_{S_N}{\tilde{a}}_{0,N}& {\gamma }_{S_N}{\tilde{a}}_{1,N}& {\gamma }_{S_N}{\tilde{a}}_{2,N}& \dots & -{\gamma }_{S_N}\end{array}\right).$ (37)

Then system (34)-(35) can be written in the form:

$\frac{\text{d}{x}_{DF}}{\text{d}\tau }=A_{DF}{x}_{DF}+f_{DF}\left(\tau \right),x_{DF}\left(0\right)=x_0.$ (38)

Let us further define, ${\tilde{S}}^{\text{*}}={\left({\tilde{S}}_0^{\text{*}},{\tilde{S}}_1^{\text{*}},\cdots ,{\tilde{S}}_N^{\text{*}}\right)}^{\text{T}}$ and $\tilde{b}={\left(0,0,\cdots ,-r_N\left(1-h\right)\right)}^{\text{T}}$.

Theorem 1 (On the existence of steady states). The disease free system (34)-(35), has two steady states: 1) the trivial steady state; given by $x_{ss}^0=\left({\tilde{S}}^{*},{\tilde{H}}_L^{*}\right)=\left(0,0\right)$ , with $0$ a $1\times \left(N+1\right)$ zero vector and 2) the non-trivial steady state given by $x_{ss}^1=\left({\tilde{S}}^{*},{\tilde{H}}_L^{*}\right)=\left({\tilde{S}}_1^{*},1\right)$ .

Proof. At equilibrium, $\frac{\text{d}{\tilde{S}}_k}{\text{d}\tau }=\frac{\text{d}{\tilde{H}}_L}{\text{d}\tau }=0$. From Equation (35), ${\tilde{H}}_L^{\text{*}}=0$ or 1. For

${\tilde{H}}_L^{\text{*}}=0$, Equation (34) gives ${\tilde{S}}_k^{\text{*}}=0$ for all $k$ , yielding the trivial steady state. For ${\tilde{H}}_L^{\text{*}}=1$, Equation (34) gives $\tilde{A}{\tilde{S}}^{\text{*}}=\tilde{b}$, where $\tilde{A}$ has negative diagonal and non-negative off-diagonal entries. By the Gershgorin Circle Theorem [33], $\tilde{A}$ is nonsingular; hence ${\tilde{S}}^{\text{*}}={\tilde{A}}^{-1}\tilde{b}$ exists, establishing the non-trivial disease-free steady state.

Theorem 2 (Stability of the steady states).

The steady states defined by Theorem 1 have the following stability properties: (i) The trivial steady state is locally unstable for all parameters, and (ii) The non-trivial steady state is globally asymptotically stable.

Proof. Stability is determined by the eigenvalues of the Jacobian at the steady states. From the block structure of the Jacobian, the characteristic equation can be written as

$\left|J_{DF}\left(x_{ss}\right)-\lambda I\right|=\left[\widehat{r}\left(1-2{\tilde{H}}_L^{\text{*}}\right)-\lambda \right]\left|A_{DF}-\lambda I\right|=0,$(39)

where $A_{DF}$ is defined in Equation (37). At the trivial steady state $x_{ss}^0$ , $\lambda =\widehat{r}>0$ , implying instability. At the non-trivial steady state $x_{ss}^1$ , $\lambda =-\widehat{r}<0$ , while $A_{DF}$ has negative diagonal and non-negative off-diagonal entries. Hence, by the Gershgorin Circle Theorem [33] [34], all eigenvalues of $A_{DF}$ have negative real parts, establishing local asymptotic stability of $x_{ss}^1$ . Moreover, the corresponding solution remains bounded as $\tau \to \infty$ [21], completing the proof. □

3.2. The Epidemiological Model with No Barrier Measures

In this subsection, our goal is to analyse the scaled epidemiological system when there are no barrier measures (that is $N=0$ ) but other non-pharmaceutical interventions such as quarantine are still in place. From our scaling and when $N=0$ , we observe that ${\gamma }_{S_0}={\widehat{\tau }}_{0,0}=1$ , $r_0=\tilde{r}$ . System (29) now reduces to:

$\begin{array}{l}\frac{\text{d}{\tilde{S}}_0}{\text{d}\tau }=\tilde{r}\left(1-h{\tilde{H}}_L\right){\tilde{H}}_L+{\tilde{\alpha }}_{0,0}{\tilde{S}}_{F0}+{\tilde{\alpha }}_{Q0}{\tilde{S}}_Q+{\tilde{\alpha }}_{R0}{\tilde{R}}_R-a_0{\lambda }_0\left(\tilde{Y}\right){\tilde{S}}_0-{\tilde{S}}_0;\\ \frac{\text{d}{\tilde{S}}_{F0}}{\text{d}\tau }={\gamma }_{S_{F0}}\left[{\tilde{\theta }}_{0,0}{\lambda }_0\left(\tilde{Y}\right){\tilde{S}}_0-{\tilde{S}}_{F0}\right];\frac{\text{d}{\tilde{S}}_Q}{\text{d}\tau }={\gamma }_{S_Q}\left[{\widehat{\theta }}_0{\lambda }_0\left(\tilde{Y}\right){\tilde{S}}_0-{\tilde{S}}_Q\right];\\ \frac{\text{d}{\tilde{I}}_{F0}}{\text{d}\tau }={\gamma }_{I_{F0}}\left[{\tilde{\beta }}_{0,0}{\tilde{S}}_{F0}-{\tilde{I}}_{F0}\right];\frac{\text{d}{\tilde{I}}_H}{\text{d}\tau }={\gamma }_{I_H}\left[f_h\left({\tilde{S}}_Q,{\tilde{S}}_{F0},{\tilde{I}}_{F0}\right)-{\tilde{I}}_H\right];\\ \frac{\text{d}{\tilde{I}}_T}{\text{d}\tau }={\gamma }_{I_T}\left[f_t\left({\tilde{S}}_Q,{\tilde{I}}_H,{\tilde{S}}_{F0},{\tilde{I}}_{F0}\right)-{\tilde{I}}_T\right];\\ \frac{\text{d}{\tilde{R}}_R}{\text{d}\tau }={\gamma }_{R_R}\left[{\tilde{r}}_{F0}{\tilde{I}}_{F0}+{\tilde{r}}_T{\tilde{I}}_T+{\tilde{r}}_H{\tilde{I}}_H-{\tilde{R}}_R\right];\\ \frac{\text{d}{\tilde{R}}_D}{\text{d}\tau }={\gamma }_{R_D}\left[{\widehat{\delta }}_{F0}{\tilde{I}}_{F0}+{\widehat{\delta }}_T{\tilde{I}}_T+{\widehat{\delta }}_H{\tilde{I}}_H-{\tilde{R}}_D\right];\\ \frac{\text{d}{\tilde{E}}_v}{\text{d}\tau }={\gamma }_{E_v}\left[{\tilde{\psi }}_{F0}{\tilde{S}}_{F0}+{\tilde{\nu }}_{F0}{\tilde{I}}_{F0}+{\tilde{\nu }}_Q{\tilde{S}}_Q+{\tilde{\nu }}_H{\tilde{I}}_H+{\tilde{\nu }}_T{\tilde{I}}_T-{\tilde{E}}_v\right].\end{array}$(40)

The total populations satisfy the scaled equations:

$\begin{array}{l}\frac{\text{d}{\tilde{H}}_L}{\text{d}\tau }=\widehat{r}\left(1-{\tilde{H}}_L\right){\tilde{H}}_L-{\tilde{\delta }}_T{\tilde{I}}_T-{\tilde{\delta }}_H{\tilde{I}}_H-{\tilde{\delta }}_{F0}{\tilde{I}}_{F0};\\ \frac{\text{d}\tilde{H}}{\text{d}\tau }=\tilde{r}\left(1-h{\tilde{H}}_L\right){\tilde{H}}_L-\tilde{H}-{\tilde{r}}_D{\tilde{R}}_D;\end{array}$(41)

and

$\begin{array}{l}{\lambda }_0\left(\tilde{Y}\right)=\frac{{\tilde{S}}_{F0}}{\tilde{H}}+{\widehat{\phi }}_{0,0}\frac{{\tilde{I}}_{F0}}{\tilde{H}}+{\widehat{\epsilon }}_0\frac{{\tilde{S}}_Q}{\tilde{H}}+{\widehat{\varphi }}_0\frac{{\tilde{I}}_T}{\tilde{H}}+{\widehat{\zeta }}_0\frac{{\tilde{I}}_H}{\tilde{H}}+{\widehat{d}}_0\frac{{\tilde{R}}_D}{\tilde{H}}+{\widehat{\chi }}_0\left(\frac{{\tilde{E}}_v}{1+{\tilde{E}}_v}\right);\\ f_h\left({\tilde{S}}_Q,{\tilde{S}}_{F0},{\tilde{I}}_{F0}\right)=\{\begin{array}{ll}\tilde{\gamma }{\tilde{S}}_Q+{\tilde{\xi }}_{F0}{\tilde{S}}_{F0}+{\tilde{\rho }}_{F0}{\tilde{I}}_{F0},\hfill & \text{if}〚{\tilde{I}}_T〛>b_{\text{max}}\hfill \\ \tilde{\beta }{\tilde{S}}_Q+{\tilde{\omega }}_{F0}{\tilde{S}}_{F0}+{\tilde{\sigma }}_{F0}{\tilde{I}}_{F0},\hfill & \text{if}〚{\tilde{I}}_T〛\le b_{\text{max}}\hfill \end{array}\\ f_t\left({\tilde{S}}_Q,{\tilde{I}}_H,{\tilde{S}}_{F0},{\tilde{I}}_{F0}\right)=\{\begin{array}{ll}0,\hfill & \text{if}〚{\tilde{I}}_T〛>b_{\text{max}}\hfill \\ \widehat{\gamma }{\tilde{S}}_Q+\tilde{\eta }{\tilde{I}}_H+{\widehat{\mathcal{Q}}}_0,\hfill & \text{if}〚{\tilde{I}}_T〛\le b_{\text{max}}\hfill \end{array},\end{array}\}$(42)

where ${\widehat{\mathcal{Q}}}_0={\widehat{\omega }}_{F0}{\tilde{S}}_{F0}+{\widehat{\sigma }}_{F0}{\tilde{I}}_{F0}$ . The dimensionless parameters in the case when $N=0$ can be deduced from the parameter groupings in Equation (28).

3.2.1. Basic Reproduction Number of the Epidemiological Model

The basic reproduction number, ${ℛ}_0$ , measures the average number of secondary infections generated in an absolutely susceptible population by a single infectious individual throughout the period within which the individual is infectious [35]. To compute ${ℛ}_0$ using the next generation matrix approach [36], we identify the vector of new infection, $ℱ$ , and the vector of transitions, $\mathcal{V}$ , and identify ${ℛ}_0$ as the dominant eigenvalue of the next generation matrix $F{V}^{-1}$ where

$F=\left[\frac{\partial {ℱ}_i}{\partial x_j}\left(E_{dfe}^{\text{*}}\right)\right]$ , $V=\left[\frac{\partial {\mathcal{V}}_i}{\partial x_j}\left(E_{dfe}^{\text{*}}\right)\right]$ , $i,j=1,2,\cdots ,10$ . Here,

$E_{dfe}^{\text{*}}=\left(1,0,0,0,0,0,0,0,0,1\right)$ is the non-trivial disease free steady state of system (40)-(41). To proceed, we identify all state variables for the infection from system (40) as, ${\tilde{S}}_{F0}$ , ${\tilde{S}}_Q$ , ${\tilde{I}}_{F0}$ , ${\tilde{I}}_H$ , ${\tilde{I}}_T$ , ${\tilde{R}}_D$ , ${\tilde{E}}_v$ , and their corresponding equations to be:

$\begin{array}{l}\frac{\text{d}{\tilde{S}}_{F0}}{\text{d}\tau }={\gamma }_{S_{F0}}\left[{\tilde{\theta }}_{0,0}{\lambda }_0\left(\tilde{y}\right){\tilde{S}}_0-{\tilde{S}}_{F0}\right];\frac{\text{d}{\tilde{S}}_Q}{\text{d}\tau }={\gamma }_{S_Q}\left[{\widehat{\theta }}_0{\lambda }_0\left(\tilde{y}\right){\tilde{S}}_0-{\tilde{S}}_Q\right];\\ \frac{\text{d}{\tilde{I}}_{F0}}{\text{d}\tau }={\gamma }_{I_{F0}}\left[{\tilde{\beta }}_{0,0}{\tilde{S}}_{F0}-{\tilde{I}}_{F0}\right];\frac{\text{d}{\tilde{I}}_H}{\text{d}\tau }={\gamma }_{I_H}\left[f_h\left({\tilde{S}}_Q,{\tilde{S}}_{F0},{\tilde{I}}_{F0}\right)-{\tilde{I}}_H\right];\\ \frac{\text{d}{\tilde{I}}_T}{\text{d}\tau }={\gamma }_{I_T}\left[f_t\left({\tilde{S}}_Q,{\tilde{I}}_H,{\tilde{S}}_{F0},{\tilde{I}}_{F0}\right)-{\tilde{I}}_T\right];\\ \frac{\text{d}{\tilde{R}}_D}{\text{d}\tau }={\gamma }_{R_D}\left[{\widehat{\delta }}_{F0}{\tilde{I}}_{F0}+{\widehat{\delta }}_T{\tilde{I}}_T+{\widehat{\delta }}_H{\tilde{I}}_H-{\tilde{R}}_D\right];\\ \frac{\text{d}{\tilde{E}}_v}{\text{d}\tau }={\gamma }_{E_v}\left[{\tilde{\psi }}_{F0}{\tilde{S}}_{F0}+{\tilde{\nu }}_{F0}{\tilde{I}}_{F0}+{\tilde{\nu }}_Q{\tilde{S}}_Q+{\tilde{\nu }}_H{\tilde{I}}_H+{\tilde{\nu }}_T{\tilde{I}}_T-{\tilde{E}}_v\right].\end{array}\}$(43)

It is known that the value of ${ℛ}_0$ obtained from the next generation approach is sensitive to the manner in which we interpret transitions and new infections [37]. Here, we can interrogate possibilities concerning the viral particles and the environment: (i) Environmental contamination seen as a transition, (ii) Environmental contamination seen as new infections.

CASE I: The contaminated environment as a transition. This pathway is common with viral infections whereby the viral particles’ population cannot maintain itself through growth in the environment. So, the shedding of viral particle terms by infectious individuals are placed in the transition state rather than the new infection state, giving:

${ℛ}_0^I=\{\begin{array}{ll}{ℛ}_{0_a}^I={ℛ}_{0_a}^h+{ℛ}_{0_a}^e;\hfill & \text{if}〚{\tilde{I}}_T〛\le b_{\text{max}}\hfill \\ {ℛ}_{0_b}^I={ℛ}_{0_b}^h+{ℛ}_{0_b}^e;\hfill & \text{if}〚{\tilde{I}}_T〛>b_{\text{max}}\hfill \end{array}.$(44)

For; ${ℛ}_{0_a}^h={\tilde{\theta }}_{0,0}{M}_1+{\widehat{\theta }}_0{M}_2$ ; ${ℛ}_{0_a}^e={\widehat{\chi }}_0\left({\tilde{\theta }}_{0,0}{M}_3+{\widehat{\theta }}_0{M}_4\right)$ , and ${ℛ}_{0_b}^h={\tilde{\theta }}_{0,0}{T}_1+{\widehat{\theta }}_0{T}_2$ ;${ℛ}_{0_b}^e={\widehat{\chi }}_0\left({\tilde{\theta }}_{0,0}{T}_3+{\widehat{\theta }}_0{T}_4\right)$ . Where:

$\begin{array}{l}{M}_1=1+{\widehat{\phi }}_{0,0}{\tilde{\beta }}_{0,0}+{\widehat{\zeta }}_0\left({\tilde{\beta }}_{0,0}{\tilde{\sigma }}_{F0}+{\tilde{\omega }}_{F0}\right)+{\widehat{\varphi }}_0\left({\tilde{\beta }}_{0,0}{\widehat{\sigma }}_{F0}+{\widehat{\omega }}_{F0}+\tilde{\eta }\left({\tilde{\beta }}_{0,0}{\tilde{\sigma }}_{F0}+{\tilde{\omega }}_{F0}\right)\right)\\ +{\widehat{d}}_0\left({\tilde{\beta }}_{0,0}\left({\widehat{\delta }}_{F0}+{\tilde{\sigma }}_{F0}{\widehat{\delta }}_H+{\widehat{\delta }}_T{\widehat{\sigma }}_{F0}+\tilde{\eta }{\widehat{\delta }}_T{\tilde{\sigma }}_{F0}\right)+{\tilde{\omega }}_{F0}\left({\widehat{\delta }}_H+\tilde{\eta }{\widehat{\delta }}_T\right)+{\widehat{\delta }}_T{\widehat{\omega }}_{F0}\right);\\ M_2={\widehat{\epsilon }}_0+\tilde{\beta }{\widehat{\zeta }}_0+{\widehat{\varphi }}_0\left(\widehat{\gamma }+\tilde{\eta }\tilde{\beta }\right)+{\widehat{d}}_0\left(\tilde{\beta }{\widehat{\delta }}_H+\widehat{\gamma }{\widehat{\delta }}_T+\tilde{\eta }\tilde{\beta }{\widehat{\delta }}_T\right);\\ M_3={\tilde{\psi }}_{F0}+{\tilde{\nu }}_{F0}{\tilde{\beta }}_{0,0}+{\tilde{\nu }}_T\left({\widehat{\omega }}_{F0}+{\widehat{\sigma }}_{F0}{\tilde{\beta }}_{0,0}+\tilde{\eta }\left({\tilde{\beta }}_{0,0}{\tilde{\sigma }}_{F0}+{\tilde{\omega }}_{F0}\right)\right)\\ +{\tilde{\nu }}_H\left({\tilde{\omega }}_{F0}+{\tilde{\sigma }}_{F0}{\tilde{\beta }}_{0,0}\right);\\ M_4=\tilde{\beta }{\tilde{\nu }}_H+{\tilde{\nu }}_Q+{\tilde{\nu }}_T\left(\tilde{\eta }\tilde{\beta }+\widehat{\gamma }\right);\\ T_1=1+{\widehat{\phi }}_{0,0}{\tilde{\beta }}_{0,0}+{\widehat{\zeta }}_0\left({\tilde{\beta }}_{0,0}{\tilde{\rho }}_{F0}+{\tilde{\xi }}_{F0}\right)+{\widehat{d}}_0\left({\tilde{\beta }}_{0,0}\left({\widehat{\delta }}_{F0}+{\tilde{\rho }}_{F0}{\widehat{\delta }}_H\right)+{\tilde{\xi }}_{F0}{\widehat{\delta }}_H\right);\\ T_2={\widehat{\epsilon }}_0+\tilde{\gamma }{\widehat{\zeta }}_0+{\widehat{d}}_0\tilde{\gamma }{\widehat{\delta }}_H;T_3={\tilde{\psi }}_{F0}+{\tilde{\nu }}_{F0}{\tilde{\beta }}_{0,0}+{\tilde{\nu }}_H\left({\tilde{\xi }}_{F0}+{\tilde{\rho }}_{F0}{\tilde{\beta }}_{0,0}\right);T_4=\tilde{\gamma }{\tilde{\nu }}_H+{\tilde{\nu }}_Q.\end{array}\}$(45)

The quantities ${ℛ}_{0_i}^h$ and ${ℛ}_{0_i}^e$ , $i\in \left\{a,b\right\}$ , defined in Equation (44), correspond, respectively, to the average number of secondary infections through human-to-human and environment-to-human transmissions caused by one infectious individual in its infectious lifetime. From the expression of ${ℛ}_0^I$ defined in Equation (44), the contribution of each transmission pathway in a viral disease outbreak is separable. Hence, the control efforts can be focused on infectious individuals or the contaminated environment depending on the amount of effort required to reduce the sum of ${ℛ}_{0_i}^h$ and ${ℛ}_{0_i}^e$ , $i\in \left\{a,b\right\}$ , to less than unity.

CASE II: Contaminated environment as a reservoir of infection. By this we mean the environment is assumed to act as a reservoir, where secondary viral particles are added into the environment solely through viral particle shedding by infectious host. By this, we mean, infectious individuals shed viral particles into the environment, these viral particles exists independently, they accumulate and decay over time. Also, these viral particles go on to infect new individuals. Intuitively, infection can happen even after infectious individuals leave, through the environment. Hence, all the shedding rates are placed in the $ℱ$ matrix. Then, after some calculations, we have:

${ℛ}_0^{II}=\{\begin{array}{l}\begin{array}{ll}{ℛ}_{0_a}^{II}=\frac{{ℛ}_{0_a}^h+\sqrt{{\left({ℛ}_{0_a}^h\right)}^2+4{ℛ}_{0_a}^e}}{2};\hfill & \text{if}〚{\tilde{I}}_T〛\le b_{\text{max}},\hfill \\ {ℛ}_{0_b}^{II}=\frac{{ℛ}_{0_b}^h+\sqrt{{\left({ℛ}_{0_b}^h\right)}^2+4{ℛ}_{0_b}^e}}{2};\hfill & \text{if}〚{\tilde{I}}_T〛>b_{\text{max}}.\hfill \end{array}\hfill \end{array}$(46)

Here, ${ℛ}_{0_i}^h$ and ${ℛ}_{0_i}^e$ , $i\in \left\{a,b\right\}$ , retain their definitions above. Thus, ${ℛ}_0\in \left\{{ℛ}_0^I,{ℛ}_0^{II}\right\}$ separates human-to-human and environmental transmission. Hence, ${ℛ}_0$ allows the contributions of quarantined individuals, non-quarantined individuals, human transmission, environmental transmission, or all pathways to be examined separately. We therefore obtain the following result for ${ℛ}_0^{II}$ in Equation (46).

Lemma 3. Let ${ℛ}_0^h\ge 0$ and ${ℛ}_0^e\ge 0$ be given. Define:

${ℛ}_0=\frac{{ℛ}_0^h+\sqrt{{\left({ℛ}_0^h\right)}^2+4{ℛ}_0^e}}{2}.$ (47)

(i) If ${ℛ}_0^h={ℛ}_0^e=0$ , then ${ℛ}_0=0$ .

(ii) If ${ℛ}_0^h=0$ , then ${ℛ}_0=\sqrt{{ℛ}_0^e}$ and if ${ℛ}_0^e=0$ , then ${ℛ}_0={ℛ}_0^h$ . Thus, $0<{ℛ}_0\le 1\iff 0<{ℛ}_0^e\le 1$ , and $0<{ℛ}_0\le 1\iff 0<{ℛ}_0^h\le 1$ .

(iii) $0<{ℛ}_0\le 1\iff 0<{ℛ}_0^h+{ℛ}_0^e\le 1$ and ${ℛ}_0>1\iff {ℛ}_0^h+{ℛ}_0^e>1$ .

Proof. The proofs of (i) and (ii) are straightforward and are therefore omitted. We prove only the first part of (iii); the second part follows by a similar argument. Suppose ${ℛ}_0^h>0$ , ${ℛ}_0^e>0$ and rationalizing ${ℛ}_0$ defined in Equation (47) we get

$\begin{array}{l}0<{ℛ}_0\le 1\iff 0<\frac{2{ℛ}_0^e}{\sqrt{{\left({ℛ}_0^h\right)}^2+4{ℛ}_0^e}-{ℛ}_0^h}\le 1\\ \iff 0<4{\left({ℛ}_0^e\right)}^2+4{ℛ}_0^e{ℛ}_0^h\le 4{ℛ}_0^e\\ \iff 0<{ℛ}_0^h+{ℛ}_0^e\le 1.\end{array}$

□

Theorem 4. Let ${ℛ}_0^I$ and ${ℛ}_0^{II}$ be given as in Equations (44) and (46) respectively. Then:

(i) ${ℛ}_0^I=1\iff {ℛ}_0^{II}=1$ .

(ii) ${ℛ}_0^I<{ℛ}_0^{II}$ whenever ${ℛ}_0^I<1$ and ${ℛ}_0^I>{ℛ}_0^{II}$ whenever ${ℛ}_0^I>1$ .

(iii) ${ℛ}_0^I={ℛ}_0^{II}={ℛ}_0^h$ whenever ${ℛ}_0^e=0$ , and ${ℛ}_0^I>{ℛ}_0^{II}$ whenever ${ℛ}_0^h=0$ and ${ℛ}_0^e>1$ .

Proof. We will prove only the first point here as the second point has been proven already in Lemma 3 and the third point is straightforward.

$\begin{array}{l}{ℛ}_0^I=1\iff {ℛ}_0^h+{ℛ}_0^e=1\iff {ℛ}_0^h=1-{ℛ}_0^e\iff {\left({ℛ}_0^h\right)}^2+4{ℛ}_0^e={\left(1+{ℛ}_0^e\right)}^2\\ \iff {ℛ}_0^{II}=1.\end{array}$

□

Remark. Throughout this article, we shall continue with ${ℛ}_0^{II}$ defined by Equation (46) and simply refer to it as ${ℛ}_0$ .

3.2.2. Existence and Stability of Steady States

We consider only realistic solutions in the sense of Definition 2.1.

Theorem 5 (Existence of steady state solutions). System (40) possesses:

(i) A trivial disease free steady state: $E_0^{\text{*}}=0=\left(0,0,0,0,0,0,0,0,0,0\right)$ , which exists whenever ${\tilde{H}}_L^{\text{*}}=0$;

(ii) a non-trivial disease free steady state: $E_{dfe}^{\text{*}}=\left(1,0,0,0,0,0,0,0,0,1\right)$ , which exists whenever ${\tilde{H}}_L^{\text{*}}=1$; and

(iii) an endemic steady state: $E_{ee}^{\text{*}}=\left({\tilde{S}}_0^{\text{*}},{\tilde{S}}_{F0}^{\text{*}},{\tilde{S}}_Q^{\text{*}},{\tilde{I}}_{F0}^{\text{*}},{\tilde{I}}_H^{\text{*}},{\tilde{I}}_T^{\text{*}},{\tilde{R}}_R^{\text{*}},{\tilde{R}}_D^{\text{*}},{\tilde{E}}_v^{\text{*}},{\tilde{H}}^{\text{*}}\right)$, which exists whenever ${\tilde{H}}_L^{\text{*}}\in \left(0,1\right)$ where ${\tilde{H}}_L^{\text{*}}$ satisfies:

$P\left({\tilde{H}}_L^{\text{*}}\right)=q_3{\tilde{H}}_L^{\text{*}3}+q_2{\tilde{H}}_L^{\text{*}2}+q_1{\tilde{H}}_L^{\text{*}}+q_0=0,\text{for}〚{\tilde{I}}_T〛\le b_{\text{max}};$ or(48)

$P_3\left({\tilde{H}}_L^{\text{*}}\right)=p_3{\tilde{H}}_L^{\text{*}3}+p_2{\tilde{H}}_L^{\text{*}2}+p_1{\tilde{H}}_L^{\text{*}}+p_0=0,\text{for}〚{\tilde{I}}_T〛>b_{\text{max}}.$(49)

The coefficients $q_3,q_2,q_1,q_0$ are as defined in Equation (64) and $p_3,p_2,p_1,p_0$ are as defined in Equation (65).

Proof. The steady states of system (40) are solutions of the algebraic system of equations obtained by setting the right hand sides to zero. Let $x^{\text{*}}=\left({\tilde{S}}_0^{\text{*}},{\tilde{S}}_{F0}^{\text{*}},{\tilde{S}}_Q^{\text{*}},{\tilde{I}}_{F0}^{\text{*}},{\tilde{I}}_H^{\text{*}},{\tilde{I}}_T^{\text{*}},{\tilde{R}}_R^{\text{*}},{\tilde{R}}_D^{\text{*}},{\tilde{E}}_v^{\text{*}},{\tilde{H}}^{\text{*}}\right)$ be a steady state solution of system (40), we get;

$\tilde{r}\left(1-h{\tilde{H}}_L^{\text{*}}\right){\tilde{H}}_L^{\text{*}}+{\tilde{\alpha }}_{0,0}{\tilde{S}}_{F0}^{\text{*}}+{\tilde{\alpha }}_{Q0}{\tilde{S}}_Q^{\text{*}}+{\tilde{\alpha }}_{R0}{\tilde{R}}_R^{\text{*}}-a_0{\lambda }_0\left({\tilde{Y}}^{\text{*}}\right){\tilde{S}}_0^{\text{*}}-{\tilde{S}}_0^{\text{*}}=0;$(50)

${\tilde{S}}_{F0}^{\text{*}}={\tilde{\theta }}_{0,0}{\lambda }_0\left({\tilde{Y}}^{\text{*}}\right){\tilde{S}}_0^{\text{*}};{\tilde{S}}_Q^{\text{*}}={\widehat{\theta }}_0{\lambda }_0\left({\tilde{Y}}^{\text{*}}\right){\tilde{S}}_0^{\text{*}};{\tilde{I}}_{F0}^{\text{*}}={\tilde{\beta }}_{0,0}{\tilde{S}}_{F0}^{\text{*}};$(51)

${\tilde{I}}_H^{\text{*}}=f_h^{\text{*}}\left({\tilde{S}}_Q^{\text{*}},{\tilde{S}}_{F0}^{\text{*}},{\tilde{I}}_{F0}^{\text{*}}\right);{\tilde{I}}_T^{\text{*}}=f_t^{\text{*}}\left({\tilde{S}}_Q^{\text{*}},{\tilde{I}}_H^{\text{*}},{\tilde{S}}_{F0}^{\text{*}},{\tilde{I}}_{F0}^{\text{*}}\right);$(52)

${\tilde{R}}_R^{\text{*}}={\tilde{r}}_{F0}{\tilde{I}}_{F0}^{\text{*}}+{\tilde{r}}_T{\tilde{I}}_T^{\text{*}}+{\tilde{r}}_H{\tilde{I}}_H^{\text{*}};{\tilde{R}}_D^{\text{*}}={\widehat{\delta }}_{F0}{\tilde{I}}_{F0}^{\text{*}}+{\widehat{\delta }}_T{\tilde{I}}_T^{\text{*}}+{\widehat{\delta }}_H{\tilde{I}}_H^{\text{*}};$(53)

${\tilde{\psi }}_{F0}{\tilde{S}}_{F0}^{\text{*}}+{\tilde{\nu }}_{F0}{\tilde{I}}_{F0}^{\text{*}}+{\tilde{\nu }}_Q{\tilde{S}}_Q^{\text{*}}+{\tilde{\nu }}_H{\tilde{I}}_H^{\text{*}}+{\tilde{\nu }}_T{\tilde{I}}_T^{\text{*}}-{\tilde{E}}_v^{\text{*}}=0;$(54)

$\tilde{r}\left(1-h{\tilde{H}}_L^{\text{*}}\right){\tilde{H}}_L^{\text{*}}-{\tilde{H}}^{\text{*}}-{\tilde{r}}_D{\tilde{R}}_D^{\text{*}}=0;$(55)

$\widehat{r}\left(1-{\tilde{H}}_L^{\text{*}}\right){\tilde{H}}_L^{\text{*}}-{\tilde{\delta }}_T{\tilde{I}}_T^{\text{*}}-{\tilde{\delta }}_H{\tilde{I}}_H^{\text{*}}-{\tilde{\delta }}_{F0}{\tilde{I}}_{F0}^{\text{*}}=0$(56)

with

${\lambda }_0\left({\tilde{Y}}^{\text{*}}\right)=\frac{{\tilde{S}}_{F0}^{\text{*}}}{{\tilde{H}}^{\text{*}}}+{\widehat{\phi }}_{0,0}\frac{{\tilde{I}}_{F0}^{\text{*}}}{{\tilde{H}}^{\text{*}}}+{\widehat{\epsilon }}_0\frac{{\tilde{S}}_Q^{\text{*}}}{{\tilde{H}}^{\text{*}}}+{\widehat{\varphi }}_0\frac{{\tilde{I}}_T^{\text{*}}}{{\tilde{H}}^{\text{*}}}+{\widehat{\zeta }}_0\frac{{\tilde{I}}_H^{\text{*}}}{{\tilde{H}}^{\text{*}}}+{\widehat{d}}_0\frac{{\tilde{R}}_D^{\text{*}}}{{\tilde{H}}^{\text{*}}}+{\widehat{\chi }}_0\left(\frac{{\tilde{E}}_v^{\text{*}}}{1+{\tilde{E}}_v^{\text{*}}}\right);$(57)

$f_h^{\text{*}}\left({\tilde{S}}_Q^{\text{*}},{\tilde{S}}_{F0}^{\text{*}},{\tilde{I}}_{F0}^{\text{*}}\right)=\{\begin{array}{ll}\tilde{\gamma }{\tilde{S}}_Q^{\text{*}}+{\tilde{\xi }}_{F0}{\tilde{S}}_{F0}^{\text{*}}+{\tilde{\rho }}_{F0}{\tilde{I}}_{F0}^{\text{*}},\hfill & \text{if}〚{\tilde{I}}_T〛>b_{\text{max}}\hfill \\ \tilde{\beta }{\tilde{S}}_Q^{\text{*}}+{\tilde{\omega }}_{F0}{\tilde{S}}_{F0}^{\text{*}}+{\tilde{\sigma }}_{F0}{\tilde{I}}_{F0}^{\text{*}},\hfill & \text{if}〚{\tilde{I}}_T〛\le b_{\text{max}}\hfill \end{array};$(58)

and

$f_t^{\text{*}}\left({\tilde{S}}_Q^{\text{*}},{\tilde{I}}_H^{\text{*}},{\tilde{S}}_{F0}^{\text{*}},{\tilde{I}}_{F0}^{\text{*}}\right)=\{\begin{array}{ll}0,\hfill & \text{if}〚{\tilde{I}}_T〛>b_{\text{max}}\hfill \\ \widehat{\gamma }{\tilde{S}}_Q^{\text{*}}+\tilde{\eta }{\tilde{I}}_H^{\text{*}}+{\widehat{\mathcal{Q}}}_0^{\text{*}},\hfill & \text{if}〚{\tilde{I}}_T〛\le b_{\text{max}}\hfill \end{array}.$(59)

From (51),

${\lambda }_0\left({\tilde{Y}}^{\text{*}}\right){\tilde{S}}_0^{\text{*}}=\frac{1}{{\tilde{\theta }}_{0,0}}{\tilde{S}}_{F0}^{\text{*}}.$(60)

We now express ${\tilde{S}}_Q^{\text{*}}$, ${\tilde{I}}_{F0}^{\text{*}}$, ${\tilde{I}}_H^{\text{*}}$, ${\tilde{I}}_T^{\text{*}}$, ${\tilde{R}}_R^{\text{*}}$, ${\tilde{R}}_D^{\text{*}}$ and ${\tilde{E}}_v^{\text{*}}$ as a function of ${\tilde{S}}_{F0}^{\text{*}}$ and consider the following cases.

Case 1

Let us now consider the case where $〚{\tilde{I}}_T〛\le b_{\text{max}}$. Here, ${\tilde{S}}_Q^{\text{*}}\left({\tilde{S}}_{F0}^{\text{*}}\right)=\left(\frac{{\widehat{\theta }}_0}{{\tilde{\theta }}_{0,0}}\right){\tilde{S}}_{F0}^{\text{*}}$; ${\tilde{I}}_{F0}^{\text{*}}\left({\tilde{S}}_{F0}^{\text{*}}\right)={\tilde{\beta }}_{0,0}{\tilde{S}}_{F0}^{\text{*}}$; ${\tilde{I}}_H^{\text{*}}\left({\tilde{S}}_{F0}^{\text{*}}\right)=A_1{\tilde{S}}_{F0}^{\text{*}}$; ${\tilde{I}}_T^{\text{*}}\left({\tilde{S}}_{F0}^{\text{*}}\right)=A_2{\tilde{S}}_{F0}^{\text{*}}$; ${\tilde{R}}_R^{\text{*}}\left({\tilde{S}}_{F0}^{\text{*}}\right)=A_3{\tilde{S}}_{F0}^{\text{*}}$;

${\tilde{R}}_D^{\text{*}}\left({\tilde{S}}_{F0}^{\text{*}}\right)=A_4{\tilde{S}}_{F0}^{\text{*}}$ and ${\tilde{E}}_v^{\text{*}}\left({\tilde{S}}_{F0}^{\text{*}}\right)=A_5{\tilde{S}}_{F0}^{\text{*}}$. Expressing ${\tilde{S}}_{F0}^{\text{*}}$ in Equation (56), ${\tilde{S}}_0^{\text{*}}$

in Equation (50), ${\tilde{H}}^{\text{*}}$ in Equation (55) each in terms of ${\tilde{H}}_L^{\text{*}}$ and substituting appropriate expressions into Equation (57), we get;

$\begin{array}{l}{\tilde{S}}_{F0}^{\text{*}}\left({\tilde{H}}_L^{\text{*}}\right)=A_8\left(1-{\tilde{H}}_L^{\text{*}}\right){\tilde{H}}_L^{\text{*}};\\ {\tilde{S}}_0^{\text{*}}\left({\tilde{H}}_L^{\text{*}}\right)=\left[A_9\left(1-{\tilde{H}}_L^{\text{*}}\right)+\tilde{r}\left(1-h{\tilde{H}}_L^{\text{*}}\right)\right]{\tilde{H}}_L^{\text{*}};\\ {\tilde{H}}^{\text{*}}\left({\tilde{H}}_L^{\text{*}}\right)=\left[\tilde{r}\left(1-h{\tilde{H}}_L^{\text{*}}\right)-{\tilde{r}}_D{A}_4{A}_8\left(1-{\tilde{H}}_L^{\text{*}}\right)\right]{\tilde{H}}_L^{\text{*}},\\ {\lambda }_0\left({\tilde{Y}}^{\text{*}}\right)\left({\tilde{S}}_{F0}^{\text{*}},{\tilde{H}}^{\text{*}}\right)=\left[\frac{A_7}{{\tilde{H}}^{\text{*}}}+\frac{{\widehat{\chi }}_0{A}_5}{1+A_5{\tilde{S}}_{F0}^{\text{*}}}\right]{\tilde{S}}_{F0}^{\text{*}},\end{array}\}$(61)

for

$\begin{array}{l}{A}_1=\frac{\tilde{\beta }{\widehat{\theta }}_0}{{\tilde{\theta }}_{0,0}}+{\tilde{\omega }}_{F0}+{\tilde{\sigma }}_{F0}{\tilde{\beta }}_{0,0};A_2=\frac{\widehat{\gamma }{\widehat{\theta }}_0}{{\tilde{\theta }}_{0,0}}+\tilde{\eta }{A}_1+{\widehat{\omega }}_{F0}+{\widehat{\sigma }}_{F0}{\tilde{\beta }}_{0,0};\\ A_3={\tilde{r}}_{F0}{\tilde{\beta }}_{0,0}+{\tilde{r}}_T{A}_2+{\tilde{r}}_H{A}_1;A_4={\widehat{\delta }}_{F0}{\tilde{\beta }}_{0,0}+{\widehat{\delta }}_T{A}_2+{\widehat{\delta }}_H{A}_1;\\ A_5=M_3+\frac{{\widehat{\theta }}_0}{{\tilde{\theta }}_{0,0}}{M}_4;A_6={\tilde{\alpha }}_{0,0}+{\tilde{\alpha }}_{R0}{A}_3+\frac{{\tilde{\alpha }}_{Q0}{\widehat{\theta }}_0}{{\tilde{\theta }}_{0,0}};\\ A_7=M_1+\frac{{\widehat{\theta }}_0}{{\tilde{\theta }}_{0,0}}{M}_2;A_8=\frac{\widehat{r}}{{\tilde{\delta }}_T{A}_2+{\tilde{\delta }}_H{A}_1+{\tilde{\delta }}_{F0}{\tilde{\beta }}_{0,0}};\\ A_9=\left(A_6-\frac{a_0}{{\tilde{\theta }}_{0,0}}\right)A_8=\frac{A_8}{{\tilde{\theta }}_{0,0}}\left({\tilde{\theta }}_{0,0}{A}_6-a_0\right).\end{array}\}$(62)

where $M_1,M_2,M_3$ and $M_4$ are as defined in Equation (45). Now, using Equations (60) and (61), we have that, either

${\tilde{S}}_{F0}^{\text{*}}\left({\tilde{H}}_L^{\text{*}}\right)=0\text{or}\left[\frac{A_7}{{\tilde{H}}^{\text{*}}\left({\tilde{H}}_L^{\text{*}}\right)}+\frac{{\widehat{\chi }}_0{A}_5}{1+A_5{\tilde{S}}_{F0}^{\text{*}}\left({\tilde{H}}_L^{\text{*}}\right)}\right]{\tilde{S}}_0^{\text{*}}\left({\tilde{H}}_L^{\text{*}}\right)-\frac{1}{{\tilde{\theta }}_{0,0}}=0.$(63)

Therefore, ${\tilde{S}}_{F0}^{\text{*}}\left({\tilde{H}}_L^{\text{*}}\right)=0$ implies that, $A_8\left(1-{\tilde{H}}_L^{\text{*}}\right){\tilde{H}}_L^{\text{*}}=0$. Which gives ${\tilde{H}}_L^{\text{*}}=0$ or ${\tilde{H}}_L^{\text{*}}=1$. Hence, we obtain the trivial and non-trivial disease free equilibria given by $E_0^{\text{*}}$ and $E_{dfe}^{\text{*}}$ respectively. Replacing ${\tilde{H}}^{\text{*}}\left({\tilde{H}}_L^{\text{*}}\right)$, ${\tilde{S}}_0^{\text{*}}\left({\tilde{H}}_L^{\text{*}}\right)$ and ${\tilde{S}}_{F0}^{\text{*}}\left({\tilde{H}}_L^{\text{*}}\right)$ in Equation (63) and simplifying, we obtain the third order degree polynomial, $P\left({\tilde{H}}_L^{\text{*}}\right)$ as defined in Equation (48) where:

$\begin{array}{l}{q}_3=A_5\left(A_9+\tilde{r}h\right)\left[\left(A_7{A}_8+{\widehat{\chi }}_0\tilde{r}h\right)-{\widehat{\chi }}_0{\tilde{r}}_D{A}_4{A}_8\right]+\frac{A_5{A}_8}{{\tilde{\theta }}_{0,0}}\left({\tilde{r}}_D{A}_4{A}_8-\tilde{r}h\right);\\ q_2=-A_5\left(A_9+\tilde{r}\right)\left[\left(A_7{A}_8+{\widehat{\chi }}_0\tilde{r}h\right)-{\widehat{\chi }}_0{\tilde{r}}_D{A}_4{A}_8\right]\\ -A_5\left(A_9+\tilde{r}h\right)\left[\left(A_7{A}_8+{\widehat{\chi }}_0\tilde{r}\right)-{\widehat{\chi }}_0{\tilde{r}}_D{A}_4{A}_8\right]\\ -\frac{A_5{A}_8}{{\tilde{\theta }}_{0,0}}\left[2{\tilde{r}}_D{A}_4{A}_8-\tilde{r}\left(1+h\right)\right];\\ q_1=A_5\left(A_9+\tilde{r}\right)\left[\left(A_7{A}_8+{\widehat{\chi }}_0\tilde{r}\right)-{\widehat{\chi }}_0{\tilde{r}}_D{A}_4{A}_8\right]-A_7\left(A_9+\tilde{r}h\right)\\ -\frac{1}{{\tilde{\theta }}_{0,0}}\left[A_5{A}_8\left(\tilde{r}-{\tilde{r}}_D{A}_4{A}_8\right)+\left({\tilde{r}}_D{A}_4{A}_8-\tilde{r}h\right)\right];\\ q_0=A_7\left(A_9+\tilde{r}\right)+\frac{1}{{\tilde{\theta }}_{0,0}}\left({\tilde{r}}_D{A}_4{A}_8-\tilde{r}\right).\end{array}\}$(64)

Case 2

This is the case when $〚{\tilde{I}}_T〛>b_{\text{max}}$. Following similar steps as in the first case, we obtain the cubic polynomial in Equation (49), where:

$\begin{array}{l}{p}_3={\tilde{\theta }}_{0,0}{N}_4{N}_5{N}_7\left(N_6+\tilde{r}h\right)+N_4{N}_5\left({\tilde{r}}_D{N}_3{N}_5-\tilde{r}h\right)\\ -{\tilde{\theta }}_{0,0}{\widehat{\chi }}_0{N}_4\left(N_6+\tilde{r}h\right)\times \left({\tilde{r}}_D{N}_3{N}_5-\tilde{r}h\right);\\ p_2={\tilde{\theta }}_{0,0}{\widehat{\chi }}_0{N}_4\left(N_6+\tilde{r}\right)\left({\tilde{r}}_D{N}_3{N}_5-\tilde{r}h\right)+N_4{N}_5\left(\tilde{r}-{\tilde{r}}_D{N}_3{N}_5\right)\\ -{\tilde{\theta }}_{0,0}{\widehat{\chi }}_0{N}_4\left(N_6+\tilde{r}h\right)\left(\tilde{r}-{\tilde{r}}_D{N}_3{N}_5\right)-{\tilde{\theta }}_{0,0}{N}_4{N}_5{N}_7\left(N_6+\tilde{r}\right)\\ -{\tilde{\theta }}_{0,0}{N}_4{N}_5{N}_7\left(N_6+\tilde{r}h\right)-N_4{N}_5\left({\tilde{r}}_D{N}_3{N}_5-\tilde{r}h\right);\\ p_1={\tilde{\theta }}_{0,0}{N}_4{N}_5{N}_7\left(N_6+\tilde{r}\right)+{\tilde{\theta }}_{0,0}{\widehat{\chi }}_0{N}_4\left(N_6+\tilde{r}\right)\left(\tilde{r}-{\tilde{r}}_D{N}_3{N}_5\right)\\ -{\tilde{\theta }}_{0,0}{N}_7\left(N_6+\tilde{r}h\right)+\left(\tilde{r}h-{\tilde{r}}_D{N}_3{N}_5\right)-N_4{N}_5\left(\tilde{r}-{\tilde{r}}_D{N}_3{N}_5\right);\\ p_0={\tilde{\theta }}_{0,0}{N}_7\left(N_6+\tilde{r}\right)+\left({\tilde{r}}_D{N}_3{N}_5-\tilde{r}\right).\end{array}\}$(65)

For

$\begin{array}{l}{N}_1=\frac{\tilde{\gamma }{\widehat{\theta }}_0}{{\tilde{\theta }}_{0,0}}+{\tilde{\xi }}_{F0}+{\tilde{\rho }}_{F0}{\tilde{\beta }}_{0,0},N_2={\tilde{r}}_{F0}{\tilde{\beta }}_{0,0}+{\tilde{r}}_H{N}_1,\\ N_3={\widehat{\delta }}_{F0}{\tilde{\beta }}_{0,0}+{\widehat{\delta }}_H{N}_1,N_4={\tilde{\psi }}_{F0}+{\tilde{\nu }}_{F0}{\tilde{\beta }}_{0,0}+\frac{{\tilde{\nu }}_Q{\widehat{\theta }}_0}{{\tilde{\theta }}_{0,0}}+{\tilde{\nu }}_H{N}_1\\ N_5=\frac{\widehat{r}}{{\tilde{\delta }}_H{N}_1+{\tilde{\delta }}_{F0}{\tilde{\beta }}_{0,0}};N_6=\left(N_8-\frac{a_0}{{\tilde{\theta }}_{0,0}}\right)N_5=\frac{N_5}{{\tilde{\theta }}_{0,0}}\left({\tilde{\theta }}_{0,0}{N}_8-a_0\right);\\ N_7=1+{\widehat{\phi }}_{0,0}{\tilde{\beta }}_{0,0}+{\widehat{\epsilon }}_0\left(\frac{{\widehat{\theta }}_0}{{\tilde{\theta }}_{0,0}}\right)+{\widehat{\zeta }}_0{N}_1+{\widehat{d}}_0{N}_3;N_8={\tilde{\alpha }}_{0,0}+{\tilde{\alpha }}_{R0}{N}_2+\frac{{\tilde{\alpha }}_{Q0}{\widehat{\theta }}_0}{{\tilde{\theta }}_{0,0}}.\end{array}\}$(66)

□

Theorem 6 (On the existence of an endemic steady state). Let ${ℛ}_0>1$ . If,

$\begin{array}{l}{ℛ}_{0_a}^h\left(\tilde{r}+A_6{A}_8\right)+{\tilde{r}}_D{A}_4{A}_8\le \left(\tilde{r}+a_0{A}_7{A}_8\right)\\ =\frac{1}{{\tilde{\theta }}_{0,0}}\left({\tilde{\theta }}_{0,0}\tilde{r}+a_0{A}_8{ℛ}_{0_a}^h\right),\text{for}〚{\tilde{I}}_T〛\le b_{\text{max}};\end{array}$(67)

$\begin{array}{l}{ℛ}_{0_b}^h\left(\tilde{r}+N_5{N}_8\right)+{\tilde{r}}_D{N}_3{N}_5\le \left(\tilde{r}+a_0{N}_5{N}_7\right)\\ =\frac{1}{{\tilde{\theta }}_{0,0}}\left({\tilde{\theta }}_{0,0}\tilde{r}+a_0{N}_5{ℛ}_{0_b}^h\right),\text{for}〚{\tilde{I}}_T〛>b_{\text{max}};\end{array}$(68)

then there exists at least one real positive solution of ${\tilde{H}}_L^{\text{*}}$ in Equation (48) and Equation (49) in the interval $\left(0,1\right)$ .

Proof. We examine two cases: (i) Case One: When $〚{\tilde{I}}_T〛\le b_{\text{max}}$. In this case, $P\left({\tilde{H}}_L^{\text{*}}\right)$, a third degree polynomial in $H_L^{\text{*}}$ as defined in Equation (48) is a continuous function on the interval $\left[0,1\right]$ .

$\begin{array}{c}P\left(0\right)=q_0=A_7\left(A_9+\tilde{r}\right)+\frac{1}{{\tilde{\theta }}_{0,0}}\left({\tilde{r}}_D{A}_4{A}_8-\tilde{r}\right)\\ =\frac{1}{{\tilde{\theta }}_{0,0}}\left[A_7{A}_8\left({\tilde{\theta }}_{0,0}{A}_6-a_0\right)+{\tilde{r}}_D{A}_4{A}_8+\tilde{r}\left({ℛ}_{0_a}^h-1\right)\right];{ℛ}_{0_a}^h={\tilde{\theta }}_{0,0}{A}_7\\ =\frac{1}{{\tilde{\theta }}_{0,0}}\left[{ℛ}_{0_a}^h\left(\tilde{r}+A_6{A}_8\right)+{\tilde{r}}_D{A}_4{A}_8-\left(\tilde{r}+a_0{A}_7{A}_8\right)\right],\end{array}$

and $P\left(0\right)\le 0$ whenever inequality (67) holds. Also, after some calculations,

$\begin{array}{c}P\left(1\right)=q_3+q_2+q_1+q_0=\frac{1}{{\tilde{\theta }}_{0,0}}\left[{\tilde{\theta }}_{0,0}{A}_5{\widehat{\chi }}_0+{\tilde{\theta }}_{0,0}{A}_7-1\right]\\ =\frac{1}{{\tilde{\theta }}_{0,0}}\left[{\tilde{\theta }}_{0,0}\times \frac{{\widehat{\chi }}_0}{{\tilde{\theta }}_{0,0}}\left({\tilde{\theta }}_{0,0}{M}_3+{\widehat{\theta }}_0{M}_4\right)+\left({\tilde{\theta }}_{0,0}{M}_1+{\widehat{\theta }}_0{M}_2\right)-1\right]\\ =\frac{1}{{\tilde{\theta }}_{0,0}}\left[\left({ℛ}_{0_a}^h+{ℛ}_{0_a}^e\right)-1\right],\end{array}$

and $P\left(1\right)>0$ whenever $\left({ℛ}_{0_a}^h+{ℛ}_{0_a}^e\right)>1$ . Thus, by the Intermediate Value Theorem, there exists at least one positive real solution of the equation $P\left(H_L^{\text{*}}\right)=0$ , for ${\tilde{H}}_L^{\text{*}}\in \left(0,1\right)$ whenever $\left({ℛ}_{0_a}^h+{ℛ}_{0_a}^e\right)>1$ .

(ii) Case Two: When $〚{\tilde{I}}_T〛>b_{\text{max}}$. In this case, $P_3\left({\tilde{H}}_L^{\text{*}}\right)$ is as defined in (49), and proceeding in a similar manner, Equation (68) holds. Thus, by the Intermediate Value Theorem, there exists at least one positive real solution of ${\tilde{H}}_L^{\text{*}}\in \left(0,1\right)$.

□

Lemma 7. The disease free steady states defined by Theorem 5 have the following stability properties:

1) The trivial steady state, $E_0^{\text{*}}$ , is locally unstable for all parameters, and

2) The non-trivial steady state, $E_{dfe}^{\text{*}}$ , is locally asymptotically stable if ${ℛ}_0<1$ and unstable if ${ℛ}_0>1$ .

Proof. The local stability of $E_{dfe}^{\text{*}}$ is determined by the eigenvalues of the Jacobian matrix evaluated at this equilibrium. By exploiting the block structure of the Jacobian and applying the Gershgorin Circle Theorem [33], sufficient conditions ensuring that all eigenvalues have negative real parts are obtained. For ${\tilde{I}}_T\le b_{\text{max}}$, these conditions simplify to:

$\begin{array}{l}{ℛ}_c^1=\left({\tilde{\theta }}_{0,0}{\widehat{\theta }}_0{\tilde{\beta }}_{0,0}{\text{Λ}}^2{\text{Λ}}_1\right)\times \left({\tilde{\omega }}_{F0}+\tilde{\beta }+{\tilde{\sigma }}_{F0}\right)\\ \times \left({\widehat{\omega }}_{F0}+\widehat{\gamma }+{\widehat{\sigma }}_{F0}+\tilde{\eta }\right)\\ \times \left({\widehat{\delta }}_{F0}+{\widehat{\delta }}_H+{\widehat{\delta }}_T\right)<1.\end{array}$(69)

where $\Lambda =1+\widehat{\epsilon }+{\widehat{\phi }}_{0,0}+{\widehat{\zeta }}_0+{\widehat{\varphi }}_0+{\widehat{d}}_0+{\widehat{\chi }}_0$ and ${\Lambda }_1={\tilde{\psi }}_{F0}+{\tilde{\nu }}_Q+{\tilde{\nu }}_{F0}+{\tilde{\nu }}_H+{\tilde{\nu }}_T$ . Similarly, for $〚{\tilde{I}}_T〛>b_{\text{max}}$, we obtain

$\begin{array}{l}{ℛ}_c^2=\left({\tilde{\theta }}_{0,0}{\widehat{\theta }}_0{\tilde{\beta }}_{0,0}{\Lambda }^2\right)\times \left({\tilde{\xi }}_{F0}+\tilde{\gamma }+{\tilde{\rho }}_{F0}\right)\\ \times \left({\widehat{\delta }}_{F0}+{\widehat{\delta }}_H+{\widehat{\delta }}_T\right)\\ \times \left({\tilde{\psi }}_{F0}+{\tilde{\nu }}_Q+{\tilde{\nu }}_{F0}+{\tilde{\nu }}_H+{\tilde{\nu }}_T\right)<1.\end{array}$(70)

Hence, $E_{dfe}^{\text{*}}$ is locally asymptotically stable whenever the corresponding threshold condition is satisfied. □

Theorem 8 (Global Stability of the non-trivial disease free steady state). The non-trivial disease free steady state of system (40) defined in Theorem 5 is globally asymptotically stable whenever ${ℛ}_0<1$ and unstable whenever ${ℛ}_0>1$ .

Proof. Using the result of Kamgang and Sallet [38], system (40) can be written in the form:

$\{\begin{array}{l}{\dot{x}}_1=A_1\left(x\right)\cdot \left(x_1-x_1^{\text{*}}\right)+A_{12}\left(x\right)\cdot x_2\\ {\dot{x}}_2=A_2\left(x\right)\cdot x_2\end{array}$(71)

where $x={\left({\tilde{S}}_0,{\tilde{S}}_{F0},{\tilde{S}}_Q,{\tilde{I}}_{F0},{\tilde{I}}_H,{\tilde{I}}_T,{\tilde{R}}_R,{\tilde{R}}_D,{\tilde{E}}_v,\tilde{H}\right)}^{\text{T}}\in {ℝ}_{+}^{10}$ is the vector of states. We then decompose this vector of states into $x_1={\left({\tilde{S}}_0,{\tilde{R}}_R,\tilde{H}\right)}^{\text{T}}\in {ℝ}_{+}^3$ which represents densities of susceptible individuals, recovered individuals and the total human population with cadavers, and $x_2={\left({\tilde{S}}_{F0},{\tilde{S}}_Q,{\tilde{I}}_{F0},{\tilde{I}}_H,{\tilde{I}}_T,{\tilde{R}}_D,{\tilde{E}}_v\right)}^{\text{T}}\in {ℝ}_{+}^7$ which represents densities of suspected freely roaming/quarantined individuals, infected/ infectious freely roaming individuals, infected/infectious individuals at homes and treatment centres, disease induced death, and the environmental contributing

compartment respectively, so that $x=\left(x_1,x_2\right)$ . Thus, ${\dot{x}}_1={\left({\dot{\tilde{S}}}_0,{\dot{\tilde{R}}}_R,\dot{\tilde{H}}\right)}^{\text{T}}$ and

${\dot{x}}_2={\left({\dot{\tilde{S}}}_{F0},{\dot{\tilde{S}}}_Q,{\dot{\tilde{I}}}_{F0},{\dot{\tilde{I}}}_H,{\dot{\tilde{I}}}_T,{\dot{\tilde{R}}}_D,{\dot{\tilde{E}}}_v\right)}^{\text{T}}$, where $(\cdot )$ represents derivative with respect to time. When $〚{\tilde{I}}_T〛\le b_{\text{max}}$, the matrices $A_1$ , $A_{12}$ and $A_2$ are given as;

$A_1\left(x\right)=\left(\begin{array}{ccc}-1-r\left(1-2h{\tilde{H}}_L^{\text{*}}\right)& {\tilde{z}}_0+{\tilde{\alpha }}_{R0}& {\tilde{z}}_{110}\\ 0& -{\gamma }_{R_R}& 0\\ {\tilde{z}}_0& {\tilde{z}}_0& -1\end{array}\right),$(72)

$A_{12}\left(x\right)=\left(\begin{array}{ccccccc}{\tilde{z}}_{12}& {\tilde{z}}_{13}& {\tilde{z}}_{14}& {\tilde{z}}_{15}& {\tilde{z}}_{16}& {\tilde{z}}_{18}& {\tilde{z}}_{19}\\ 0& 0& {\gamma }_{R_R}{\tilde{r}}_{F0}& {\gamma }_{R_R}{\tilde{r}}_H& {\gamma }_{R_R}{\tilde{r}}_T& 0& 0\\ {\tilde{z}}_0& {\tilde{z}}_0& {\tilde{z}}_0& {\tilde{z}}_0& {\tilde{z}}_0& -{\tilde{r}}_D& 0\end{array}\right),$ (73)

and

$\begin{array}{l}{A}_2\left(x\right)=\left(\begin{array}{ccccccc}{\tilde{z}}_{22}& {\tilde{z}}_{23}& {\tilde{z}}_{24}& {\tilde{z}}_{25}& {\tilde{z}}_{26}& {\tilde{z}}_{28}& {\tilde{z}}_{29}\\ {\tilde{z}}_{32}& {\tilde{z}}_{33}& {\tilde{z}}_{34}& {\tilde{z}}_{35}& {\tilde{z}}_{36}& {\tilde{z}}_{38}& z_{39}\\ {\gamma }_{I_{F0}}{\tilde{\beta }}_{0,0}& 0& -{\gamma }_{I_{F0}}& 0& 0& 0& 0\\ z_{52}& z_{53}& z_{54}& -{\gamma }_{I_H}& 0& 0& 0\\ z_{62}& z_{63}& z_{64}& z_{65}& -{\gamma }_{I_T}& 0& 0\\ 0& 0& z_{84}& z_{85}& z_{86}& -{\gamma }_{R_D}& 0\\ z_{92}& z_{93}& z_{94}& z_{95}& z_{96}& 0& -{\gamma }_{E_v}\end{array}\right),\hfill \end{array}$ (74)

where the ${\tilde{z}}_{i,j}$ ’s and the $z_{i,j}$ ’s are matrix entries expressed in terms of the corresponding scaled parameters. Let us now consider the bounded set $\mathcal{D}$ defined by

$\mathcal{D}=\left\{\left(x_1,x_2\right)\in {ℝ}_{+}^{10}\subset {ℝ}_{+}^{3N+10},x_1\ne 0\right\}.$ (75)

We write the disease free steady state $E_{dfe}^{\text{*}}={\left(1,0,0,0,0,0,0,0,0,1\right)}^{\text{T}}$ as $x^{\text{*}}=\left(x_1^{\text{*}},0\right)$ with $x_1^{\text{*}}={\left(1,0,1\right)}^{\text{T}}\in {ℝ}_{+}^3$ and $0={\left(0,0,0,0,0,0,0\right)}^{\text{T}}\in {ℝ}_{+}^7$ .

Using these notations and following the approach in [38], one can easily verify that conditions H₁ - H₅ are fully satisfied by our system in the bounded set, $\mathcal{D}$ as defined in Equation (75). This therefore establishes the fact that the non-trivial steady state $E_{dfe}^{\text{*}}=\left(1,0,0,0,0,0,0,0,0,1\right)$ is GAS for system (40). □

4. Sensitivity Analysis and Numerical Simulations

In this section, we carry out simulations in the original parameters of the system when the treatment centres are not full, that is $〚I_T〛\le B_{\text{max}}$ . Thus system (22) with no barrier measures is given by:

$\begin{array}{l}\frac{\text{d}{S}_0}{\text{d}t}=r_h\left(1-\frac{H_L}{K_h}\right)H_L+{\alpha }_{0,0}{S}_{F0}+{\alpha }_{Q0}{S}_Q+{\alpha }_{R0}{R}_R-\left({\lambda }_{0_0}+\mu \right)S_0;\\ \frac{\text{d}{S}_{F0}}{\text{d}t}={\theta }_{0,0}{\lambda }_{0_0}{S}_0-\left({\xi }_{F0}+{\omega }_{F0}+{\beta }_{0,0}+{\alpha }_{0,0}+\mu \right)S_{F0};\\ \frac{\text{d}{S}_Q}{\text{d}t}=\left(1-{\theta }_{0,0}\right){\lambda }_{0_0}{S}_0-\left(\beta +\gamma +{\alpha }_{Q0}+\mu \right)S_Q;\\ \frac{\text{d}{I}_{F0}}{\text{d}t}={\beta }_{0,0}{S}_{F0}-\left({\sigma }_{F0}+{\rho }_{F0}+r_{F0}+{\delta }_{F0}+\mu \right)I_{F0};\\ \frac{\text{d}{I}_H}{\text{d}t}=\beta S_Q+{\xi }_{F0}{S}_{F0}+{\rho }_{F0}{I}_{F0}-\left(r_H+\eta +{\delta }_H+\mu \right)I_H;\\ \frac{\text{d}{I}_T}{\text{d}t}=\gamma S_Q+{\omega }_{F0}{S}_{F0}+{\sigma }_{F0}{I}_{F0}+\eta I_H-\left(r_T+{\delta }_T+\mu \right)I_T;\\ \frac{\text{d}{R}_R}{\text{d}t}=r_{F0}{I}_{F0}+r_T{I}_T+r_H{I}_H-\left({\alpha }_{R0}+\mu \right)R_R;\\ \frac{\text{d}{R}_D}{\text{d}t}={\delta }_{F0}{I}_{F0}+{\delta }_T{I}_T+{\delta }_H{I}_H-r_D{R}_D;\\ \frac{\text{d}{E}_v}{\text{d}t}={\psi }_{F0}{S}_{F0}+{\nu }_{F0}{I}_{F0}+{\nu }_Q{S}_Q+{\nu }_H{I}_H+{\nu }_T{I}_T-d{E}_v.\end{array}\}$ (76)

The equations for the total living human populations

$\frac{\text{d}{H}_L}{\text{d}t}=r_h\left(1-\frac{H_L}{K_h}\right)H_L-\mu H_L-{\delta }_T{I}_T-{\delta }_H{I}_H-{\delta }_{F0}{I}_{F0};$ (77)

$\frac{\text{d}H}{\text{d}t}=r_h\left(1-\frac{H_L}{K_h}\right)H_L-\mu H-\left(r_D-\mu \right)R_D;$ (78)

and

${\lambda }_{0_0}=\frac{1}{H}\left[{\tau }_{0,0}{S}_{F0}+{\phi }_{0,0}{I}_{F0}+{\epsilon }_0{S}_Q+{\varphi }_0{I}_T+{\zeta }_0{I}_H+d_0{R}_D\right]+{\chi }_0\left(\frac{B_v{E}_v}{A_v+E_v}\right).$ (79)

4.1. Sensitivity Analysis

Sensitivity analysis identifies parameters that most influence ${ℛ}_0$ and hence disease transmission. It may be local or global; here, the Forward Sensitivity Analysis is employed to quantify parameter effects through derivatives of ${ℛ}_0$ at a fixed point [39]. The corresponding parameter set is therefore considered in the original parameterization as follows:

$\begin{array}{l}\overline{\mathscr{H}}=\{{\tau }_{0,0},{\phi }_{0,0},{\epsilon }_0,{\zeta }_0,{\varphi }_0,d_0,{\chi }_0,B_v,A_v,{\sigma }_{F0},{\rho }_{F0},r_{F0},{\delta }_{F0},K_h,r_h,\mu ,\gamma ,\beta ,{\alpha }_{Q0},\\ {\xi }_{F0},{\omega }_{F0},{\alpha }_{0,0},{\beta }_{0,0},\eta ,{\psi }_{F0},{\nu }_{F0},{\nu }_Q,{\nu }_H,{\nu }_T,d,r_H,r_T,{\theta }_{0,0},r_D,{\delta }_H,{\delta }_T,{\alpha }_{R0}\}.\end{array}$

According to Chitnis et al. [40], the normalized forward sensitivity index of ${ℛ}_0$ with respect to a factor, say $x\in \overline{\mathscr{H}}$ is the ratio of the relative change in ${ℛ}_0$

to the relative change in $x$ and represented as: ${\Upsilon }_x^{{ℛ}_0}=\frac{\partial {ℛ}_0}{\partial x}.\frac{x}{{ℛ}_0}$ . For example,

using the baseline parameter values defined on Tables 2-4, we obtained specific values of the normalized local sensitivity coefficients, summarized in Table 5 and illustrated in Figure 3. With this set of values, we got ${ℛ}_0\approx 2.53$ .

Figure 3. Graph of sensitivity indices of ${ℛ}_0$ with respect to the model parameters.

A positive sensitivity index indicates that increasing (decreasing) a parameter leads to an increase (decrease) in ${ℛ}_0$ , while a negative index implies the opposite effect. For example, ${\theta }_{0,0}=+0.336277$ (Figure 3, Table 5) implies that a 10% increase (decrease) in ${\theta }_{0,0}$ results in a 3.36277% increase (decrease) in ${ℛ}_0$ , thereby enhancing (reducing) disease transmission. Biologically, ${\theta }_{0,0}$ represents the fraction of susceptible individuals who become suspected and freely roaming in the community, suggesting that reducing this rate can help control disease spread.

Table 5. Local sensitivity indices of the threshold parameter ${ℛ}_0$ using the parameters’ baseline values on Tables 2-4.

Parameter	Sen. Index	Parameter	Sen. Index
${\theta }_{0,0}$	+0.336277	$\mu$	−0.183905
${\chi }_0$	+0.201837	${\alpha }_{R0}$	0.000000
$d$	−0.201837	$K_h$	+0.201837
$\gamma$	−0.009693	$r_H$	−0.091473
${\alpha }_{0,0}$	−0.238855	${\nu }_H$	+0.034792
$\beta$	+0.002696	${\nu }_Q$	+0.009601
${\zeta }_0$	+0.091413	${\psi }_{F0}$	+0.115214
$r_D$	−0.072646	$d_0$	+0.072646
${\epsilon }_0$	+0.037839	$r_T$	−0.062352
${\tau }_{0,0}$	+0.283790	${\nu }_{F0}$	+0.023996
$r_{F0}$	−0.030942	${\nu }_T$	+0.018234
${\phi }_{0,0}$	+0.056741	${\varphi }_0$	+0.053897
${\beta }_{0,0}$	−0.012682	${\rho }_{F0}$	−0.012016
$\eta$	−0.010151	${\sigma }_{F0}$	−0.017267
${\delta }_{F0}$	−0.020474	${\delta }_H$	−0.024526
${\delta }_T$	−0.009748	${\xi }_{F0}$	−0.106333
${\omega }_{F0}$	−0.040885	$B_v$	0.000000
${\alpha }_{Q0}$	−0.040400	$\gamma$	−0.009693

It is important to note that model parameters influence ${ℛ}_0$ differently. To illustrate the effects of parameter variations, Figure 4 represents 2-D, and 3-D plots of ${ℛ}_0$ as a function of: (i) ${\phi }_{0,0}\in \left(0.0002,1.12\right)$ and ${\nu }_{F0}\in \left(0,0.5\right)$ sub-plots 4(a)-(c); (ii) ${\delta }_{F0}\in \left(0,0.2\right)$ and $d\in \left(0.0625,1.275\right)$ sub-plot 4(b)-(d) respectively with all other parameters as given on Tables 2-4 fixed.

(a) (b)

(c) (d)

Figure 4. 2-D, and 3-D plots of ${ℛ}_0$ .

4.2. Numerical Simulation

In this subsection, we perform some numerical simulations on the system (76), using the parameter values on Tables 2-4. In all choices, we used the initial conditions:

$\begin{array}{l}{S}_0\left(0\right)=970000;S_{F0}\left(0\right)=0;S_Q\left(0\right)=0;I_{F0}\left(0\right)=10;I_H\left(0\right)=0;\\ I_T\left(0\right)=0;R_R\left(0\right)=0;R_D\left(0\right)=0;E_v\left(0\right)=10.\end{array}$ (80)

Scenario I: The impact of effective contact rate. Using the baseline parameter values in Tables 2-4, we obtain ${ℛ}_0^h\approx 1.89$ , ${ℛ}_0^e\approx 1.62$ , and ${ℛ}_0\approx 2.53>1$ , with the corresponding numerical trajectories shown in Figure 5(a), Figure 5(c), Figure 5(e) and Figure 5(g). Reducing the effective contact rate to 10% (${\phi }_{0,0}=0.018$ ), while keeping all other parameters unchanged, yields ${ℛ}_0^h\approx 0.189$ , ${ℛ}_0^e\approx 0.162$ , and ${ℛ}_0\approx 0.508<1$ , as illustrated in Figure 5(b), Figure 5(d), Figure 5(f) and Figure 5(h). These results demonstrate the substantial impact of reducing effective contact on disease transmission.

(a) (b)

(c) (d)

(e) (f)

(g) (h)

Figure 5. Time series curves for system (76).

Scenario II: Estimation of cumulative infection burden and peak sensitivity.

To quantify the cumulative disease burden over the 300-day simulation period, we compute the time-integrated infectious prevalence

$A\left(\nu \right)={\displaystyle {\int }_0^{300}{I}_{F0}\left(t;\nu \right)\text{d}t},$ (81)

where $\nu$ denotes a parameter of interest. The quantity $A\left(\nu \right)$ represents the cumulative infectious person-time (person-days), thereby accounting for both the magnitude and duration of infection. Since $I_{F0}\left(t;\nu \right)$ is obtained numerically, the integral is evaluated using the composite Simpson’s rule implemented in scipy.integrate.simpson, with 2001 uniformly spaced points over $\left[0,300\right]$ days, corresponding to a step size of $\text{Δ}t=0.15$ as illustrated in Figure 6(a) and Figure 6(c).

(a) (b)

(c) (d)

Figure 6. Cumulative burden and peak infectious to parameter variations over 300 days for system (76).

To further assess the influence of key epidemiological parameters, a numerical sensitivity analysis was performed over the same 300-day period. Each parameter range was discretized into 1000 equally spaced values, and the corresponding peak values of $I_{F0}$ and $E_v$ were computed. These peaks represent, respectively, the maximum number of confirmed infectious freely roaming individuals and the maximum density of viral particles shed into the environment during the simulation period as illustrated in Figure 6(b) and Figure 6(d).

5. Discussion and Conclusions

We developed a comprehensive mathematical framework for the long-term transmission dynamics of viral diseases that incorporates human behavioural dynamics, non-pharmaceutical interventions, quarantine, re-infection, and environmental viral contamination. The model distinguishes individuals according to their adoption of barrier measures and allows behavioural responses to evolve across disease stages. Treatment-centre capacity is incorporated through two occupancy-dependent regimes. When the treatment-centre occupancy reaches the prescribed capacity $B_{\text{max}}$ , newly arising cases that cannot be accommodated are redirected to home isolation. The full-capacity and non-full-capacity regimes are therefore formulated and analysed separately, providing the appropriate framework for describing the resulting regime switching at the capacity boundary. Fundamental properties of the model, including existence, uniqueness, non-negativity, and boundedness of solutions, were established. The model was subsequently reparametrized and non-dimensionalized to facilitate its mathematical analysis and numerical investigation.

The analysis of the scaled system considered both the general disease-free model with $N$ barrier measures and the epidemiologically relevant case without barrier measures. The disease-free model admits a non-trivial globally stable equilibrium, with the scaled total human population remaining bounded. For the epidemiological model, the existence of disease-free and endemic equilibria was established, and the basic reproduction number ${ℛ}_0$ was derived using two transmission scenarios corresponding to treatment centres being non-full or full. In both cases, ${ℛ}_0$ naturally decomposes into human-to-human and environmental contributions, allowing the relative roles of quarantined and non-quarantined individuals and the different transmission pathways to be assessed.

Numerical simulations using SARS-CoV-2-related parameter values demonstrate the strong influence of effective contact rates on disease transmission. In particular, reducing the effective contact rate between susceptible and freely roaming infectious individuals substantially decreases ${ℛ}_0$ and can drive the system towards disease control when ${ℛ}_0<1$ . This highlights the importance of effective contact tracing, rapid identification of infectious individuals, and appropriate quarantine enforcement. However, when contact rates remain sufficiently high, quarantine alone may be inadequate to eliminate transmission and may instead support persistent endemic infection. The model further demonstrates that environmental contamination can make a non-negligible contribution to viral transmission, particularly during the early stages of an outbreak. Incorporating environmental viral particles therefore provides a more comprehensive representation of transmission pathways and enhances the potential applicability of the framework to broader One Health settings. The results also emphasize that accurate estimation of ${ℛ}_0$ is essential for effective policy design, since both overestimation and underestimation may result in inappropriate levels of intervention and associated epidemiological or socio-economic consequences.

Overall, the findings demonstrate that human behavioural dynamics, contact tracing, isolation, quarantine, and environmental exposure play important roles in shaping viral disease transmission. Although the basic epidemiological model incorporates intervention-related mechanisms such as contact tracing, isolation, and quarantine, the present epidemiological threshold and numerical analyses primarily consider the no-barrier case. Consequently, the implications regarding adherence to multiple barrier measures should be interpreted prospectively. The adherence-transition rates are treated as fixed parameters within each simulation and represent behavioural changes that may reflect prevailing disease risk, public awareness, beliefs, and policy sensitization. Thus, the framework captures behavioural responses through adherence stratification without imposing an explicit prevalence-dependent feedback mechanism, while providing a flexible basis for investigating how adherence, reductions in effective contact rates, and reduced environmental exposure may influence disease transmission.

A limitation of the present study is that the numerical simulations are scenario-based rather than calibrated predictions, as several transition, contact, and viral-shedding parameters are assumed or weighted using available information rather than estimated from a specific population-level dataset. Explicit analysis of interactions among multiple barrier measures, vaccination and other pharmaceutical interventions, as well as stochastic effects, remains an important direction for future work. In addition, incorporating well-defined epidemiological datasets and appropriate statistical parameter-estimation methods will improve model calibration and predictive capacity, thereby providing a more comprehensive representation of viral disease dynamics.

Acknowledgements

JAN, GAN and MIT-E acknowledge the component of Canada’s International Development Research Center Grant No. 109559-001, to the University of Buea in 2020, when aspects of this work were conceived. GAN and JAN also acknowledge support from the Cameroonian Ministry of Higher Education through the initiative for the modernisation of research in Cameroon’s Higher Education for the 2023, 2024 and 2025 granting seasons. JAN acknowledges support from the Erasmus+ student mobility grants enabling him to make significant progress on this work at Cotbus university in Germany. He additionally acknowledges support from the CDC-IMU programme for Graduate Assistantship in Developing Countries (GRAID) for the 2021, 2022 and 2023 granting seasons. All Authors acknowledge the sponsorship of the Commission for Developing Countries (CDC) in conjunction with the International Mathematics Union (IMU) through the CDC-ADMP (African Diaspora Mathematicians Program) grant that made it possible for interactive collaborative work during the 2018 and 2019 grant sponsored visits to the University of Buea by MIT-E during which period the international collaboration was forged.

Author Contributions

Conceptualization, the model was jointly developed by all three authors; methodology, all three authors agreed on the methodology and the analyses presented; formal analysis, was conducted by James A. Njongwe under the surpervisor of G.A. Ngwa and M.I. Teboh-Ewungkem; writing, review and editing, was done by all three authors in a concensus framework; visualization, James A. Njongwe did the simulations; supervision, was done by G.A. Ngwa and M.I. Teboh-Ewungkem; project administration, G.A. Ngwa; funding acquisition, G.A. Ngwa and M.I. Teboh-Ewungkem. All authors have read and agreed to the published version of the manuscript.