Qualitative behavior of a highly non-linear Cutaneous Leishmania epidemic model under convex incidence rate with real data

1.   Introduction

Leishmaniasis, caused by Leishmania type parasites, represents a significant public health concern. Extensive research has focused on four distinct types of Leishmaniasis ^[1,2]. Among these, Cutaneous Leishmaniasis (CL) emerges as the most prevalent form, primarily attributed to L. tropica and L. major. The condition manifests as skin lesions, often affecting exposed body regions like the extremities and face. CL endemic regions encompass parts of South America, Europe, and Asia. The World Health Organization (WHO) classifies CL as a highly perilous global disease, given its rapid annual incidence growth, surpassing one million new cases. Ten countries collectively contribute to 70–75 percent of the global CL incidence: Peru, Costa Rica, North Sudan, Ethiopia, Syria, Iran, Brazil, Colombia, Algeria, and Afghanistan ^[3]. The transmission of CL primarily occurs through blood-sucking insects, with various fly species serving as the principal vectors. Worldwide, approximately 700 species of such vectors have been identified, including 37 in Pakistan. Notably, CL has recently surged in northwestern Pakistan, resulting in significant mortality rates, particularly in the Khyber Pakhtunkhwa (KPK) province. A previous CL outbreak in 1997 within this province was linked to an Afghan refugee camp, where migrants from Kabul, Afghanistan, afflicted by CL, facilitated the epidemic ^[4].

The vectors responsible for parasite transmission belong to the Phlebotomus genus. These sand flies exhibit a broad habitat range, from arid deserts to tropical rainforests. They exhibit multiple hosts, including dogs, chickens, humans, mammals, livestock, and vertebrates ^[5]. These sand flies typically measure 2–3mm in length and have a latent period ranging from three to seven days ^[6]. They can reach heights of up to 2.51 meters (8.3 feet) above ground level ^[7]. Utilizing animal blood for sand fly feeding has been shown to reduce the risk of accidental transmission of human blood-borne pathogens, rendering it a cost-effective alternative to maintaining and preparing animals as blood sources ^[8]. Extreme temperatures have been observed to adversely affect the fecundity of these species ^[9]. Furthermore, recent research has delved into various epidemiological models of Cutaneous Leishmaniasis, contributing to a more comprehensive understanding of the disease dynamics ^{[10,11,12,13,14,15]}.

The utilization of mathematical models stands as a fundamental tool in the comprehensive analysis of the transmission dynamics and management strategies for infectious diseases. The process of formulating a mathematical model inherently involves delineating key assumptions, defining relevant parameters, and specifying the variables that drive the system's behavior. Notably, Chaves et al. ^[16] have made significant contributions to this field by presenting a mathematical model designed to elucidate the dynamics of American Cutaneous Leishmaniasis. Within their research, they delved into the establishment of threshold conditions that govern the persistence of infection within the population, shedding light on critical factors influencing disease prevalence. Building upon this foundation, Chaves et al. have further extended their investigations to explore the intricate interplay between deforestation and vector-borne diseases in recent work ^[5,16]. Their research endeavors have illuminated the profound impact of environmental alterations on disease transmission dynamics, providing valuable insights for disease control strategies. In a distinct line of inquiry, Das et al. ^[17] have undertaken a detailed study focusing on the temporal dimension by examining the effect of time delays on the dynamics of Cutaneous Leishmaniasis. Their investigations have revealed how temporal factors can significantly influence disease patterns and persistence. Calzada et al. ^[18] contributed to our understanding of vector-borne diseases by investigating the compositional changes that occur in sandflies following the application of thermal fogging, a method often employed for vector control. This research has provided critical information on the efficacy and potential ecological impacts of vector control interventions. Furthermore, Chaves ^[19] and Bacaer and Guernaoui ^[20] have engaged in the study of seasonal models of Cutaneous Leishmaniasis. By considering the seasonal variations in disease transmission and vector populations, their work has advanced our understanding of how environmental and climatic factors can influence disease dynamics, thus informing the development of seasonally-tailored control strategies. In light of these diverse research endeavors, it is evident that a multitude of strategies are under consideration for the control and eventual eradication of various infectious diseases. These strategies encompass a broad spectrum of interventions, ranging from environmental management and vector control to the development of mathematical models that offer valuable insights into disease dynamics and guide evidence-based decision-making in public health.

In this research endeavor, we introduce a more intricate framework, building upon the foundation laid in ^[21]. We acknowledge the significant role of disease incubation and, consequently, introduce an 'exposed' class for both human and vector populations. This entails accounting for individuals who are susceptible to the disease but are in a latent phase of infection. In essence, we develop a model for CL denoted as $\mathsf{S}_h \mathsf{E}_h \mathsf{I}_h \mathsf{R}_h \mathsf{S}_v \mathsf{E}_v \mathsf{I}_v$, featuring a convex incidence rate. The core focus of our investigation centers on unraveling the dynamic behavior of the CL model. Leveraging center manifold theory, we unveil a forward bifurcation phenomenon within the model, leading to the emergence of multiple endemic equilibria. On an alternative note, when a unique endemic equilibrium prevails, we establish its global stability using the geometric approach pioneered by Li and Muldowney ^[22]. The significance of our findings lies in the demonstration of global asymptotic endemism in a system denoted as $\mathsf{S}_h \mathsf{E}_h \mathsf{I}_h \mathsf{R}_h \mathsf{S}_v \mathsf{E}_v \mathsf{I}_v$, demonstrating a convex incidence rate concerning the state variables $\mathsf{I}_h$ and $\mathsf{I}_v$. This achievement is particularly captivating, given that such functional characteristics do not typically satisfy the sufficient conditions for global stability as expounded in ^[23]. Subsequently, we perform a sensitivity analysis to discern the parameters most susceptible to influence, employing the direct differentiation method.

