Mathematical model and stability analysis on the transmission dynamics of skin sores

Abstract

In this study, a non-linear deterministic model for the transmission dynamics of skin sores (impetigo) disease is developed and analysed by the help of stability of differential equations. Some basic properties of the model including existence and positivity as well as boundedness of the solutions of the model are investigated. The disease-free and endemic equilibrium were investigated, as well as the basic reproduction number, R₀, also calculated using the next-generation matrix approach. When R₀ < 1, the model's stability analysis reveals that the system is asymptotically stable at disease-free critical point globally as well as locally. If R₀ > 1, the system is asymptotically stable at disease-endemic equilibrium both locally and globally. The long-term behaviour of the skin sores model's steady-state solution in a population is investigated using numerical simulations of the model.

Introduction

Skin sores (impetigo) is a highly contagious, superficial skin infection that most commonly affects children [1], which increases in prevalence in late summer [2]. Streptococcus pyogenes (group A Streptococcus) is one of the most important bacterial causes of skin and soft tissue infections worldwide [3]. Impetigo most typically affects exposed parts, such as the face and extremities. The lesions are well-localised, however they are frequently numerous and might be bullous or non-bullous [4].

Staphylococcus aureus and Streptococcus pyogenes produce this epidermal infection [5]. Although Staphylococcus aureus is the most common cause of impetigo in developed countries, S. pyogenes is still a major cause in developing countries [6]. Infections of the skin are caused primarily by group A streptococci (GAS) [7]. GAS infections have the potential to become invasive and can lead to serious post-infectious consequences of glomerulonephritis and rheumatic heart disease [8].

The age-specific prevalence of skin sores among children of up to 14 years of age is high, but remarkably consistent across the age groups studied [9]. Impetigo has increased occurrence in close contact, warm and humid environments [10]. In vitro investigations are used to assess the efficacy of prospective novel antimicrobial agents, with in vivo studies in animal models and/or people following if they are successful [11]. According to recent estimates, impetigo affects anywhere from 111 million children in impoverished nations to 140 million people worldwide at any given moment [12].

Impetigo is divided into non-bullous (also known as impetigo contagiosa) and bullous types. Bullous impetigo is caused by a staphylococcal toxin and does not require a host response to manifest clinical disease, whereas non-bullous impetigo is an infection-related host response [1] (Figs 1 and 2).

Fig. 1. Non-bullous impetigo.

Fig. 2. Bullous impetigo.

In recent years, mathematical modelling is becoming key techniques in the area of infectious disease propagation and control [13–18]. The studies on the mathematical analysis of human diseases and epidemic models have been conducted, which combine numerical studies with dynamic system methodologies such as stability analysis, LaSalles invariance principle, Routh–Hurwitz criterion and Lyapunov function [19].

Tanaka and Ono [20] provide an overview of how mathematical modelling can improve understanding of skin diseases, considering both their benefits and challenges. They discuss the significance of strong linkages between experimental, mathematical modelling and data analysis in order to effectively utilise the modelling technique to improve skin research in the post-genomic era.

Lydeamore et al. [7] constructed a stochastic representation of the susceptible-infectious-susceptible model. They have provided the first model-based estimates of skin sores infection duration (between 12 and 20 days), infection strength and baseline reproductive ratio in three different settings. By performing the estimation in a modelling framework, the interval-censored nature of the data was taken into account.

Motivated by the works of [7], for this manuscript, a non-linear deterministic model to study and analyse the transmission of skin sores is proposed. Moreover, we subdivide the total population into susceptible (J(t)), infected (K(t)) and recovered (K(t)) classes.

The manuscript is organised by four sections. In the first section, introduction to the transmission disease skin sores is presented. We construct a mathematical model for the dynamics of skin sore transmission in section ‘Model formulation’. Additionally, we go over the fundamental characteristics of the model and stability analysis of the disease-free and endemic equilibria. The topic of numerical simulation is covered in section ‘Numerical simulation’. At last, the conclusion is provided in section ‘Conclusion’.

Model formulation

In this section, we develop a new mathematical model for skin sores transmission dynamics by dividing the total population N(t) into three compartments: susceptible (J(t)), those who are at risk of becoming infected by the disease; infected (K(t)), those who have the pathogen in their organism and can transmit it; and recovered (L(t)), those who have recovered from the disease. The model also assumes that there is a positive recruitment rate b into susceptible class (J(t)) and positive natural death rate ρ for all time under the study. The susceptible individuals can become infected at rate (ψJ/N) through contact with infected individuals which are depending on time, where ψ is effective contact rate. The remaining model parameters are likewise defined as follows: σ is the recovery rate and ω is the rate of susceptibility of recovered individuals (Fig. 3, Table 1).

