Existence, finiteness and non-negativity of solutions

Positivity of solution

Theorem 1: If $\Omega = \{ {( {J, \;K, \;L} ) \in R_ + ^3 \colon J( 0 ) > 0, \;K( 0 ) \ge 0, \;L( 0 ) \ge 0} \}$, then the solution set {J(t), K(t), L(t)} of the system of equation (1) is non-negative for all t ≥ 0.

Proof: We assume that $\tau = \sup \{ {t \ge 0\;\colon \;J_0( u ) > 0, \; \;K_0( u ) \ge 0, L_0( u ) \ge 0, \;\;\forall u\in [ {0, \;\;t} ] } \}$. Since J ₀(t) > 0, K ₀(t) ≥ 0 and L ₀(t) ≥ 0, then τ ≥ 0. If τ < 0, then automatically J ₀(t) or K ₀(t) or L ₀(t) = 0 at τ. Now, use the system of equation (1)

(3)$$\displaystyle{{dJ} \over {dt}} = b-\displaystyle{{\psi JK} \over N} + \omega L-\rho J.$$

The answer to equation (3) at time t is thus given by applying a variant of the constant formula.

(4)$$\eqalign{J( \tau ) &= K( 0 ) e^{-\int _0^T ( {( \psi K/N) + \rho } ) ( a ) dJ} + \int _0^t ( {b + \omega L} ) \cr&\quad\times e^{-\int _J^t ( {( \psi K/N) + \rho } ) ( u ) du}dJ > 0.}$$

Furthermore, because of the non-negative in [0, τ], J(τ) > 0. Similar evidence demonstrates that K(τ) ≥ 0, and L(τ) ≥ 0, this leads to contradiction. Thus, τ = ∞. As a result, for t ≥ 0, the solutions are positive.

Invariant region

Lemma 1: Suppose that the following system's initial conditions are all positive in $R_ + ^3$.

(5)$$\Omega = \left\{{( {J, \;\;K, \;\;L} ) \in R_ + ^3 ; \;\;0 \le J( t ) + K( t ) + L( t ) \;\colon \;N( t ) \le \displaystyle{b \over \rho }} \right\}, \;$$

where Ω is essentially invariant.

Proof: Finding the derivatives of N(t) = J(t) + K(t) + L(t) by considering time as independent variable and using equation (1), we have

(6)$$\displaystyle{{dN} \over {dt}} = b-\rho N.$$

After some simplification,

(7)$$N( t ) = \displaystyle{b \over \rho } + \left({N( 0 ) -\displaystyle{b \over \rho }} \right)e^{-\rho t}.$$

Then, $\lim \mathop {\sup }\limits_{t\to \infty } N( t) \le ( b/\rho )$. Hence, Ω is positively invariant. Therefore, all solutions of the system given in (1) including initial condition are in Ω.

Disease-free critical points

Disease-free critical points (DFE), E ₀ are solutions that reach a stable state, in which the population is free of disease. In the absence of disease, we have K = 0. Then, the disease-free equilibrium points (DFE), E ₀, are given by:

(8)$$E_0 = ( {J^0, \;\;K^0, \;\;L^0} ) = \left({\displaystyle{b \over \rho }, \;\;0, \;\;0} \right).$$

The basic reproduction number (R ₀): The basic reproduction number R ₀ of mathematical model of the system of equation (1) is determined using next-generation matrix technique [Reference van den Driessche and Watmough21]. Now, let x = (K, L, J), then the system of equation (1) can be rewritten as:

(9)$$\displaystyle{{dx} \over {dt}} = F( x ) -V( x ) $$

where

(10)$$F( x ) = \left({\matrix{ {\displaystyle{{\psi JK} \over N}} \cr 0 \cr 0 \cr } } \right), \;$$

and

(11)$$V( x ) = \left({\matrix{ {( {\sigma + \rho } ) K} \hfill \cr {-\sigma K + ( {\omega + \rho } ) L} \hfill \cr {-b + \displaystyle{{\psi JK} \over N}-\omega L + \rho J} \hfill \cr } } \right).$$

The Jacobian matrices of F(x) and V(x) at disease-free equilibrium point, E ₀ is given by:

(12)$$F = \left({\matrix{ \psi & 0 & 0 \cr 0 & 0 & 0 \cr 0 & 0 & 0 \cr } } \right), \;$$

and