Our article is organized as follows: Section 1 offers an in-depth exploration of the CL epidemic model. In Section 2, we formulate the ACL epidemic model and ascertain the threshold value $R_0$ through the next-generation method. Sections 3 and 4 impose specific conditions on the threshold value to elucidate our proposed model's local and global stability characteristics. Section 5 delves into a comprehensive bifurcation analysis of the model. The process of parameter estimation and data fitting is undertaken in Section 6. Section 7 focuses on sensitivity analysis, pinpointing the parameters with the greatest influence. Numerical simulations and discussions of the solutions are presented in Section 8. Finally, Section 9 provides a brief conclusion summarizing the findings and future directions of the study.

2.   Model formulation

In this section, we introduce a comprehensive model to depict the dynamics of the CL epidemic recently formulated in ^[21]. The model comprises seven distinct classes, further categorized into four human population subclasses denoted as $\mathsf{S}_h(t)$, $\mathsf{E}_h(t)$, $\mathsf{I}_h(t)$, and $\mathsf{R}_h(t)$, representing susceptible, exposed, infected, and recovered individuals, respectively. Additionally, three vector population subclasses are included, namely $\mathsf{S}_v(t)$, $\mathsf{E}_v(t)$, and $\mathsf{I}_v(t)$, signifying susceptible, exposed, and infected vectors. The total human population, denoted as $N_h(t)$, is defined as the sum of these subclasses, i.e., $N_h(t) = \mathsf{S}_h(t) + \mathsf{E}_h(t) + \mathsf{I}_h(t) + \mathsf{R}_h(t)$. The mathematical model formulated in ^[21] is given as follows:

An intriguing example of a nonlinear incidence rate employed in this model is the convex incidence rate function denoted as $f(\mathsf{S}, \mathsf{I}) = \beta \mathsf{S} \mathsf{I}(1+\vartheta \mathsf{I})$, where $\beta$ and $\vartheta$ are positive constants. This incidence rate signifies an elevated infection rate over a short time span, arising from dual exposures. The term $\beta IS$ reflects infection transmission from a single contact, while the component $\beta \vartheta \mathsf{I}^2 \mathsf{S}$ accounts for the generation of new infections resulting from a double exposure. Notably, prior research has predominantly conducted stability analyses of epidemic models employing geometric approaches with nonconvex incidence rates with respect to the infected population $I$. However, for convex incidence rates, this approach often imposes overly restrictive conditions ^[24]. The parameter $\Gamma_h$ represents the influx of new individuals into the susceptible human class $\mathsf{S}_h$. Infection of susceptible humans occurs when they are bitten by infected sandflies at a rate denoted as $\alpha \Psi_1 \mathsf{I}_v(t) \mathsf{S}_h(t)[1+\vartheta \mathsf{I}_v(t)]$, where $\Psi_1$ represents the transmission probability of ACL from sandflies to humans, and $\alpha$ denotes the sandfly biting rate ^[25]. The parameter $\theta$ characterizes the rate at which uninfected exposed humans transition into the infectious state, while $\mathcal{K}_1$ signifies the rate at which exposed humans become infectious. Additionally, $\beta$ denotes the rate of natural recovery of humans from the infected class.

Infected humans, following bites from sandflies, transmit the infection to the sandflies at a rate specified as $\alpha \Phi_1 \mathsf{I}_h(t) \mathsf{S}_v(t)[1+\vartheta \mathsf{I}_h(t)]$. In this context, $\mathfrak{U}_h$ represents the natural death rate in humans, whereas $\mathfrak{U}_v$ designates the natural death rate in sandflies. The probability of CL transmission from humans to sandflies is encapsulated by the parameter $\Phi_1$. Furthermore, $\mathcal{K}_2$ characterizes the rate at which exposed vector classes transition to the infected state. The total vector population $N_v$ is defined as $N_v(t) = \mathsf{S}_v(t) + \mathsf{E}_v(t) + \mathsf{I}_v(t)$. Keeping in view of the above discussion, we reformulate model (2.1) with a convex incidence rate. The complete mathematical representation of the model is given by

with

2.1.   Basic properties of the model

2.1.1.   Invariant region analysis

The total population dynamics are described by

In the context of this study, we define a biologically feasible region, denoted as $\Delta$, which is described by the conditions

Here, $\Delta$ encompasses a set of non-negative real values for the seven variables representing different population classes. Importantly, this region is bounded by the constraints that the total number of vectors ($N_v$) should not exceed $\frac{\Gamma_v}{\mathfrak{U}_v}$, and the total number of humans ($N_h$) should not exceed $\frac{\Gamma_h}{\mathfrak{U}_h}$.

Furthermore, based on the solutions to Eqs (2.4) and (2.5), it can be established that

This observation indicates that as time progresses towards infinity, the total human population ($N_h$) tends to stabilize at $\frac{\Gamma_h}{\mathfrak{U}_h}$, and similarly, the total vector population ($N_v$) stabilizes at $\frac{\Gamma_v}{\mathfrak{U}_v}$. Consequently, it can be concluded that the model is well-posed, and the region $\Delta$ is a positively invariant domain. In other words, the populations remain within this biologically plausible region over time, ensuring the model's validity and stability.

Lemma 1. The orthant $\mathbb{R}^7_+$ is positively invariant for the system described by (2.2).

Proof. Let $X$ represent all the state variables from $\mathsf{S}_h$ to $\mathsf{I}_v$, and assume

System (2.2) is expressed in the form

where

In this context, we observe that the matrix denoted as $L$ exhibits properties characteristic of a Metzler matrix, characterized by nonnegative elements in its off-diagonal entries, while matrix $B$ satisfies the condition $B \geq 0$. Consequently, we deduce that the dynamical system described by Eq (2.2) maintains positive invariance within the state space $\mathbb{R}^7_+$.