Fig. 3. Flow chart of the model.

Table 1. Parameters of the model

Based on the assumption and flow chart of the model, we have the following system of non-linear ordinary differential equations:

(1)$$\eqalign{& \displaystyle{{dJ} \over {dt}} = b-\psi \displaystyle{{JK} \over N} + \omega L-\rho J, \cr & \displaystyle{{dK} \over {dt}} = \psi \displaystyle{{JK} \over N}-{\rm (}\sigma + \rho )K, \cr & \displaystyle{{dL} \over {dt}} = {\rm }\sigma K-{\rm (}\omega + \rho {\rm )}L,} $$

with

(2)$$J( 0 ) = J_0 > 0, \;\;K( 0 ) = K_0 \ge 0, \;\;L( 0 ) = L_0 \ge 0.$$

Existence, finiteness and non-negativity of solutions

Existence of solution

The mathematical model of our system is well-posed. In actuality, Picard's (or Cauchy–Lipschitz) theorem states that its solution exists, is unique and is constantly dependent on the initial data. Because we are dealing with human population, we must guarantee that our solutions are positive and bounded.

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 [21]. 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 [22], 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 [23], 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 [22], 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 [24], 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

Numerical simulation

The numerical results for the system of (1) for various parameter values are presented. The simulation is done by using Matlab ODE45. The values of parameters are given in Table 3. We assumed that the total population and recruitment rate are related by b = ψN. The initial conditions given below serve as the foundation for the simulations and analyses: J(0) = 10 000, K(0) = 500 and L(0) = 2000.

Table 3. Parameter values of the model

From Figure 4, we observe that susceptible and recovered individuals asymptotically increase to the disease-free equilibrium point, while the infected individuals decrease asymptotically to the disease-free equilibrium point. Such conditions exist due to the fact that R ₀ = 0.8603 < 1. This supports that the stability of the disease-free equilibrium point exists when R ₀ < 1, which means if R ₀ < 1, for the duration of the disease period, a single infected individual on average produces less than one newly infected individual.

Fig. 4. Time series plot of state variables for R ₀ = 0.8603 < 1.

From Figure 5, we observe that susceptible and recovered individuals are decreased due to the influence of infected individuals, then they become infected; as a result, the infected individuals are increased. Therefore, infected individuals are increased and the disease endemic equilibrium point exists and stable. The existence of this condition is due to the fact that R ₀ = 2.0492 which is greater than one. This supports the theorem that the stability of disease endemic equilibrium point exists when R ₀ > 1, this means that when R ₀ > 1, on average, each infected person creates more than one other contaminated person, then infected individuals will be able to spread in the given society.

Fig. 5. Time series plot of state variables for R ₀ = 2.0492 > 1.

Figure 6 shows that as effective contact rate, ψ grows, the number of infected people rises while the number of susceptible people falls as a result of the influence of infected people.

Fig. 6. Variations of susceptible population J(t) and infected population K(t) w.r.t. time t for different values of ψ.

Furthermore from Figure 7, we observe that as the recovery rate σ increases, infected individuals decrease while the recovered individuals are increased.

Fig. 7. Variations of infected population K(t) and recovered population L(t) w.r.t. time t for different values of σ.

Also from Figure 8, we see that as the rate of susceptibility of recovered individuals decreases, recovered individuals are increased while the susceptible individuals are decreased.

Fig. 8. Variations of susceptible population J(t) and recovered population K(t) w.r.t. time t for different values of ω.

Conclusion

A non-linear deterministic mathematical model for skin sores transmission dynamics was constructed in this work. We first established that the model is well-posed epidemiologically and mathematically. The basic reproduction number (R ₀), as well as endemic and the disease-free critical points are determined. The value of (R ₀) is used to assess whether equilibrium points are stable or not. When R ₀ < 1, the disease-free critical point (E ₀) will be both globally and locally asymptotically stable. That is, on average, one infected individual produces less than one newly infected individual over the course of its disease period. As a result, the infected people are gradually removed from the society. Also, the disease endemic critical point (E ₁) is globally and locally asymptotically stable if R ₀ is greater than 1. That is, each infected individual produces on average more than one new infected individual, then infected individuals will be able to spread in the given society. Consequently, the disease persists in the society. The computational simulations and analyses demonstrate that when the effective contact rate declines and as the recovery rate increases, infected individuals decrease from the society. The model we adapted did not provide optimal control and cost-effectiveness of various intervention strategies that can be explored in the future to find out which strategy is the top in controlling skin sores.