(13)$$V = \left({\matrix{ {\sigma + \rho } & 0 & 0 \cr {-\sigma } & {\omega + \rho } & 0 \cr \psi & {-\omega } & \rho \cr } } \right).$$

Using the method of next-generation matrix, the fundamental reproduction number, R ₀ is the spectral radius of FV ⁻¹ or the dominant eigenvalue of FV ⁻¹ and thus, the fundamental reproduction number R ₀ is given by:

(14)$$R_0 = \displaystyle{\psi \over {\sigma + \rho }}.$$

Endemic critical point

Disease endemic critical point, E ₁, is said to be steady state if the disease continues in the population. If E ₁ = (J*, K*, L*) is disease equilibrium point, it satisfies the following algebraic equations:

(15)$$\eqalign{& b-\displaystyle{{\psi J^\ast K^\ast } \over {N^\ast }} + \omega L^\ast{-}\rho J^\ast { = } \,0, \;\cr & \displaystyle{{\psi J^\ast K^\ast } \over {N^\ast }}-( {\sigma + \rho } ) K^\ast{ = } \, 0, \;\cr & \sigma K^\ast{-}( {\omega + \rho } ) L^\ast{ = } \, 0.} $$

In the system (15) of the second equation, we obtain:

(16)$$J^\ast{ = } \displaystyle{b \over {\rho R_0}}.$$

Similarly, in system (15) of the third equation, we get:

(17)$$L^\ast{ = } \displaystyle{\sigma \over {\omega + \rho }}K^\ast.$$

Substituting equations (16) and (17) in system (15) of the first equation, we have:

(18)$$K^\ast{ = } \displaystyle{{b( {\omega + \rho } ) ( {R_0-1} ) } \over {\rho R_0( {\omega + \sigma + \rho } ) }}.$$

Then, after substituting equation (18) in equation (17), we obtain:

(19)$$R^\ast{ = } \displaystyle{{\sigma b( {R_0-1} ) } \over {\rho R_0( {\omega + \sigma + \rho } ) }}.$$

Local stability of the disease-free equilibrium point

The linearised form of the system of equation (1) at the steady state can be used to discuss the local stability of the disease-free equilibrium.

Theorem 3: If R ₀ < 1, the equilibrium solution of the non-linear system (1) is locally asymptotically stable.

Proof: At the point of equilibrium without disease, the Jacobian matrix of system given in (1) is:

(20)$$\eqalign{J (E_0) & = J (J^0,\;K^0,\;L^0{\rm )} = J\left( {\displaystyle{b \over \rho },\;0,\;0} \right) \cr & = \left( {\matrix{ {-\rho } & {-\psi } & \omega \cr 0 & {\psi -\left( {\sigma + \rho } \right)} & 0 \cr 0 & \sigma & {-\left( {\omega + \rho } \right)} \cr } } \right).} $$

The characteristic equation of (20) at disease-free equilibrium point E ₀ is

$$\vert {J( {E_0} ) -\lambda I_3} \vert = \left\vert {\matrix{ {-\rho -\lambda } & {-\psi } & \omega \cr 0 & {( {\psi -( {\sigma + \rho } ) } ) -\lambda } & 0 \cr 0 & \sigma & {-( {\omega + \rho } ) -\lambda } \cr } } \right\vert = 0, \;$$

(21)$$( {-\mu -\lambda } ) ( ( {\psi -( {\sigma + \rho } ) ) -\lambda } ) ( {-( {\omega + \psi } ) -\lambda } ) = 0.$$

Consequently,

(22)$$\lambda _1 = {-}\rho < 0, \;\lambda _2 = {-}( {\omega + \rho } ) < 0, \;$$

(23)$$\lambda _3 = \psi -( {\sigma + \rho } ) = ( {\sigma + \rho } ) \left({\displaystyle{\psi \over {\sigma + \rho }}-1} \right) = ( {\sigma + \rho } ) ( {R_0-1} ).$$

Here, λ ₃ = (σ + ρ)(R ₀ − 1) < 0 if R ₀ < 1 implying that all the eigenvalues are negative. Thus, disease-free equilibrium point E ₀ is locally asymptotically stable.

Global stability of the disease-free equilibrium point

Theorem 4: Suppose that R ₀ < 1, then the disease-free equilibrium point E ₀ is globally asymptotically stable.

Proof: To prove the global stability of the disease-free equilibrium point E ₀, we need to construct the following Lyapunov function:

(24)$$G( {J, \;\;K, \;\;L} ) = \displaystyle{1 \over 2}( {( {J-J^0} ) + ( {K-K^0} ) } ) ^2.$$