Lemma 2. If solutions exist, they are positive for all $t > 0$ under the initial conditions.

Proof. Let us assume that solutions exist in the interval $I$, for all $t\in I \subset [0, \infty)$. Now, consider the second equation of system (2.2) and its solution. This takes the following form:

We also take the third equation, and it has a solution as follows:

Clearly, it is evident from the preceding solutions that these variables exhibit strictly positive behavior. Similarly, it is possible to establish that $\mathsf{S}_h, \mathsf{R}_h, \mathsf{S}_v, \mathsf{E}_v$, and $\mathsf{I}_v$ exhibit solutions that remain non-negative.

2.2.   Basic reproductive number $R_0$

The disease-free equilibrium point of the system represented by Eq (2.2) is designated as $E^0$, which is given by

Considering the compartment comprising individuals in the states $(\mathsf{E}_h, \mathsf{I}_h, \mathsf{E}_v, \mathsf{E}_v)$ as the infected cohort, it can be deduced from the dynamics described by system (2.2) that

Utilizing the approach based on the next-generation matrix, we derive the Jacobian matrix denoted as $J$ for the aforementioned system when it is at the equilibrium point representing the absence of the disease. This matrix is expressed as

Now, breaking down matrix $J$ into its constituent components represented by $F$ and $V$, denoted as $J = F-V$, we obtain:

and

Hence, we have

Then, the characteristic equation becomes

The dominant eigenvalue gives us $R_0$, i.e.,

3.   Local stability

In this section, we establish the local stability properties of the system described by Eq (2.2) at both the disease-free equilibrium (DFE) point denoted as $E^0$ and the endemic equilibrium point denoted as $E^*$.

3.1.   Disease-free equilibrium point

Theorem 3.1. The system represented by Eq (2.2) exhibits asymptotic stability at the DFE point $E^0$ when $R_0 < 1$.

Proof. Within the framework of system (2.2), the Jacobian matrix at the equilibrium point $E^0$ is expressed as:

Subsequently, the characteristic equation pertaining to the matrix $J^{\mid0\mid}$ is formulated as follows:

This equation can be organized in the following manner:

where

Certainly, it is evident that three of the eigenvalues exhibit negativity, specifically, $\lambda_1 = -\mathfrak{U}_h$, $\lambda_2 = -\mathfrak{U}_h$, and $\lambda_3 = -\mathfrak{U}_v$. To analyze the remaining eigenvalues, we have:

We utilize the Routh-Hurwitz (RH) criterion to establish that the roots of the aforementioned equation have negative real parts, specifically satisfying condition $(H_1): a_i > 0$ where $i$ ranges from 1 to 4, and ensuring that $a_1a_2a_3 > a_3^2 + a_1a_4^2$. It is clear that $a_1 > 0, a_2 > 0, a_3 > 0$, and $a_4 > 0$, under the condition:

This inequality implies:

and subsequently, $a_1a_2a_3 > a^2_3+a^2_1a_4$. Hence, it can be observed that $(H_1)$ holds under the condition $R_0 < 1$. Consequently, the real parts of all eigenvalues are negative if and only if $R_0 < 1.$

3.2.   Endemic equilibrium point

System (2.2) is rearranged to express $\mathsf{S}^{*}_h, \mathsf{E}^{*}_h, \mathsf{I}^{*}_h, \mathsf{R}^{*}_h, \mathsf{S}^{*}_v,$ and $\mathsf{E}^{*}_v$ in terms of $\mathsf{I}^{*}_v$. Thus,

For $\Gamma_v \mathcal{K}_2 > (\mathfrak{U}_v + \mathcal{K}_2)\mathfrak{U}_v \mathsf{I}^*_v$, we can ensure the positivity of Eq (3.5). Moreover, through the substitution of $\mathsf{E}^*_v$ and $\mathsf{S}^*_v$ into

we derive

where

Using Descartes' Rule of Signs, one can easily come to the conclusion that Eq (3.7) exhibits two roots with distinct signs when the conditions $a_0 < 0$ and $a_1 > 0$ are met.

implies that

the term $\bigg(\frac{\mathcal{K}_1 \alpha \Psi_1 \Gamma_h}{(\mathcal{K}_1+\theta+\mathfrak{U}_h)(\beta+\mathfrak{U}_h) }\bigg)$ is derived from $R_0$. Therefore, in the case of $R_0 > 1$, a solitary positive equilibrium point can be identified within the model.

Theorem 3.2. If $R_0 \in \bigg(1, \frac{\sqrt{\Gamma_v}}{\sqrt{ \mathsf{S}_h^* \mathsf{S}_v^* (1+\vartheta \mathsf{I}_h^*)(1+2\vartheta \mathsf{I}^*_v)}}\bigg)$, then the endemic equilibrium $E^*$ of the system (2.2) is locally asymptotically stable.

Proof. For system (2.2), the Jacobian matrix at $E^*$ is given by

where $\mathcal{A}_{11} = (a \Psi_1 \mathsf{I}^*_v(1+\vartheta \mathsf{I}^*_v)), \; \; \mathcal{A}_{17} = a \Psi_1 \mathsf{S}_h^* (1+2\vartheta \mathsf{I}^*_v), \; \; \mathcal{A}_{22} = (\mathcal{K}_1+\theta+\mathfrak{U}_h)$, and $\mathcal{A}_{55} = (a \Phi_1 \mathsf{I}^*_h(1+\vartheta \mathsf{I}^*_h))$.

Indeed, a single eigenvalue of the Jacobian matrix $J^{\mid*\mid}$ is $\lambda_1 = -\mathfrak{U}_h$. The resulting matrix simplifies to:

After simplification, we get

where

The eigenvalues of $J^{\mid*\mid}_2$ take the form

and

and $\lambda_7$ is negative if

4.   Global asymptotic stability

In this section, we examine the global stability analysis of system (2.2) concerning both the disease-free (DFE) and endemic equilibria. We apply Castillo Chavez's method ^[26] to establish global stability at the disease-free equilibrium, and we employ a generalized form of Lyapunov theory ^[22] to analyze the global stability of the endemic equilibrium.

4.1.   Disease free equilibrium point

For the model described in Eq (2.2), global stability at the disease-free equilibrium point is established by employing the Castillo Chavez approach ^[26]. We shall now outline the Castillo Chavez approach ^[26]. Consider the two-dimensional system

Here, $\mathfrak{R}_1 = (\mathsf{S}_h, \mathsf{S}_v, \mathsf{R}_h)\in R^{3}$ and $\mathfrak{R}_2 = (\mathsf{I}_h, \mathsf{I}_v, \mathsf{E}_h, \mathsf{E}_v)\in R^{4}$.

$1)$ If $\frac{d\mathfrak{R}_1}{dt} = G(\mathfrak{R}_1, 0)$, then we can conclude that the DFE point, denoted as $\mathfrak{R}^{0}_1$, exhibits global asymptotic stability.

$2)$ $H(\mathfrak{R}_1, \mathfrak{R}_2) = B\mathfrak{R}_2-\bar{H}(\mathfrak{R}_1, \mathfrak{R}_2)$, where $\bar{H}(\mathfrak{R}_1, \mathfrak{R}_2)\geq0$ for $(\mathfrak{R}_1, \mathfrak{R}_2)\in\vartheta,$ and $B_{\mathfrak{R}_2}$ is an $M$-matrix.

Lemma 3. Let $R_0 < 1$. Then, if the aforementioned conditions are met, the equilibrium point $E^0 = (\mathfrak{R}^{0}_1, 0)$ of system (2.2) is globally asymptotically stable ^[26].

We now apply the above techniques to establish the global stability of the model represented by Eq (2.2) at the disease-free equilibrium. Consequently, we arrive at the following stability result.

Theorem 4.1. If $R_0 < 1$, then system (2.2) is globally asymptotically stable at the disease-free equilibrium point $E^0$, and unstable otherwise.

Proof: Let

By employing the system described in Eq (2.2), we obtain

For $\mathsf{S}_h = S^0_h$, $\mathsf{S}_v = S^0_v$, and $G(\mathfrak{R}_1, 0) = 0$ we obtain

Hence, as time approaches infinity ($t\rightarrow \infty$), it follows from Eq (4.4) that $\mathfrak{R}_1$ converges to $\chi^{0}_1$. Therefore, $\mathfrak{R}_1 = \chi^{0}_1$ represents a state of global asymptotic stability. Now,

As the total population is confined within the bounds of $S^0_h$ and $S^0_v$, i.e., $\mathsf{S}_h \leq \mathsf{S}^0_h$ and $\mathsf{S}_v \leq S^0_v$, it follows that $\alpha \Psi_1 \mathsf{S}^0_h \mathsf{I}^0_v \geq \alpha \Psi_1 \mathsf{I}_v \mathsf{S}_h [1+\vartheta \mathsf{I}_v]$ and $\alpha \Phi_1 \mathsf{S}^0_v \mathsf{I}^0_h \geq \alpha \Phi_1 \mathsf{I}_h \mathsf{S}_v [1+\vartheta \mathsf{I}_h]$. Consequently, it can be inferred that $\bar{H}(\mathfrak{R}_1, \mathfrak{R}_2) \geq 0$. Clearly, matrix $B$ satisfies the properties of an $M$-matrix. Therefore, both conditions are satisfied. As per Lemma 3, it can be concluded that the DFE point $E^0$ is globally asymptotically stable.

4.2.   Global stability of endemic equilibrium

For the assessment of global stability at the endemic equilibrium point $E^*$ in system (2.2), a geometric approach, as described in reference ^[22], is employed. The methodology can be summarized as follows: To explore the conditions under which $E^*$ attains global asymptotic stability, consider a differential equation

In this equation, $U \subset \mathbb{R}^{n}$ is an open, simply connected set, and $f: U \rightarrow \mathbb{R}^{n}$ is a function that is continuously differentiable $C^{1}(U)$. Assuming $f(x^{*}) = 0$ represents a solution of Eq (4.6), and for $x(t, x_0)$, the following conditions hold:

$3)$ There exists a compact absorbing set $K\in U$.

$4)$ System (4.6) has a unique equilibrium.

The equilibrium point denoted as $x^{*}$ is considered to be globally asymptotically stable within the set $U$ when it exhibits local asymptotic stability, and all trajectories within the set $U$ ultimately converge to the equilibrium point $x^{*}$. For $n\geq 2$, it is necessary to impose a condition on the function $f$ to avoid the presence of non-constant periodic solutions in Eq (4.6). This condition is commonly referred to as the Bendixson criteria. The traditional Bendixson criteria, given by $div f(x) < 0$ for $n = 2$, maintains its validity even under $C^{1}$ conditions ^[27]. Additionally, a point $x_0$ within the set $U$ is defined as "wandering" for Eq (4.6) if there exists a neighborhood $N$ of $x_0$ and a positive time value $\tau$, such that $N \cap x(t, N)$ is empty for all $t > \tau$. Hence, the following principle of global stability is established for autonomous systems in any finite dimension:

