Positivity and boundedness
Since the system (5) represents the population in each class and all model parameters are positive, it lies in a region \(\Omega\) defined by
$$\begin{aligned} \Omega&= \Omega _{_{H}}\times \Omega _{_{P}} \in \mathbb {R}^{7}_{_{+}}\times \mathbb {R}^{2}_{_{+}} = \mathbb {R}^{9}_{_{+}} , ~~\text {where}\\ \Omega _{_{H}}&= \left\{ \begin{bmatrix} S(t)\\ E(t)\\ I_{_{T}}(t)\\ I_{_{C}}(t)\\ I_{_{CT}}(t)\\ T(t)\\ R(t) \end{bmatrix} \in \mathbb {R}^{7}_{_{+}}\left| \begin{matrix} S(t) > 0,\\ E(t)\ge 0,\\ I_{_{T}}(t)\ge 0,\\ I_{_{C}}(t)\ge 0,\\ I_{_{CT}}(t)\ge 0,\\ T(t)\ge 0,\\ R(t)\ge 0. \end{matrix}\right. \right\} ,~~ \hspace{5mm} \text {and} \hspace{5mm} \Omega _{_{P}} = \left\{ \begin{bmatrix} B_T(t)\\ B_C(t) \end{bmatrix} \in \mathbb {R}^{2}_{_{+}}\left| \begin{matrix} B_T(t) \ge 0,\\ B_C(t) \ge 0. \end{matrix}\right. \right\} . \end{aligned}$$
For the model (5) to be mathematically and biologically meaningful, it is necessary to prove that its solutions are all positive if their initial values are non-negative. In this subsection, we investigate the positivity of the solutions of the model and establish their bounds.
Positivity of the solution of the model
Theorem 1
The solutions of the model (5) are positive in the region \(\Omega\) if its initial values are non-negative for all \(t\ge 0\).
Consider the ordinary differential equation for susceptible human \(S\) of the model (5),
$$\begin{aligned} \dfrac{dS}{dt} = \Lambda +\alpha R – (\lambda +\mu _{1})S\ge -(\lambda +\mu _{1})S, \end{aligned}$$
and observe that
$$\begin{aligned} \dfrac{dS}{dt} \ge -(\lambda +\mu _{1})S. \end{aligned}$$
(6)
Upon separating variables and integrating (6), we obtain
$$\begin{aligned} S(t) \ge S_(0)e^{-\int _{_{_{0}}}^{^{t}}{(\lambda (t)+\mu _{1})dt}}\ge 0. \end{aligned}$$
By the same process, we obtain
$$\begin{aligned} E(t)&\ge E(0)e^{-(\tau + \mu _{1})t}\ge 0,~~I_{T}(t) \ge I_{T}(0) e^{- (r_{1}+k_{1}+\Psi _{1}+\mu _{1})t}\ge 0,~~\\&I_{C}(t) \ge I_{C}(0)e^{-\int _{_{_{0}}}^{^{t}}{(\mu _{1}+r_{3}+k_{3}+\Psi {3})t}} \ge 0, ~~I_{CT}(t) \ge I_{CT}(0) e^{- (r_{2}+\sigma +k_{2}+\mu _{1})t}\ge 0,\\&~~T(t) \ge T(0) e^{- (\Psi _{2}+\mu _{1})t}\ge 0,~~R(t) \ge R(0)e^{- \mu _{1}t}\ge 0, ~~B_{T}(t) \ge B_{T}(0)e^{- \mu _{3}t}\ge 0~~\\ \text {and} \hspace{3mm} B_{C}(t)&\ge B_{C}(0)e^{- \mu _{2}t}\ge 0. \end{aligned}$$
This implies that all solutions of the model (5) are positive for all \(t\ge 0\). This ends the proof of Theorem 1.
Boundedness of the model’s solutions
The boundedness of the solutions of the model (5) in the region \(\Omega\) is demonstrated by considering the total rate of change of the human population, \(\dfrac{dN}{dt}\), governed by
$$\begin{aligned} \dfrac{dN}{dt} = \dfrac{dS}{dt} + \dfrac{dE}{dt} + \dfrac{dI_{T}}{dt}+\dfrac{dI_{C}}{dt}+ \dfrac{dI_{CT}}{dt}+ \dfrac{dT}{dt}+\dfrac{dR}{dt} = \Lambda -\mu _{_{1}}N, \end{aligned}$$
and observe that
$$\begin{aligned} \dfrac{dN}{dt} \le \Lambda -\mu _{_{1}}N. \end{aligned}$$
(7)
Using Birkhoff and Rota’s theorem17 on differential inequalities, along with the method of separation of variables applied to the inequality (7), and after substituting the initial conditions, we obtain
$$\begin{aligned} N(t) \le \dfrac{\Lambda }{\mu _{_{1}}} + (N(0)-\dfrac{\Lambda }{\mu _{_{1}}})e^{-\mu _{_{1}}t}. \end{aligned}$$
(8)
As \(t\rightarrow \infty\), inequality (8) reduces to
$$\begin{aligned} N(t) \le \dfrac{\Lambda }{\mu _{_{1}}}, \end{aligned}$$
Since human population is non-negative, implies that
$$\begin{aligned} 0 \le N(t) \le \dfrac{\Lambda }{\mu _{_{1}}}. \end{aligned}$$
Similarly, it follows that the upper bound for the bacterial concentration of salmonella Typhi is bounded above by
$$\begin{aligned} \frac{dB_{_{T}}}{dt}= g_{{T}}B_{{T}} \left( 1-{\frac{B_{{T}}}{K_{{T}}}} \right) +g_{{1}}I_{{T}}+g_{{2}}I_{{{ CT}}}-\mu _{{3} }B_{{T}} \le g_{_{T}} B_{_{T}} \left( 1 – \frac{B_{_{T}}}{K_{_{T}}}\right) + (g_{_{1}} + g_{_{2}}) \frac{\Lambda }{\mu _{_{1}}} – \mu _{_{3}} B_{_{T}} \end{aligned}$$
(9)
From inequality (Eq. 7),
If \(B_{_{T}} \ge \frac{(g_{_{1}} + g_{_{2}}) \Lambda }{\mu _{_{1}}}\), then
$$\begin{aligned} \frac{dB_{_{T}}}{dt} \le (g_{_{T}} – \mu _3) B_{_{T}} – \frac{g_{_{T}}}{K_{_{T}}} (B_{_{T}})^2 + B_{_{T}} \end{aligned}$$
(10)
After differentiate (Eq. 10)
$$\begin{aligned} (g_{_{T}} – \mu _{_{3}}+1) B_{_{T}} \left( 1 – \frac{g_{_{T}} B_{_{T}}}{K_{_{T}} (g_{_{T}} – \mu _{_{3}} + 1)}\right) \end{aligned}$$
(11)
The constant
$$\begin{aligned} \frac{K_{_{T}} (g_{_{T}} – \mu _{_{3}} + 1)}{g_{_{T}}} \end{aligned}$$
(12)
is the upper bound for the differential inequality (11) since (11) is the logistic growth model with carrying capacity (12). For some \(t \ge 0\), \(\frac{(g_{1} + g_{2}) \Lambda }{\mu _{1}}\) is an upper bound for \(B_{T}\) where (10) is false, while \(B_{T}\) is bounded above by (12) for the rest of the time points in the domain \(B_{T}\) if ( 10) is true. The constant \(\frac{(g_{1} + g_{2}) \Lambda }{\mu _{1}}\) is the maximum shedding rate from the Typhoid fever infected individuals and dually infected individuals. In both cases
$$\begin{aligned} B_{_{T}} \le \max \left\{ \frac{K_{_{T}} (g_{_{T}} – \mu _{_{3}} + 1)}{g_{_{T}}}, \frac{(g_{_{1}} + g_{_{2}}) \Lambda }{\mu _{_{1}}}\right\} \end{aligned}$$
(13)
Within the feasible region
$$\begin{aligned} \Omega&= \{S, E, I_{_{T}}, I_{_{C}}, I_{_{CT}}, T, R, B_{_{T}}, B_{_{C}}\} \in \mathbb {R}_+^9 \, | \, 0 \le N \le \frac{\Lambda }{\mu _{_{1}}}, \nonumber \\ B_{_{T}}&\in \left[ 0, \max \left\{ \frac{K_{_{T}} (g_{_{T}} – \mu _{_{3}} + 1)}{g_{_{T}}}, \frac{(g_{_{1}} + g_{_{2}}) \Lambda }{\mu _{_{1}}}\right\} \right] , \nonumber \\ B_{_{C}}&\in \left[ 0, \max \left\{ \frac{K_{_{C}} (g_{_{C}} – \mu _{_{2}} + 1)}{g_{_{C}}}, \frac{(g_{_{3}} + g_{_{4}}) \Lambda }{\mu _{_{1}}}\right\} \right] \end{aligned}$$
(14)
Hence, we can conclude that the solutions of the developed model are bounded in the positive region \(\Omega \in \mathbb {R}^{9}_{+}\). In this case, the developed model is mathematically and biologically meaningful, and it is considered for further analysis.
Equilibrium points and the effective reproduction number \(\mathcal {R}\)
Disease free equilibrium point (DFE)
A disease-free equilibrium point refers to a state in which a community, population, or system is entirely free from COVID-19 and typhoid fever. To analyze the disease-free equilibrium point, it is necessary to examine both the conditions required for its formation and the factors that contribute to the absence of disease transmission.
The DFE is established by setting the right-hand side of the model system (5) to zero and substituting zero values for the infectious states and the force of infection (i.e., \(E = I_{T} = I_{C} = I_{CT} = T = R = B_{T} = B_{C} = \lambda = 0\)).
Upon solving and simplifying the resulting system, we get
$$\begin{aligned} \text {DFE}=(S,E,I_{T},I_{_{C}},I_{_{CT}},T,R,B_{_{T}},B_{_{C}})=\left( \dfrac{\Lambda }{\mu _{_{1}}},0,0,0,0,0,0,0,0\right) ,~\text {respectively.} \end{aligned}$$
The formulated DFE is used to compute a state where there are no causative agents of COVID-19, Typhoid fever, or both. Furthermore, the DFE is used to compute the effective reproduction number \(\mathcal {R}\)of the co-infection of COVID-19 and Typhoid fever, which describes the persistence or extinction of the disease(s) within the population.
The effective reproductive number \(\mathcal {R}\)
The effective reproductive number, denoted as \(\mathcal {R}\), is a crucial metric in the study of epidemiological models for disease management. By estimating \(\mathcal {R}\), epidemiologists can assess the potential effects of an outbreak and predict its progression. Specifically, if \(\mathcal {R} , it indicates that each infected individual, on average, causes fewer than one secondary infection, leading to the eventual extinction of the disease within the population. When \(\mathcal {R} = 1\), the average number of secondary infections remains stable, it means that the disease will neither spread nor die out. Conversely, if \(\mathcal {R} > 1\), it suggests that a single infected individual generates more than one secondary infection resulting in the disease potentially spreading throughout the entire population.
Adopting the next-generation matrix approach by18,19 to formulate the formula for \(\mathcal {R}\) of the model (5) is given by
$$\begin{aligned} \mathcal {R} = l_{11} = R_{hh} + R_{ch} + R_{hc} + R_{hb} + R_{hp} + R_{cb} + R_{cc}, \end{aligned}$$
where
$$R_{_{hh}} = {\dfrac{\alpha _{_{1}}\Lambda _{_{1}}}{\mu _{_{1}} \left( r_{_{1}}+m+\mu _{_{1}} \right) }}+{\dfrac{\alpha _{_{2}}\Lambda _{_{1}}\beta _{_{1}}}{\mu _{_{1}} \left( r_{_{1}}+m+\mu _{_{1}} \right) \left( \mu _{_{3}}-g\right) }},~R_{_{ch}} = {\dfrac{\kappa _{_{2}}\Lambda _{_{2}}\beta _{_{1}}}{\mu _{_{2}} \left( r_{_{1}}+m+\mu _{_{1}} \right) \left( \mu _{_{3}}-g \right) }},$$
$$R_{_{hc}} = {\dfrac{\alpha _{_{2}}\Lambda _{_{1}}\beta _{_{2}}}{\mu _{_{1}} \left( r_{_{2}}+\phi +\mu _{_{2}} \right) \left( \mu _{_{3}}-g\right) }}+{\dfrac{\alpha _{_{3}}\Lambda _{_{1}}\theta }{\mu _{_{1}} \left( r_{_{2}}+\phi +\mu _{_{2}} \right) \mu _{_{4}}}},~R_{_{hb}} ={\dfrac{\alpha _{_{2}}\Lambda _{_{1}}}{\mu _{_{1}} \left( \mu _{_{3}}-g\right) }},~R_{_{hp}}= {\dfrac{\alpha _{_{3}}\Lambda _{_{1}}}{\mu _{_{1}}\mu _{_{4}}}},$$
$$R_{_{cb}} = {\dfrac{\kappa _{{2}}\Lambda _{_{2}}}{\mu _{_{2}} \left( \mu _{_{3}}-g\right) }}~ \text {and}~R_{_{cc}} = {\dfrac{\kappa _{_{1}}\Lambda _{_{2}}}{\mu _{_{2}} \left( r_{_{2}}+\phi +\mu _{_{2}} \right) }}+{\dfrac{\kappa _{_{2}}\Lambda _{_{2}}\beta _{_{2}}}{\mu _{_{2}} \left( r_{_{2}}+\phi +\mu _{_{2}} \right) \left( \mu _{_{3}}-g \right) }}.$$
By inspection, it is observed that \(\mathcal {R}\) depends on sub-effective reproductive numbers: \(R_{hh} + R_{ch} + R_{hc} + R_{hb} + R_{hp} + R_{cb}\) and \(R_{cc}\) These respectively define the expected number of newly infected individuals produced by one infected person with Typhoid fever, COVID-19, and co-infection; the secondary number of infected Typhoid fever patients generated by ingesting Salmonella Typhi shed by one infected Typhoid fever individual; and the secondary number of infected COVID-19 patients generated by ingesting SARS-CoV-2 shed by one infected COVID-19 individual.
Global stability of disease free equilibrium (DFE)
The disease-free equilibrium is considered stable if it can maintain its state despite perturbations or initial conditions. This concept, known as global stability, ensures that a system remains at the disease-free equilibrium regardless of initial circumstances. Global stability analysis aims to assess whether a population can stay disease-free without external influences, such as infections or changes in dynamics. Additionally, it offers insights into the effectiveness of control measures in preventing disease recurrence and replication.
Theorem 2
DFE is globally asymptotically stable whenever \(\mathcal {R} as \(t \rightarrow \infty\) and unstable otherwise.
Let the infectious state of the system (5) be denoted by
$$\begin{aligned} P=\left( E(t), I_{T}(t), I_{C}(t), I_{CT}(t), T(t), B_{T}(t), B_{C}(t) \right) ^{\text {T}}. \end{aligned}$$
Adopting the comparison principle by18 we get:
$$\begin{aligned} \dot{P} \le (F-V)P, \end{aligned}$$
(15)
where matrices \(F\) and \(V\) are evaluated at disease-free equilibrium in \(\dot{P}=(E^{‘},I^{‘}_{T},I^{‘}_{C},I^{‘}_{CT},T^{‘},B^{‘}_{T},\) \(B^{‘}_{C})^{\text {T}}\) is a vector containing the rate of change of infected states, where
$$\begin{aligned} \left. \begin{array}{lll} E^{‘} & =\lambda \,S- \left( \mu _{{1}}+\tau \right) E,\\ I^{‘}_{T} & =\tau \,nE+\sigma \,h_{{1}}I_{{{ CT}}}- \left( k_{{1}}+\Psi _{{1}}+\mu _{{1}}+r_{{1}} \right) I_{{T}},\\ I^{‘}_{C} & =\tau \,mE+\sigma \,h_{{2}}I_{{{ CT}}}- \left( \mu _{{1}}+r_{{3}}+k_{{3}}+\Psi _{{3}} \right) I_{{C}},\\ I^{‘}_{CT} & = \left( 1-(m+n) \right) \tau \,E- \left( \sigma \,h_{{1}}+\sigma \,h_{{2}}+k_{{2}} \left( 1-(h_{{1}}+h_{{2}})+\mu _{{1}}+r_{{2}} \right) \right) I_{{{ CT}}},\\ T^{‘} & =k_{{2}} \left( 1-(h_{{1}}+h_{{2}})\sigma \right) I_{{{ CT}}}+k_{{1}}I_{{T}}+k_{{3}}I_{{C}}- \left( \Psi _{{2}}+\mu _{{1}} \right) T,\\ B^{‘}_{T} & =g_{{T}}B_{{T}} \left( 1-{\frac{B_{{T}}}{K_{{T}}}} \right) +g_{{1}}I_{{T}}+g_{{2}}I_{{{ CT}}}-\mu _{{3} }B_{{T}},\\ B^{‘}_{C} & =g_{{C}}B_{{C}} \left( 1-{\dfrac{B_{{C}}}{K_{{C}}}} \right) +g_{{3}}I_{{{ CT}}}+g_{{4}}I_{{C}}-\mu _{{2}}B_{{C}}, \end{array}\right\} \end{aligned}$$
(16)
Applying the Perron-Frobenius theorem19 for Metzler matrices \(F\) and \(V^{-1}\), we observe that \(V^{-1}F\) has a dominant eigenvalue \(\mathcal {R} = \rho (V^{-1}F) = \rho (FV^{-1})\) that corresponds to a non-negative vector \(\textbf{a}\). It follows that
$$\begin{aligned} \textbf{a}^{\text {T}}V^{-1}F = \mathcal {R}\textbf{a}^{\text {T}}. \end{aligned}$$
Motivated by20, we use the Lyapunov function:
$$\begin{aligned} f(t) = \textbf{a}^{\text {T}}V^{-1}P. \end{aligned}$$
(17)
Differentiating the function (17) with respect to \(t\) along with the infected states of the system (5), and using Eq. (15), we get:
$$\begin{aligned} \dot{f} = \textbf{a}^{\text {T}}V^{-1}\dot{P} \le \textbf{a}^{\text {T}}V^{-1} (F-V)P = (\mathcal {R}-1)\textbf{a}^{\text {T}}P, \quad \text {where} \quad \textbf{a}^{\text {T}}P \ge 0. \end{aligned}$$
The stability criteria \(\dot{f} = (\mathcal {R}-1)\textbf{a}^{\text {T}}P \le 0\) is attained only for \(\mathcal {R} .
If \(\dot{f} = 0\) for \(\mathcal {R} , then \(\textbf{a}^{\text {T}}P = 0\). Because the Perron-Frobenius eigenvector \(\textbf{a}\) has non-negative elements, then
$$\begin{aligned} P = (E, I_{T}, I_{C}, I_{CT}, T, B_{T}, B_{C}) = (0, 0, 0, 0, 0, 0, 0). \end{aligned}$$
(18)
Substituting the coordinates of (18) into the model (5), we get the DFE \(E_{0} = \left( \frac{\Lambda }{\mu _{1}}, 0, 0, 0, 0, 0, 0, 0, 0 \right)\).
When \(\dot{f} , it implies that the infected variables \(I_{T}, I_{C}, I_{CT}, T, B_{T}\), and \(B_{C}\) lose their energy to transmit the disease to others, and therefore the system turns to DFE as \(t \rightarrow \infty\). The two scenarios suggest that every solution of the system (5) converges to DFE as \(t \rightarrow \infty\) for \(\mathcal {R} . By LaSalle’s invariant principle21, DFE is globally asymptotically stable whenever \(\mathcal {R} .
Global stability of endemic equilibrium \(E^{*}\)
The endemic equilibrium \(E^{*}\) is a state where Salmonella typhi and SARS-Cov-2 infections persist among individuals. During the endemic state, the infected individuals and pathogens transmit the diseases to susceptible individuals, and therefore, the number of infected individuals remains relatively constant over time. That is, the forces of infections \(\lambda _{1}\), \(\lambda _{2}\), \(\lambda _{3}\), and \(\lambda _{4} \ne 0\). The endemic equilibrium is said to be globally stable if, in the short or long run, the disease becomes endemic (does not die out). During this state, any invariant set in \(\Omega\) close to \(E^{*}\) remains close to it and eventually converges to \(E^{*}\) over an indefinite time.
Definition 1
A function \(V_{x}\) is called a Lyapunov function if it is continuous and differentiable for all values of \(x \in U\) satisfying the conditions:
-
i)
\(V(x) = 0\)
-
ii)
\(V(x) > 0 \forall x \in U \setminus \{0\}\)
-
iii)
\(V'(x) \le 0 \forall x \in U \setminus \{0\}\)
Theorem 3
The endemic equilibrium \(E^{*}\) of the model system (5) is globally asymptotically stable if the effective reproduction number \(\mathcal {R} > 1\). To prove theorem (3.3), we adopt the approach in Osman et al. (2020) by considering the non-linear function.
$$\begin{aligned} \begin{aligned} X&= S^{*}\left[ \dfrac{S}{S^{*}}-\ln {\left( \dfrac{S}{S^{*}}\right) }\right] +E^{*}\left[ \dfrac{E}{E^{*}}-\ln {\left( \dfrac{E}{E^{*}}\right) }\right] + I_T^{*}\left[ \dfrac{I_T}{I_T^{*}}-\ln {\left( \dfrac{I_T}{I_T^{*}}\right) }\right] \\&\quad + I_C^{*}\left[ \dfrac{I_C}{I_C^{*}}-\ln {\left( \dfrac{I_C}{I_C^{*}}\right) }\right] + I_{CT}^{*}\left[ \dfrac{I_{CT}}{I_{CT}^{*}}-\ln {\left( \dfrac{I_{CT}}{I_{CT}^{*}}\right) }\right] + T^{*}\left[ \dfrac{T}{T^{*}}-\ln {\left( \dfrac{T}{T^{*}}\right) }\right] \\&\quad + R^{*}\left[ \dfrac{R}{R^{*}}-\ln {\left( \dfrac{R}{R^{*}}\right) }\right] + B_T^{*}\left[ \dfrac{B_T}{B_T^{*}}-\ln {\left( \dfrac{B_T}{B_T^{*}}\right) }\right] + B_C^{*}\left[ \dfrac{B_C}{B_C^{*}}-\ln {\left( \dfrac{B_C}{B_C^{*}}\right) }\right] \end{aligned} \end{aligned}$$
(19)
It can be observed that the function \(X\) in (i) in the definition above is valid for all \(x \ge 0\) where \(x = \{ S, E, I_T, I_C, T, R, B_T, B_C \}\) and \(x^{*}\) is the endemic equilibrium. Hence, \(X\) is a Lyapunov function.
To prove if condition (iii) holds, we differentiate \(X\) along the solution of the model (5) to obtain:
$$\begin{aligned} \begin{array}{lll} \dfrac{dX}{dt} & = \left( 1-\dfrac{S^{*}}{S}\right) \dfrac{dS}{dt}+\left( 1-\dfrac{E^{*}}{E}\right) \dfrac{dE}{dt}+ \left( 1-\dfrac{I_T^{*}}{I_T}\right) \dfrac{dI_T}{dt} \\ & \quad + \left( 1-\dfrac{I_C^{*}}{I_C}\right) \dfrac{dI_C}{dt} + \left( 1-\dfrac{I_{CT}^{*}}{I_{CT}}\right) \dfrac{dI_{CT}}{dt} + \left( 1-\dfrac{T^{*}}{T}\right) \dfrac{dT}{dt} \\ & \quad + \left( 1-\dfrac{R^{*}}{R}\right) \dfrac{dR}{dt} + \left( 1-\dfrac{B_T^{*}}{B_T}\right) \dfrac{dB_T}{dt} \quad + \left( 1-\dfrac{B_C^{*}}{B_C}\right) \dfrac{dB_C}{dt}. \end{array} \end{aligned}$$
(20)
Substituting system (5) into Eq. (21) yields
$$\begin{aligned} \begin{aligned} \frac{dX}{dt}&= \left( 1-\frac{S^{*}}{S}\right) (\alpha \,R-\lambda \,S-\mu _{1}S+\Lambda ) \\&\quad +\left( 1-\frac{E^{*}}{E}\right) (\lambda \,S- (\mu _{1}+\tau ) E) \\&\quad +\left( 1-\frac{I_T^{*}}{I_T}\right) (\tau \,nE+\sigma \,h_1I_{{CT}}- (k_1+\Psi _1+\mu _{1}+r_1) I_{T}) \\&\quad +\left( 1-\frac{I_C^{*}}{I_C}\right) (\tau \,mE+\sigma \,h_2I_{{CT}}- (\mu _{1}+r_3+k_3+\Psi _3) I_{C}) \\&\quad +\left( 1-\frac{I_{CT}^{*}}{I_{CT}}\right) (1-(m+n))\tau \,E- (\sigma \,h_1+\sigma \,h_2+k_2(1-(h_1+h_2)+\mu _{1}+r_2)) I_{{CT}}) \\&\quad +\left( 1-\frac{T^{*}}{T}\right) (k_2(1-(h_1+h_2)\sigma ) I_{{CT}}+k_1I_{T}+k_3I_{C}- (\Psi _2+\mu _{1}) T) \\&\quad +\left( 1-\frac{R^{*}}{R}\right) \Psi _2T+\Psi _3I_{C} +\Psi _1I_{T})-({\alpha +\mu _{{1}}})R \\&\quad +\left( 1-\frac{B_T^{*}}{B_T}\right) (g_T B_T (1-{\frac{B_T}{K_T}}) +g_1I_{T}+g_2I_{{CT}}-\mu _3B_T) \\&\quad +\left( 1-\frac{B_C^{*}}{B_C}\right) (g_C B_C (1-{\frac{B_C}{K_C}}) +g_3I_{{CT}}+g_4I_{C}-\mu _2B_C). \end{aligned} \end{aligned}$$
(21)
Now substitute
\(S = S – S^{*}, E – E^{*}, I_T – I_T^{*}, I_C – I_C^{*}, I_{CT} – I_{CT}^{*}, T – T^{*}, R – R^{*}, B_T – B_T^{*}, B_C – B_C^{*}\) we get
$$\begin{aligned} \begin{array}{lll} \dfrac{dX}{dt} & =\left( \dfrac{S – S^{*}}{S}\right) ^2(-\lambda – \mu _1) + \Lambda \left( \dfrac{(S – S^{*})}{S}\right) + \alpha \left( \dfrac{(R – R^{*})(S – S^{*})}{S}\right) \\ & \quad + \lambda \,(S – S^{*})- \left( \mu _{{1}}+\tau \right) (E – E^{*}) + \tau \,n(E – E^{*})+\sigma \,h_{{1}}(I_{{{ CT}}} – I^{*}_{{{ CT}}})\\ & \quad – \left( k_{{1}}+\Psi _{{1}}+\mu _{{1}}+r_{{1}} \right) (I_{{T}} – I^{*}_{{T}})\\ & \quad +\tau \,m(E – E^{*})+\sigma \,h_{{2}}(I_{{{ CT}}} – I^{*}_{{{ CT}}})- \left( \mu _{{1}}+r_{{3}}+k_{{3}}+\Psi _{{3}} \right) (I_{{C}} – I^{*}_{{C}})\\ & \quad + \left( 1-(m+n)) \right) \tau \,(E – E^{*})- \left( \sigma \,h_{{1}}+\sigma \,h_{{2}}+k_{{2}} \left( 1-(h_{{1}}+h_{{2}})+\mu _{{1}}+r_{{2}} \right) \right) (I_{{{ CT}}} – I^{*}_{{{ CT}}})\\ & \quad + k_{{2}} \left( 1-(h_{{1}}+h_{{2}})\sigma \right) (I_{{{ CT}}} – I^{*}_{{{ CT}}}) + k_{{1}}(I_T – I^{*}_T)+k_{{3}}(I_{{C}} – I^{*}_{{C}}) – \left( \Psi _{{2}}+\mu _{{1}} \right) (T – T^{*})\\ & \quad -(\alpha +\mu _{{1}})(R – R^{*})+\Psi _{{2}}(T – T^{*}) + \Psi _{{3}}(I_{{C}} – I^{*}_{{C}}) + \Psi _{{1}}(I_T – I^{*}_T) \\ & \quad +g_{{T}}(B_T – B^{*}_T) \left( 1-{\frac{(B_T – B^{*}_T)}{K_{{T}}}} \right) +g_{{1}}(I_T – I^{*}_T) + g_{{2}}(I_{{{ CT}}} – I^{*}_{{{ CT}}})-\mu _{{3}}(B_T – B^{*}_T)\\ & \quad +g_{{C}}(B_C – B^{*}_C) \left( 1-{\frac{(B_C – B^{*}_C)}{K_{{T}}}} \right) +g_{{3}}(I_{CT} – I^{*}_{CT}) + g_{{4}}(I_C – I^{*}_C)-\mu _{{2}}(B_C – B^{*}_C). \end{array} \end{aligned}$$
(22)
Now, collecting positive and negative terms together in the system (22):
$$\begin{aligned} \begin{array}{lll} \dfrac{dX}{dt}&= Q – Z \end{array} \end{aligned}$$
(23)
Let,
$$\begin{aligned} \begin{array}{lll} Q & = \Lambda \left( \dfrac{(S – S^{*})}{S}\right) + \alpha \left( \dfrac{(R – R^{*})(S – S^{*})}{S}\right) \\ & \quad + \lambda \,(S – S^{*}) + \tau \,n(E-E^{*})+\sigma \,h_{{1}}(I_{{{ CT}}} – I^{*}_{{{ CT}}})\\ & \quad +\tau \,m(E – E^{*})+\sigma \,h_{{2}}(I_{{{ CT}}} – I^{*}_{{{ CT}}})\\ & \quad + \left( 1-(m+n)) \right) \tau \,(E – E^{*})\\ & \quad + k_{{2}} \left( 1-(h_{{1}}+h_{{2}})\sigma \right) (I_{{{ CT}}} – I^{*}_{{{ CT}}}) + k_{{1}}(I_T – I^{*}_T)+k_{{3}}(I_{{C}} – I^{*}_{{C}})\\ & \quad + \Psi _{{3}}(I_{{C}} – I^{*}_{{C}}) + \Psi _{{1}}(I_T – I^{*}_T) \\ & \quad +g_{{T}}(B_T – B^{*}_T) +g_{{1}}(I_T – I^{*}_T) + g_{{2}}(I_{{{ CT}}} – I^{*}_{{{ CT}}})\\ & \quad +g_{{C}}(B_C – B^{*}_C) + g_{{4}}(I_C – I^{*}_C) + g_{{3}}(I_{CT} – I^{*}_{CT}). \end{array} \end{aligned}$$
(24)
and
$$\begin{aligned} \begin{array}{lll} Z & = \left( \dfrac{S – S^{*}}{S}\right) ^2(\lambda + \mu _1)\\ & \quad + \left( \mu _{{1}} + \tau \right) (E – E^{*}) + \left( k_{{1}}+\Psi _{{1}}+\mu _{{1}}+r_{{1}} \right) (I_{{T}} – I^{*}_{{T}})\\ & \quad + \left( \mu _{{1}}+r_{{3}}+k_{{3}}+\Psi _{{3}} \right) (I_{{C}} – I^{*}_{{C}})\\ & \quad + \left( \sigma \,h_{{1}}+\sigma \,h_{{2}}+k_{{2}} \left( 1-(h_{{1}}+h_{{2}})+\mu _{{1}}+r_{{2}} \right) \right) (I_{{{ CT}}} – I^{*}_{{{ CT}}})\\ & \quad + \left( \Psi _{{2}}+\mu _{{1}} \right) (T – T^{*}) (\alpha +\mu _{{1}})(R – R^{*})\\ & \quad +g_{{T}}\frac{\left( B_T – B^{*}_T\right) ^2}{K_T} +g_{{C}}\frac{\left( B_C – B^{*}_C\right) ^2}{K_C} + \mu _{{3}}(B_T – B^{*}_T) +\mu _{{2}}(B_C – B^{*}_C).\\ \end{array} \end{aligned}$$
(25)
Then,
$$\begin{aligned} \begin{array}{lll} \dfrac{dX}{dt}&
(26)
\(S = S^{*}, E = E^{*}, I_T = I_T^{*}, I_C = I_C^{*}, I_{CT} = I_{CT}^{*}, T = T^{*}, R = R^{*}, B_T = B_T^{*},\) and \(B_C = B_C^{*}\)
Therefore, the maximum compact invariant set in \(\{ S, E, I_T, I_C, I_{CT}, T, R, B_T, B_C \in \Omega :\frac{dX}{dt} = 0 \}\) is a singleton \(\{ E^{*} \}\), which is the endemic equilibrium of the model system (5). Then, by LaSalle’s invariant principle, it implies that \(E^{*}\) is globally asymptotically stable in the interior of \(\Omega\) if \(Q .
Global sensitivity of model parameters on the state variables using LHS–PRCC
In this study, we utilize the Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficient (PRCC) methods to conduct uncertainty and global sensitivity analyses22. The LHS method samples model parameters while accounting for uncertainty by treating them as probabilistic variables that are uniformly distributed across specified ranges. The PRCC method calculates the correlation coefficients between model outputs (state variables) and model inputs (parameters), assessing the strength of their relationship. A positive PRCC value indicates that as the input parameter increases, the output state variable also tends to increase. The PRCC of the parameters with respect to the state variables \(S, E, I_T, I_C, I_{CT}, T, R, B_T,\) and \(B_C\) is computed using the basic model (5), LHS parameters, a large sample size of \(N = 1000\), initial conditions for simulating the basic model, and an algorithm from23. This relationship is illustrated in the following figures.
PRCC values for LHS parameters of model (5) in the typhoid infected human population.
The results in Fig. 2 indicate that \(\alpha\), \(\tau\), \(n\), \(\Lambda\), and \(g_1\) are positively correlated to the infected Typhoid fever humans \(I_T\), meaning that they are responsible for increasing the number of infected humans \(I_T\). All these parameters should be reduced throughout the epidemic of Typhoid fever, with the exception of shedding behavior and loss of immunity, which should be reduced during the first twenty (20) days. Through various control measures such as sanitation practices, boosting individuals’ immune systems, educational campaigns, and safe environmental and food management systems, efforts should be made to reduce the number of humans infected by Typhoid fever. On the other hand, \(\Psi _{1}\), \(k_1\), and \(\mu _{1}\) are inversely proportional to \(I_T\), meaning that they lower the number of infected humans \(I_T\) whenever they increase.
PRCC values for LHS parameters of model (5) in the COVID-19 infected human population.
Figure 3 indicates that \(\alpha\), \(\tau\), \(g_C\), \(m\), \(\beta _3\), and \(\Lambda\) are significantly positively correlated to \(I_C\). These parameters are responsible for the increase in the number of SARS-CoV-2 infections in the human population when they increase. The influence of parameters \(\alpha\), \(\tau\), \(g_C\), \(m\), \(\beta _3\), and \(\Lambda\) is high at the beginning of the epidemic and decreases slightly, stabilizing after 20 days. The influence of \(\Lambda\) increases with time because the more individuals are recruited, the more interactions occur, resulting in more physical transmission among humans. On the other hand, parameters such as \(k_3\), \(\mu _1\), \(g_3\), and \(k_C\) are significantly negatively correlated with \(I_C\), implying that they are responsible for lowering the infections among the human population when they increase. Control measures that speed up recovery rates among humans are ideal for reducing human infections. The parameters such as \(r_3\), \(g_4\), and \(\mu _3\) are observed to have little impact on the dynamics of infected humans.
PRCC values for LHS parameters of model (5) in the typhoid, and COVID-19 co-infected human population.
Figure 4 shows that \(\alpha\), \(\tau\), \(\beta _3\), \(\Lambda\), and \(g_C\) are directly proportional to the number of \(I_{CT}\). The increase or decrease of these parameters at specific times influences the increase or decrease in the number of co-infection cases. All these parameters should be reduced throughout the co-infection epidemics, except for the contact rate of COVID-19 through human-to-human transmission \(\beta _3\) and the transition rate \(\tau\), which should be reduced during the first 18 days and 30 days, respectively. To achieve this, control measures such as vaccinations, water sanitation, and treatment should be introduced to reduce the spread of co-infection diseases. On the other hand, \(m\), \(\mu _1\), \(\mu _2\), \(k_2\), and \(K_C\) are inversely proportional to \(I_{CT}\). This implies that an increase or decrease in these parameters results in a decrease or increase in co-infection cases.