Clearly, (G(J, K, L)) ≥ 0 at disease-free equilibrium point and equal to zero at J = J ⁰ and K = K ⁰. Then, the derivative of equation (24) with respect to time t becomes:

(25)$$\displaystyle{d \over {dt}}( {G( {J, \;\;K, \;\;L} ) } ) = ( {J-J^0 + K-K^0} ) \left({\displaystyle{{dJ} \over {dt}} + \displaystyle{{dJ} \over {dt}}} \right).$$

Substituting the values of (dJ/dt) and (dK/dt) from the system of equation (1), we have

(26)$$\eqalign{ \displaystyle{d \over {dt}}{\rm (}G{\rm (}J,\;K,\;L)) & = (J-J^0 + K-K^0)(b + \omega L-\rho J-{\rm (}\sigma + \rho )K),\; \cr & = {\rm }-{\rm (}J-J^0 + K-K^0{\rm ) (}D-E{\rm )}.} $$

Clearly, (d/dt)(G(J, K, L)) ≤ 0 if and only if D > E, where D = ρJ + (σ + ρ)K and E = b + ωL. Furthermore, (d/dt)(G(J, K, L)) = 0 if and only if J = J ⁰ and K = K ⁰. Thus, by the invariance principle of LaSalle [Reference Lasalle22], the disease-free equilibrium point E ₀ is globally asymptotically stable.

Local stability of the endemic equilibrium point

Here, to illustrate the local stability of the endemic disease equilibrium state, the Jacobian stability method is used.

Theorem 5: The system of (1) has locally asymptotically stable equilibrium solution E ₁, when R ₀ > 1.

The Jacobian matrix of the system given in (1), at the disease endemic equilibrium point, is:

(27)$$J( {E_1} ) = \left({\matrix{ {-\left({\displaystyle{{\psi K^\ast } \over N} + \rho } \right)} & {-\displaystyle{{\psi J^\ast } \over N}} & \omega \cr {\displaystyle{{\psi K^\ast } \over N}} & {-\left({\sigma + \rho -\displaystyle{{\psi J^\ast } \over N}} \right)} & 0 \cr 0 & \sigma & {-( {\omega + \rho } ) } \cr } } \right).$$

The characteristic equation of (27) at disease endemic equilibrium point E ₁

$$|J(E_1)-\lambda I_3| = \left| {\matrix{ {-\left( {\displaystyle{{\psi K^*} \over N} + \rho } \right)-\lambda } & {-\displaystyle{{\psi J^*} \over N}} & \omega \cr {\displaystyle{{\psi K^*} \over N}} & {-\left( {\sigma + \rho -\displaystyle{{\psi J^*} \over N}} \right)-\lambda } & 0 \cr 0 & \sigma & {-(\omega + \rho )-\lambda } \cr } } \right|.$$

That is

(28)$$\eqalign{& \left({-\left({\displaystyle{{\psi K^\ast } \over N} + \rho } \right)-\lambda } \right)\left({-\left({\sigma + \rho -\displaystyle{{\psi J^\ast } \over N}} \right)-\lambda } \right)( {-( {\omega + \rho } ) -\lambda } ) \cr & \quad + \left({\displaystyle{{\psi J^\ast } \over N}} \right)\left({\displaystyle{{\psi K^\ast } \over N}} \right)( {-( {\omega + \psi } ) -\lambda } ) + \omega \sigma \left({\displaystyle{{\psi K^\ast } \over N}} \right) = 0.} $$

Now, equation (28) can be simplified as

(29)$$\lambda ^3 + a_1\lambda ^2 + a_2\lambda + a_3 = 0, \;$$

where

$$\eqalign{& a_1 = 3\rho + \sigma + \omega + \displaystyle{\psi \over {R_0}}\left( {\displaystyle{{ (\omega + \psi )(R_0-1{\rm )}} \over {\omega + \sigma + \rho }}-1} \right),\; \cr & a_2 = \displaystyle{{\psi (\omega + \psi )(2\mu + \sigma + \omega )(R_0-1)} \over {R_0 (\omega + \sigma + \rho )}} + \rho (\omega + \rho ),\; \cr & a_3 = \displaystyle{{\rho \psi } \over {R_0}} (\omega + \rho {\rm ) (}R_0-1{\rm )}.} $$