Lemma 4. Given the fulfillment of conditions (3) and (4) and the adherence to the Bendixson criterion in the context of (DFE) (4.6), wherein these conditions exhibit resilience under $C^{1}$ local perturbations of the function $f$ at all non-equilibrium and non-wandering points within the domain of (DFE) (4.6), we establish that the equilibrium point $x^{*}$ achieves global asymptotic stability within the set $U$, provided that it initially demonstrates local stability.

Next, we define a matrix-valued function $P$ on $U$ as

Equation (4.7) represents a matrix-valued function on $U$. We further assume that $P^{-1}$ exists and is continuous for $x\in K$. Now, we define a quantity $\bar{q}$ as

where $J^{[2]}$ is the second additive compound matrix of $J$, i.e., $J(x) = Uf(x)$, and $B = P_fP^{-1}+PJ^{[2]}P^{-1}$. Let $\ell(B)$ be the Lozinskii measure of matrix $B$ with respect to the norm $\|.\|$ in $R^{n}$, defined as

Therefore, if $\bar{q} < 0$, this implies the absence of any orbit that gives rise to a simple closed rectifiable curve, such as periodic orbits and heterocyclic cycles.

Lemma 5. If $U$ is simply connected and conditions (3) and (4) are satisfied, then the unique equilibrium $x^{*}$ of (DFE) (4.6) is globally asymptotically stable in $U$ if $\bar{q} < 0$ ^[22].

According to ^[27], $\bar{q}$ in Lemma 3 can be evaluated as

where $D_+$ represents the right-hand derivative. Now, we apply the above techniques to demonstrate the global stability of model (2.2) at its endemic equilibrium. Thus, we have the following stability analysis.

Theorem 4.2. Assuming $R_0 > 1$, the sole endemic equilibrium denoted as $E^*$ achieves global asymptotic stability within the region $\Delta$ when the following inequalities hold true:

Proof: In order to establish the global asymptotic stability of the presented model (2.2) at the endemic equilibrium $E^*$, we will focus on the subsystem within (2.2), which can be described as:

Commencing with the Jacobian matrix $J$ associated with system (4.10),

The matrix representing the second additive compound is expressed as

where

Now, choose a function $Q(\chi)$ such that

The matrix $Q^{-1}(\chi)$ is given by

Subsequently, when we calculate the time derivative, denoted as $Q_f(\chi)$, we obtain

Hence,

and

Now, using

matrix $B = Q_fQ^{-1}+QJ^{\mid2\mid}Q^{-1}$ becomes

where

We consider the following norm on $\mathbb{R}^6$ ^[28]. Furthermore, if we consider $z$ as belonging to $\mathbb{R}^6$ with components $\mathfrak{M}_i$, where $i$ ranges from 1 to 6, then

and

Next, we have to find some $\tau > 0$ so that $D_+\|z\| \leq -\tau z$. If $z$ is satisfies this inequality, then, by linearity, $-z$ will also satisfy this inequality. As mentioned in ^[29], we subdivide our proof into all sixteen different, possible cases, depending upon the definition of the norm (4.19). The following inequalities are incorporated for accomplishing the required analysis:

and

Case 1: $U_1 > U_2, \mathfrak{M}_1, \mathfrak{M}_2, \mathfrak{M}_3 > 0$, and $| \mathfrak{M}_1| > | \mathfrak{M}_2| + | \mathfrak{M}_3|$. Then,

so that

Case 2: $U_1 > U_2, \; \; \; \; \; \mathfrak{M}_1, \mathfrak{M}_2, \mathfrak{M}_3 > 0$, and $| \mathfrak{M}_1| < | \mathfrak{M}_2| + | \mathfrak{M}_3|$. Then,

so that

Case 3: $U_1 > U_2, \; \; \; \; \; \; \mathfrak{M}_1 < 0, \mathfrak{M}_2, \mathfrak{M}_3 > 0$, and $| \mathfrak{M}_1| > | \mathfrak{M}_2|$. Then,

so that

Since $| \mathfrak{M}_1| > | \mathfrak{M}_2|$, we then have

Case 4: $U_1 > U_2, \; \; \; \; \; \mathfrak{M}_1 < 0, \mathfrak{M}_2, \mathfrak{M}_3 > 0$, and $| \mathfrak{M}_1| < | \mathfrak{M}_2|$. Then,

so that

Case 5: $U_1 > U_2, \; \; \; \; \; \mathfrak{M}_1, \mathfrak{M}_2 > 0, \mathfrak{M}_3 < 0$, and $| \mathfrak{M}_2| > | \mathfrak{M}_1| + | \mathfrak{M}_3|$. Then,

so that

Case 6: $U_1 > U_2, \; \; \; \; \; \; \mathfrak{M}_1, \mathfrak{M}_2 > 0, \mathfrak{M}_3 < 0$, and $| \mathfrak{M}_2| < | \mathfrak{M}_1| + | \mathfrak{M}_3|$. Then,

so that

Since $-\bigg(a \Phi_1 \mathsf{I}_h(1+\vartheta \mathsf{I}_h)\bigg)| \mathfrak{M}_2| < 0$, then

Case 7: $U_1 > U_2, \; \; \; \; \; \mathfrak{M}_1, \mathfrak{M}_3 > 0, \mathfrak{M}_2 < 0$, and $| \mathfrak{M}_1| > max \{| \mathfrak{M}_2|, | \mathfrak{M}_3|\}$. Then,

so that

Case 8: $U_1 > U_2, \; \; \; \; \; \mathfrak{M}_1, \mathfrak{M}_3 > 0, \mathfrak{M}_2 < 0$, and $| \mathfrak{M}_2| > max \{| \mathfrak{M}_1|, | \mathfrak{M}_3|\}$. Then,

so that

Case 9: $U_1 > U_2, \; \; \; \; \; \mathfrak{M}_1, \mathfrak{M}_3 > 0, \mathfrak{M}_2 < 0$, and $| \mathfrak{M}_3| > max \{| \mathfrak{M}_1|, | \mathfrak{M}_2|\}$. Then,

so that

Since $-(a \Phi_1 \mathsf{I}_h(1+\vartheta \mathsf{I}_h))| \mathfrak{M}_2| < 0$, we then have

Case 10: $U_1 < U_2, \; \; \; \; \; \mathfrak{M}_4, \mathfrak{M}_5, \mathfrak{M}_6 > 0$. Then,

so that

Since $| \mathfrak{M}_2|+| \mathfrak{M}_3| < U_1 < U_2 < | \mathfrak{M}_4|+| \mathfrak{M}_5|+| \mathfrak{M}_6|$, we then have

Case 11: $U_1 < U_2, \; \; \; \; \; \mathfrak{M}_4, \mathfrak{M}_5 > 0, \mathfrak{M}_6 < 0$, and $| \mathfrak{M}_5| > | \mathfrak{M}_6|$. Then,

so that

Since $| \mathfrak{M}_2|+| \mathfrak{M}_3| < U_1 < U_2 < | \mathfrak{M}_4|+| \mathfrak{M}_5|$, we then have

Case 12: $U_1 < U_2, \; \; \; \; \; \mathfrak{M}_4, \mathfrak{M}_5 > 0, \mathfrak{M}_6 < 0$, and $| \mathfrak{M}_5| < | \mathfrak{M}_6|$. Then,

so that

Since $| \mathfrak{M}_2| < U_1 < U_2 < | \mathfrak{M}_4|+| \mathfrak{M}_6|$, we then have

Case 13: $U_1 < U_2, \; \; \; \; \; \mathfrak{M}_4, \mathfrak{M}_6 > 0, \mathfrak{M}_5 < 0$, and $| \mathfrak{M}_5| > | \mathfrak{M}_4|+| \mathfrak{M}_6|$. Then,

so that

Since $-\bigg(a \Psi_1 \mathsf{I}_v(1+\vartheta \mathsf{I}_v)\frac{ \mathsf{E}_v}{ \mathsf{S}_v}\bigg)| \mathfrak{M}_3| < 0$ and $-\bigg(a \Phi_1 \mathsf{I}_h \mathsf{S}_v(1+\vartheta \mathsf{I}_h)\frac{ \mathsf{E}_v}{ \mathsf{S}_v}\bigg)| \mathfrak{M}_4| < 0$, we then have

Case 14: $U_1 < U_2, \; \; \; \; \; \mathfrak{M}_4, \mathfrak{M}_6 > 0, \mathfrak{M}_5 < 0$, and $| \mathfrak{M}_5| < | \mathfrak{M}_4|+| \mathfrak{M}_6|$. Then,

so that

Since $| \mathfrak{M}_2| < U_1 < U_2 < | \mathfrak{M}_4|+| \mathfrak{M}_6|$, we then have

Case 15: $U_1 < U_2, \; \; \; \; \; \mathfrak{M}_5, \mathfrak{M}_6 > 0, \mathfrak{M}_4 < 0$, and $| \mathfrak{M}_5| < | \mathfrak{M}_4|$. Then,

so that

Since $-\bigg(a \Psi_1 \mathsf{I}_v(1+\vartheta \mathsf{I}_v)\frac{ \mathsf{E}_v}{ \mathsf{S}_v}\bigg)| \mathfrak{M}_2| < 0$, thus we have

Case 16: $U_1 < U_2, \; \; \; \; \; \mathfrak{M}_5, \mathfrak{M}_6 > 0, \mathfrak{M}_4 < 0$ and $| \mathfrak{M}_5| > | \mathfrak{M}_4|$. Then,

so that

Since $| \mathfrak{M}_3| < U_1 < U_2 < | \mathfrak{M}_5|+| \mathfrak{M}_6|$ and $| \mathfrak{M}_4| < | \mathfrak{M}_5|$, it can be deduced that

Now, we examine the subsystem within the framework of system (2.2), where

Now, let us reformulate system (4.38) as

The integrating factors for the system can be expressed as $e^{t(\beta+\mathfrak{U}_h)}$, $e^{t(\mathfrak{U}_h)}$, and $e^{t(\mathfrak{U}_v)}$. We utilize these integrating factors to solve system (4.39). As time $t$ becomes significantly large, i.e., as $t\rightarrow \infty$, $\mathsf{I}_h$ tend to $\mathsf{I}^*_h$, $\mathsf{R}_h$ tends to $R^*_h$, and $\mathsf{I}_v$ approaches $\mathsf{I}^*_v$. This behavior demonstrates that the endemic equilibrium point $E^*$ is indeed globally asymptotically stable.

5.   Bifurcation analysis

Controlling and mitigating epidemic diseases, bifurcation analysis plays a pivotal role. In epidemiological modeling, it is a classical requirement to ensure the disease's eradication by imposing a condition on the basic reproduction number, denoted as $R_0$, which must be less than unity ^[30,31]. Furthermore, researchers have extensively studied the bifurcation analysis of dynamical systems ^{[32,33,34,35,36,37,38,39,40,41]}. However, it is noteworthy that if a bifurcation occurs, the existence of an endemic equilibrium is guaranteed even when $R_0 < 1$. This highlights the significance of conducting bifurcation analysis in the fields of epidemiology and preventive medicine. We demonstrate that system (2.2) exhibits bifurcation by using central manifold theory ^[42]. Specifically, we set $R_0 = 1$ to determine $\Psi_1^*$ as

We calculate the Jacobian matrix at the DFE as