Thus, by Routh–Hurwitz criterion, the system of equation (1) is locally asymptotically stable if a ₁ > 0, a ₂ > 0 and a ₃ > 0. Hence, E ₁ is locally asymptotically stable if R ₀ > 1.

Global stability of disease endemic equilibrium point

For examining the global asymptotic stability of the disease endemic equilibrium point E ₁, the following model is used.

(30)$$\eqalign{& b = \displaystyle{{\psi J^\ast K^\ast } \over N}-\omega L^\ast{ + } \rho J^\ast , \;\cr & \displaystyle{{\psi J^\ast K^\ast } \over N} = ( {\sigma + \rho } ) K^\ast , \;\cr & \sigma K^\ast{ = } ( {\omega + \rho } ) L^\ast.} $$

Theorem 6: If R ₀ > 1, then the model given in (1) is globally asymptotically stable at E ₁ when ω = 0.

Proof: Using the method proposed by [Reference Fantaye and Birhanu23], we define the subsequent Lyapunov function for equation (1):

(31)$$V( t ) = J-J^\ast{-}J^\ast \ln \displaystyle{J \over {J^\ast }} + K-K^\ast{-}K^\ast \ln \displaystyle{K \over {K^\ast }} + \displaystyle{{\psi J^\ast } \over {\sigma N}}\left({L-L^\ast{-}L^\ast \ln \displaystyle{L \over {L^\ast }}} \right).$$

After differentiating equation (31) with respect to time t, we have

(32)$$\displaystyle{{dV} \over {dt}} = \left({1-\displaystyle{{J^\ast } \over J}} \right)\displaystyle{{dJ} \over {dt}} + \left({1-\displaystyle{{K^\ast } \over K}} \right)\displaystyle{{dK} \over {dt}} + \displaystyle{{\psi K^\ast } \over {\sigma N}}\left({1-\displaystyle{{L^\ast } \over L}} \right)\displaystyle{{dL} \over {dt}}.$$

Now,

(33)$$\eqalign{& \left( {1-\displaystyle{{J^{\ast}} \over J}} \right)\displaystyle{{dJ} \over {dt}} = \left( {1-\displaystyle{{J^{\ast}} \over J}} \right)\left[ {b-\displaystyle{{\psi JK} \over N} + \omega L-\psi J} \right],\; \cr & = \left( {1-\displaystyle{{J^{\ast}} \over J}} \right)\left[ {\displaystyle{{\psi J^{\ast}K^{\ast}} \over N}-\omega L^{\ast} + \rho J^{\ast}-\displaystyle{{\psi JK} \over N} + \omega L-\rho J} \right],\; \cr & = \left( {1-\displaystyle{{J^{\ast}} \over J}} \right)\left[ {\displaystyle{{\psi J^{\ast}K^{\ast}} \over N}\left( {1-\displaystyle{{JK} \over {J^{\ast}K^{\ast}}}} \right)-\omega L^{\ast}\left( {1-\displaystyle{L \over {L^{\ast}}}} \right)-\rho J\left( {1-\displaystyle{{J^{\ast}} \over J}} \right)} \right],\; \cr & \le \left[ {\displaystyle{{\psi J^{\ast}K^{\ast}} \over N}\left( {1-\displaystyle{{JK} \over {J^{\ast}K^{\ast}}}} \right)\left( {1-\displaystyle{{J^{\ast}} \over J}} \right)-\omega L^{\ast}\left( {1-\displaystyle{L \over {L^{\ast}}}} \right)\left( {1-\displaystyle{{J^{\ast}} \over J}} \right)} \right],\; \cr & = \displaystyle{{\psi J^{\ast}K^{\ast}} \over N}\left[ {1-\displaystyle{{J^{\ast}} \over J}-\displaystyle{{JK} \over {J^{\ast}K^{\ast}}} + \displaystyle{K \over {K^{\ast}}}} \right].} $$

(34)$$\eqalign{ \left( {1-\displaystyle{{K^{\ast}} \over K}} \right)\displaystyle{{dK} \over {dt}} & = \left( {1-\displaystyle{{K^{\ast}} \over K}} \right)\left[ {\displaystyle{{\psi JK} \over N}-{\rm (}\sigma + \rho )K} \right],\; \cr & = \left( {1-\displaystyle{{K^{\ast}} \over K}} \right)\left[ {\displaystyle{{\psi JK} \over N}-\displaystyle{{\psi J^{\ast}K} \over N}} \right],\; \cr & = \displaystyle{{\psi J^{\ast}K^{\ast}} \over N}\left[ {1 + \displaystyle{{JK} \over {J^{\ast}K^{\ast}}}-\displaystyle{J \over {J^{\ast}}}-\displaystyle{K \over {K^{\ast}}}} \right].} $$