Consequently, the characteristic equation of $J^{0}(\Psi_1^*)$ becomes $(\lambda+\mathfrak{U}_h)(\lambda+\mathfrak{U}_h)(\lambda+\mathfrak{U}_v)(\lambda^4+a_1\lambda^3+a_2\lambda^2+a_3\lambda) = 0,$ where

The existence of a zero eigenvalue within the spectrum serves as a confirmation of the occurrence of a bifurcation phenomenon. Following the approach outlined in ^[42], we proceed to determine the left eigenvectors associated with the matrix $J^{0}(\Psi_1^*)$:

Furthermore, we compute the right eigenvectors in a manner consistent with established methodologies:

When we substitute the following values into Eq (2.2)

● $\mathsf{S}_h$ for $x_1$

● $\mathsf{E}_h$ for $x_2$

● $\mathsf{I}_h$ for $x_3$

● $\mathsf{R}_h$ for $x_4$

● $\mathsf{S}_v$ for $x_5$

● $\mathsf{E}_v$ for $x_6$

● $\mathsf{I}_v$ for $x_7$

we obtain

Considering both the left and right eigenvectors, we can compute coefficients $a$ and $b$ in accordance with the methods described in references ^[43] and ^[24]. Consequently, we have

This relationship suggests that

Next, it is evident that

We also find

By this relation, we have

Hence, it follows that

Given that $a < 0$ and $b > 0$, a forward bifurcation is present in our proposed system, as described in Eq (2.2).

6.   Data fitting and numerical simulations

Case study: Kohat District, Khyber Pakhtunkhawa (Pakistan)

The Kohat district in the KPK province has been significantly affected by Cutaneous Leishmaniasis (CL), largely due to the influx of internally displaced people from other districts. A brief study conducted by Hussain et al. ^[3] not only delineates the prevalence of CL from April 2015 to May 2016 but also identifies the specific species responsible for CL spread among the local population and internally displaced persons (IDPs) in Kohat. The total population of the Kohat district in 2015 was approximately 0.2 million ^[3]. Our focus is on Ghamkol, a location in Kohat with a population of 3256 people.

In this section, we utilized the least squares curve fitting method to analyze reported CL instances in the Kohat district from April 2015 to September 2015. The system's (2.2) estimated parameters are based on Pakistan's overall data on confirmed events and fatalities. The ordinary least squares solution (OLS) is employed to reduce the error terms of daily reports, and the goodness of fit test uses the related relative error.

where the reported total number of infected is \(\mathsf{I}_{i}\), and the simulated total number of infected is \(\hat{ \mathsf{I}}_{i}\). The estimated values of the parameters are provided in Table 1, while a comparison of the model with the reported cases is shown in Figure 2.

In Figure 2, the total cases of Leishmania virus are depicted from April 2015 to September 2015, showcasing the peak period of Leishmania cases in the Kohat district. We conducted a data fitting procedure where we employed our Leishmania model's infected human population class. This analysis clearly demonstrates the suitability of our model in capturing the observed dynamics of the infected human population in real-world data.

7.   Sensitivity analysis

The identification of influential parameters in mitigating the spread of infectious diseases is a crucial aspect addressed through sensitivity analysis in epidemiological modeling. In particular, forward sensitivity analysis holds significant importance, although its application can be computationally demanding, especially in intricate biological models. Within this context, the sensitivity analysis of the basic reproduction number ($R_0$) has garnered substantial attention from the ecological and epidemiological communities. To formally define the concept, we introduce the notion of the normalized forward sensitivity index for $R_0$, which depends differentiably on a parameter denoted as $\omega$.

Definition 1. The normalized forward sensitivity index ($S_\omega$) of $R_0$, with respect to a parameter $\omega$, is defined as

In practice, three distinct methods are conventionally employed to compute these sensitivity indices:

● (ⅰ) Direct Differentiation: Calculating $S_\omega$ directly as $S_\omega = \frac{\omega}{R_0} \frac{\partial R_0}{\partial \omega}$, where $\omega$ represents a parameter.

● (ⅱ) Latin Hypercube Sampling: Employing a Latin hypercube sampling technique.

● (ⅲ) Linearization of the System: Linearizing the underlying system and subsequently solving the resultant set of linear algebraic equations.

In our study, we choose to employ the direct differentiation method because of its ability to provide analytical expressions for these sensitivity indices ^[44,45]. These indices play a dual role by shedding light on the influence of various factors related to the transmission of infectious diseases and supplying valuable data regarding the relative alterations between $R_0$ and various parameters. Consequently, this knowledge assists in the formulation of effective disease control strategies.

An examination of the results presented in Table 2 reveals that certain parameters, specifically $a$, $\Phi_1$, $\mathfrak{U}_h$, $\Gamma_v$, $\Psi_1$, $\mathcal{K}_1$, and $\mathcal{K}_2$, exert a positive influence on the basic reproduction number ($R_0$). This implies that a 10% increase or decrease in these parameters would correspondingly result in a 10 percent increase or decrease in the value of $R_0$. More specifically, these parameters exhibit sensitivity indices of 10%, 4.9%, 1.3%, 5.0%, 5.0%, 1.0%, and 1.75%, respectively. Conversely, parameters such as $\theta$, $\beta$, $\mathfrak{U}_v$, and $\Gamma_h$ are found to have a negative influence on the reproduction number $R_0$. In practical terms, this means that increasing these parameters by 10 percent would lead to a decrease in the value of $R_0$. The sensitivity indices for these parameters reveal the extent of this impact, with values of –4.79%, –2.27%, –0.17% and –5.0%, respectively. These findings provide valuable insights into the critical parameters affecting disease transmission and can inform targeted intervention strategies.

8.   Numerical simulations and dynamic analysis

In this investigation, we have conducted a comprehensive analysis of the global dynamics and sensitivity of an Anthroponotic Cutaneous Leishmania epidemic model. Our primary focus has centered on the basic reproduction number ($\mathcal{R}_0$), a fundamental parameter crucial for determining the overall stability of the proposed model. Key findings are presented in Theorems 4.1 and 4.2. Theorem 4.1 establishes global stability concerning the disease-free equilibrium when the threshold parameter ($\mathcal{R}_0$) remains less than or equal to unity. Conversely, Theorem 4.2 establishes global stability regarding the endemic equilibrium point when the threshold parameter surpasses the value of one. Global dynamics are examined using the Castillo-Chavez method and a geometrical approach, providing a comprehensive understanding of the Anthroponotic Cutaneous Leishmania epidemic model's intricate behavior. To validate our theoretical results, we conducted simulations using the parameters outlined in Table 1. Notably, the calculated value of the basic reproduction number ($\mathcal{R}_0 = 0.000000162315016773084$) falls below the critical threshold of one, affirming the global asymptotic stability of the disease-free equilibrium, as per Theorem 2.1. Additionally, disease compartments $\mathsf{S}_h(t)$, $\mathsf{E}_h(t)$, $\mathsf{I}_h(t)$, $\mathsf{R}_h(t)$, $\mathsf{S}_v(t)$, $\mathsf{E}_v(t)$, and $\mathsf{I}_v(t)$ converge toward their respective disease-free equilibrium points, as visually depicted in Figure 4.

This section is dedicated to elucidating the dynamics predicted by the Anthroponotic Cutaneous Leishmania model (2.2) through graphical representations. The primary objective is to scrutinize the influence of crucial parameters on the dynamics and assess the potential for controlling the epidemic. The model (2.2) is numerically solved using an iterative scheme based on the standard RK-4 method, and the results are presented in Figures 4–6. Figure 4(a) illustrates the diminishing trajectory of susceptible human individuals ($\mathsf{S}_h$) over time, aligning with expectations as infected individuals transition out of the susceptible class. In Figure 4(b), the population of exposed human individuals ($\mathsf{E}_h$) displays a gradual decrease, accelerating as the interaction duration between individuals is shortened. Similarly, the infected individuals ($\mathsf{I}_h$) depicted in Figure 4(c) show a gradual increase in the early weeks, corresponding to the decline in the exposed class. Subsequently, the infected class gradually decreases, approaching the equilibrium point and achieving stability. With the decline in the number of infected individuals, there is a rapid increase in the recovered human population ($\mathsf{R}_h$), as shown in Figure 4(d). The decline in susceptible vectors ($\mathsf{S}_v$) due to viral exposure is demonstrated in Figure 4(e), aligning with epidemiological models' expectations. Figure 4(g) portrays the population of exposed vectors ($\mathsf{E}_v$), exhibiting steady and rapid growth attributed to susceptible vectors becoming infected during the initial weeks of the outbreak. Finally, Figure 4(g) delineates the trajectory of infected vectors ($\mathsf{I}_v$), indicating a substantial portion exiting the exposed class within a few weeks of exposure. This observation underscores the dynamic evolution of the Anthroponotic Cutaneous Leishmania epidemic, providing insights into infection progression in both human and vector populations. The graphical representations illustrating variations in all compartments in response to changes in the parameter $\alpha$, representing the sandfly biting rate, are presented in Figure 5. As the sandfly biting rate increases, there is a discernible upward trend in the counts of individuals and vectors in classes $\mathsf{I}_h$ and $\mathsf{I}_v$, as depicted in Figures 5(c), (g). This observation aligns with expectations, as an elevated sandfly biting rate results in more susceptible individuals being exposed to ACL infections. Additionally, the influence of changes in the natural death rate ($\mu_h$) of individuals demonstrates a continuous increase in the population, consistent with the trends depicted in Figure 6.

9.   Conclusions

In summary, this study has delved into the qualitative behavior of an Anthroponotic Cutaneous Leishmania epidemic model characterized by a convex incidence rate. Through the application of advanced mathematical and epidemiological techniques, we have significantly advanced our comprehension of the dynamics inherent in the transmission of this disease. Our investigation has led to several key insights and potential directions for future research:

● Epidemiological Modeling Significance: This study underscores the crucial role of mathematical modeling in understanding and managing infectious diseases. Drawing parallels with the COVID-19 pandemic highlights the importance of proactive modeling efforts in anticipating and addressing emerging health challenges.

● Stability Analysis Contributions: Meticulous exploration of the local and global stability conditions of the Disease-Free Equilibrium (DFE) and endemic equilibrium states, contingent upon $R_0$, provides critical insights into the potential persistence or eradication of Cutaneous Leishmaniasis.

● Parameter Sensitivity Insights: The sensitivity analysis identifies key parameters influencing the threshold number ($R_0$), a crucial determinant of disease spread. Understanding the sensitivity of these parameters is essential for devising effective control strategies.

● Geometric Approach for Global Stability: The extension of Lyapunov theory using a geometric approach, incorporating a secondary additive compound matrix to establish global stability at the endemic equilibrium, represents an innovative contribution to the toolbox of methods for analyzing disease dynamics.

● Practical Validation Through Simulations: The theoretical findings have been validated through numerical simulations, demonstrating the practical applicability of our results.

To further broaden the impact of this research, future investigations could consider reformulating an epidemic model for visceral leishmaniasis, incorporating a harmonic mean-type incidence rate. This expansion would open avenues for exploring stability analysis, parameter sensitivity, bifurcation analysis, and optimal control strategies, significantly enhancing our ability to combat visceral leishmaniasis and similar infectious diseases.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors are thankful to the Deanship of Scientific Research at Najran University for funding this work under the General Research Funding program grant code (NU/RG/SERC/12/22).

Conflict of interest

The authors declare there is no conflict of interest.