(35)$$\eqalign{ \displaystyle{{\psi J^{\ast}} \over {\sigma N}}\left( {1-\displaystyle{{L^{\ast}} \over L}} \right)\displaystyle{{dL} \over {dt}} & = \displaystyle{{\psi J^{\ast}} \over {\sigma N}}\left( {1-\displaystyle{{L^{\ast}} \over L}} \right){\rm [}\sigma K- (\omega + \rho )C],\; \cr & = \displaystyle{{\psi J^{\ast}} \over {\sigma N}}\left( {1-\displaystyle{{L^{\ast}} \over L}} \right)\left[ {\sigma K-\sigma K^{\ast}\displaystyle{L \over {L^{\ast}}}} \right],\; \cr & = \displaystyle{{\psi J^{\ast}K^{\ast}} \over N}\left[ {1 + \displaystyle{K \over {K^{\ast}}}-\displaystyle{L \over {L^{\ast}}}-\displaystyle{{KL^{\ast}} \over {K^{\ast}L}}} \right].} $$

When the outcomes of equations (33)–(35) are substituted to equation (32), we obtain

(36)$$\eqalign{ \displaystyle{{dV} \over {dt}} & = -\omega L^{\ast}\left( {1-\displaystyle{L \over {L^{\ast}}}} \right)\left( {1-\displaystyle{{J^{\ast}} \over J}} \right)\quad \cr & + \displaystyle{{\psi J^{\ast}K^{\ast}} \over N}\left[ {3-\displaystyle{{J^{\ast}} \over J}-\displaystyle{J \over {J^{\ast}}}-\displaystyle{K \over {K^{\ast}}}\left( {\displaystyle{L \over {L^{\ast}}}-1} \right)-\displaystyle{L \over {L^{\ast}}}} \right].} $$

Here, (dV/dt) ≤ 0 if [3 − (J*/J) − (J/J*) − (K/K*)((L/L*) − 1) − (L/L*)] ≤ 0. Therefore, using [Reference Lasalle22], E ₁ is globally asymptotically stable whenever R ₀ > 1.

Sensitivity analysis of model parameters

Sensitivity analysis is commonly used for verifying and identifying parameters that can influence the basic reproduction number, R ₀, which determines the robustness of model forecasts. It indicates the importance of each parameter for the transmission of disease. According to [Reference Rodrigues, Monteiro and Torres24], sensitivity indices authorise us to quantify how much a variable varies when a parameter is changed. When the variable is a differentiable function of parameters, partial derivatives can also be used for constructing the sensitivity index.

Definition: The normalised forward sensitivity index of an R ₀, which depends differentially on a parameter, x _i, is given as:

(37)$$\Pi _{x_i}^{R_0} = \displaystyle{{\partial R_0} \over {\partial x_i}}.\displaystyle{{x_i} \over {R_0}}, \;$$

where x _i represents all the basic parameters and R ₀ = (ψ/σ + ρ).

(38)$$\eqalign{& \Pi _\psi ^{R_0} = \displaystyle{{\partial R_0} \over {\partial \psi }}.\displaystyle{\psi \over {R_0}} = 1 > 0, \;\cr & \Pi _\sigma ^{R_0} = \displaystyle{{\partial R_0} \over {\partial \sigma }}.\displaystyle{\sigma \over {R_0}} = {-}\displaystyle{\sigma \over {\sigma + \rho }} < 0, \;\cr & \Pi _\rho ^{R_0} = \displaystyle{{\partial R_0} \over {\partial \rho }}.\displaystyle{\rho \over {R_0}} = {-}\displaystyle{\rho \over {\sigma + \rho }} < 0.} $$

Interpretation of sensitivity indices

The sensitivity indices of R ₀ with respect to the key parameters are shown on Table 2. The positive indices ψ show that it has a significant impact on the spread of the disease with increasing values. Since as their values rise, the basic reproduction number increases, and hence the average number of secondary cases of infection rises as well. In addition, those parameters with negative sensitivity indices σ and ρ have the effect of reducing illness pain when their values increase while the others unchanged. Furthermore, as their values increase, the basic reproduction number decreases, which results in the reduction of the disease's endemic areas.

Table 2. Sensitivity indices of the model's parameters