Existence theory and numerical solution of leptospirosis disease model via exponential decay law

1.   Introduction

Leptospirosis is a disease caused by infection with leptospira bacteria. These bacteria can be found worldwide in soil and water. There are many strains of leptospira bacteria that can cause disease. Leptospirosis is a zoonotic disease, which means it can be spread from animals to people. Infection in people can cause flu-like symptoms and can cause liver or kidney disease. In the United States, most cases of human leptospirosis result from recreational activities involving water. Infection resulting from contact with an infected pet is much less common, but it is possible. Leptospirosis is more common in areas with warm climates and high annual rainfall but it can occur anywhere ^[1].

Signs of leptospirosis may include fever, shivering, muscle tenderness, reluctance to move, increased thirst, changes in the frequency or amount of urination, dehydration, vomiting, diarrhea, loss of appetite, lethargy, jaundice (yellowing of the skin and mucous membranes), or painful inflammation within the eyes. The disease can cause kidney failure with or without liver failure ^[2]. Leptospirosis is generally treated with antibiotics and supportive care. When treated early and aggressively, the chances for recovery are good but there is still a risk of permanent residual kidney or liver damage ^[3,4].

There are so many models proposed for the dynamics of humans ^{[5,6,7,8,9,10,11]} as well as vector populations ^[12,13]. To analyze the nature of leptospirosis, Pongsumpon ^[14] develop a mathematical model where they show the rate of change of both humans and rat populations. The population of humans is further split into two major groups: The juveniles and the adults. A dynamical model for leptospirosis transmission was considered by Triampo et al. ^[15]. In their work, they considered a number of leptospiroses infections in Thailand and worked on their numerical simulations. In ^[16], Zaman took the real data and studied the behavior of the theory optimal control.

The generalisation of the classical calculus is called fractional calculus. Mathematical models for the epidemic disease with integer order derivative was recently studies by the researchers ^{[17,18,19,20]}. But mathematical modelings with variable order derivatives give an advance knowledge of a system. Fractional order models are very capable to simulate the memory effects and crossover natures and gives a big degree of performance. Fractional mathematical modelings provides more insights about a epidemic under given circumstances ^{[21,22,23,24]}. Various non-classical type operators with non-singular and singular type kernels have been exemplified in the literature ^[25,26]. The successful applications of such non-classical operators can be seen in early literature and papers therein ^[27,28]. Early, only some leptospirosis models defined by variable order operators are given. The researchers in ^[29] proposed a mathematical structure for cancer and hepatitis co-dynamics in non-classical derivative and studied its outcomes.

Nowadays, the mathematical models involving fractional order derivative were given noticeable importance because they are more accurate and realistic as compared to the classical order models ^{[30,31,32,33,34,35]}. "Motivated by the advancement of fractional calculus, many researchers have focused to investigate the solutions of nonlinear differential equations with the fractional operator by developing quite a few analytical or numerical techniques to find approximate solutions ^{[36,37,38,39,40]}. These differential equations involves several fractional differential operators like Riemann-Liouville, Caputo, Hilfer etc." ^[41,42,43].

However, "these operators possess a power law kernel and have limitations in modeling physical problems. To overcome this difficulty, recently an alternate fractional differential operator having a kernel with exponential decay has been introduced by Caputo and Fabrizio ^[44]. This novel approach of fractional derivative is known as the Caputo-Fabrizio (C-F) operator which has attracted many research scholars due to the fact that it has a non-singular kernel. Also the C-F operator is most appropriate for modeling some class of real-world problem which follows the exponential decay law" ^[44]. With the passage of time, developing a mathematical model using the C-F fractional order derivative became a remarkable field of research. In recent times, several mathematicians were busy in development and simulation of CFFDE.

We develop a vector-host epidemic model in which the total population of the host (human) at time $(t)$ is divided by susceptible $S_h$, infected $I_h$ and recovered $R_h$ humans. The population of vectors at time $(t)$ is classified by susceptible $S_v$ and infected $I_v$ classes. The total amount of host is exemplified by $N_h$ and the total amount of vectors are denoted by $N_v$. Here, $N_h = S_h+I_h+R_h$ and $N_v = S_v+I_v$. The mathematical representation of the model is given by

Here, $b_1$ is the rate of recruitment in humans size, susceptible humans can get infection by two roots, one is direct transmission or via infectious individuals where $\beta_{1}$ and $\beta_{2}$ are the mediate transmission rates, $\mu_{h}$ is the natural mortality rate for hosts where $\gamma_{h}$ is the rate of recovery of humans from the infection, $\delta_{h}$ disease death rate of infected humans due to disease, $\lambda_{h}$ is the constant rate of humans immune, $b_{2}$ recruitment rate for human density, $\delta_{v}$ is the mortality rate of vectors due to disease, $\beta_{3}$ represent disease carrying due to susceptible vectors per host per unit time and $\gamma_{v}$ is the natural mortality rate of vector.

Motivated by the above discussion, we modify the above model (1.1) by using Caputo-Fabrizio fractional operator, and the modified model is given as follow:

To exemplify the existence and uniqueness of equilibrium solutions of the given system, the number of novel results are utilized with the help of Sumudu transform. By equilibrium solutions existence, we are concluding the well-posedness of the given system.

Section 2 is devoted to some basic definitions regarding fractional operators. In section 3, the uniqueness of the solution is derived. Stability analysis of the given algorithm by the application of fixed point theorem is carried out in Section 4. Section 5 deals with the derivation of special solutions and numerical results for the model. The sensitivity of different parameters corresponding to basic reproduction is discussed in Section 6. Section 7 deals with the conclusion and future work.

2.   Preliminaries

Here we recall some definitions alongwith Caputo-Fabrizio operator of non-classical order which is defined with exponentially decay type kernel ^[44,45].

Definition 2.1. Assume $\phi\in G^{1}(a, b)$, $b > a$, $\alpha\in (0, 1)$, then the definition of Caputo-Fabrizio derivative of variable order is given by

where $N(\alpha)$ is the normalization function with $N(0) = N(1) = 1$. Also, when the function does not exist in $G^{1}(a, b)$, then the given operator can be defined by

Remark 2.1. If we fix $\epsilon = \frac{1-\alpha}{\alpha}\in [0, \infty)$, $\alpha = \frac{1}{1+\alpha}\in [0, 1]$, then the definition of given Caputo-Fabrizio variable order derivative is redefined by

In addition,

The related fractional Caputo-Fabrizio integral of the given derivative is simulated by Nieto and Losada ^[46] which is as follow:

Definition 2.2. For $\alpha\in (0, 1)$, the integral of fractional order $\alpha$ of the function $f$ is expressed by

Remark 2.2. It is to be noted that, according to the above definition, the function of order $0 < r\leq1$, which is fractional integral of Caputo type is just average among mapping $\psi$ and the integral of it of order $1$. Therefore, that can be written as

The above expression will assumes the form as

Nieto and Losada ^[46] rewrote that this Caputo-Fabrizio derivative of order $0 < \alpha < 1$ can be reformulated as

Definition 2.3. The Sumudu transform of a mapping $\phi(t)$ defined on a set specified by

is given by

Definition 2.4. For the Caputo-type fractional operator of any function $\phi(t)$, the Sumudu transform is defined by

Similarly, the Sumudu transform for any function $\phi(t)$ in the sense of Caputo-Fabrizio operator is given as follow:

Definition 2.5. For any function $\phi(t)$ where the C-F derivative is defined, the Sumudu transform in a sense of C-F is expressed by

3.   Special solution derivation

In this section, we apply Sumudu transform to both sides of the model (1.2) using Picard's iterative mathod ^[22]:

Rearranging, we get

By operating every side of Eq (3.2) by the inverse of Sumudu transform, we obtain

The recursive formula for given model can be exemplified by

The solution of (3.4) is expressed as

For simulating the given model (3.4), we use following values $S_h = 1.2$, $I_h = 2.8$, $R_h = 1.3$, $S_v = 2.8$, $I_v = 1.15$, $b_1 = 0.121$, $\mu_h = 0.0121$, $\beta_1 = 0.14$, $\beta_2 = 0.02$, $\beta_3 = 0.10$, $\lambda_h = 1$, $\delta_h = 0.8$, $\gamma_h = 0.25$, $b_2 = 0.002$, $\gamma_v = 0.23$, $\delta_v = 0.001$. Then

4.   Stability analysis of given algorithm by the application of fixed-point theorem

Suppose a Banach space $(X, \lVert\cdot\lVert)$ and a self-map $K^{*}$ of $X$. Now consider a random recursive algorithm of the formate, $x_{n+1} = g(K^{*}, x_{n})$. Assume $F(K^{*})$ be the set of fixed-point of $K^{*}$ which is non-empty and a point $k^{*}\in F(K^{*})$ where $x_{n}$ converges. Adopt $\{{y^{*}_{n}\subseteq}X\}$ and derive $e_{n} = \lVert y^{*}_{n+1}-g(K^{*}, y^{*}_{n})\lVert$. The numerical algorithm $x_{n+1} = g(K^{*}, x_{n})$ is called $K^{*}$-stable if $\lim_{n\rightarrow \infty}e^{n} = 0$ implies that $lim_{n\rightarrow \infty}y_{n}^{*} = k^{*}$. Analogously, we consider that the sequence $\{y^{*}_{n}\}$ exists an upper bound, else the sequence will be not convergent. If all such constraints are satisfied for $x_{n+1} = K^{*}x_{n}$ which is known as Picard's iteration as given in ^[47], consequently the iteration will be $K^{*}$-stable. We then state the following theorem.

Theorem 4.1. Consider a Banach space $(X, \lVert\cdot\lVert)$ and self-map $K^{*}$ of $X$ agreeing

$\forall x, y\in X$, here, $0\leq C, 0\leq c < 1$. Assume that $K^{*}$ is Picard $K^{*}$-stable. Now consider the recursive algorithm from (3.4) related to (1.2).

where $\frac{(1-\alpha+\alpha u)}{N(\alpha)}$ is the non-classical Lagrange multiplier.

Theorem 4.2. Consider the self-map $H$ given by

is $H$-stable in $L^{1}(a, b)$ if

Proof. Initially we start to proving that $H$ exists a fixed-point. For this achievement, we simulate the given below terms for all $(m, n)\in M\times M$.

Analysing the initial equation of (4.2) and take norm on both-sides,

Applying the triangular inequality, we receive

By further simplification, (4.4) becomes as

Since the same role is forming by both solutions, so we adopt in that case

By fixing it in Eq (4.5), we will get the following result:

Since $S_{hm}$, $S_{vn}$, $I_{hm}$, $I_{hn}$, $I_{vn}$ are convergent sequence, so all are bounded, then, there exist five different positive constants: $J$, $Q$, $P_{1}$, $P$ and $J_{1}$ for all $t$ such that

Now by Eqs (4.6) and (4.7), we get

where $f$ and $g$ are functions from $\mathcal{S}^{-1} \bigg[\mathcal{S}(\frac{1-\alpha+\alpha u}{N(\alpha)})\bigg].$

Similarly, we receive

where

which gives that the non-linear $H$-self map exists a fixed point. Now we prove that $H$ agrees with all the restrictions of Theorem 4.1. Let (4.8) and (4.9) hold, so using

The above expression shows that all conditions of Theorem 4.1 exist for the nonlinear map-$H$. So every restrictions mentioned in Theorem 4.1 are agreed for the given non-linear map-$H$. Hence $H$ is Picard $H$-stable, which finishes the proof of Theorem 4.2.

5.   Special solution uniqueness

In this portion, we are going to show the uniqueness of the special solution of Eq (1.2) applying the numerical algorithm. Here first we assume an exact solution for Eq (1.2), for which the special solution converges for a large quantity $m$. Here we take the Hilbert space $H = L((a, b)\times (0, T))$ which can be formulated as follow:

Now, we adopt the given operator

The target of this portion is to prove that the inner product of

where $(Z_{11}-Z_{12})$, $(Z_{21}-Z_{22})$, $(Z_{31}-Z_{32})$, $(Z_{41}-Z_{42})$ and $(Z_{51}-Z_{52})$ are special solution of system. However,

Here, we calculate the initial equation showing in the model without forgetting the general procedure:

Because the play of both solutions is nearly same, so we can consider that

Then Eq (5.2) becomes

By the relation of inner product and norm, we receive the given results:

Again repeating the similar manner, from the second equation of the model (5.1), we get the outputs

Following similar manner, the third equation of the model (5.1) gives that

The fourth one is

And the fifth one is

Upon putting Eqs (5.3)–(5.7) in Eq (5.1), we get

But, for sufficiently bigger value of $m_i,$ where $i = 1, 2, 3, 4, 5$ both solutions converge to the exact solution, applying the consequences of topological results, here three very tiny positive parameters exist $l_{m_{1}}, l_{m_{2}}, l_{m_{3}}, l_{m_{4}}$and $l_{m_{5}}$ s.t.

Now putting the exact solution to the right-side of Eq (5.8) and using the triangular inequality with taking $M = \max(m_{1}, m_{2}, m_{3}, m_{4}, m_{5}), \; l = \max(l_{m_{5}}, l_{m_{4}}, l_{m_{3}}, l_{m_{2}}, l_{m_{1}}),$ we get

As $l$ is very very small positive value, so on the basis of topological result, we get

But, it is obvious that

Therefore, we have

which yields that

This completes the proof of uniqueness.

6.   Numerical simulations and discussions

This portion gives the numerical simulation results based on the numerical algorithm as Eq (3.5) for leptospirosis disease model (1.2). The numerical method employed on Eq (1.2) hinged with (3.5). The numerical simulations are performed by adopting the parameter numerical values, $b_1 = 0.121$, $\mu_h = 0.0121$, $\beta_1 = 0.14$, $\beta_2 = 0.02$, $\beta_3 = 0.10$, $\lambda_h = 1$, $\delta_h = 0.8$, $\gamma_h = 0.25$, $b_2 = 0.002$, $\gamma_v = 0.23$, $\delta_v = 0.001$. In this practical work, the integer order cases are plotted in the company of Figure 1. Here we have just checked that how the given model is behaving at the classical sense. After following it, in the assembly of Figure 2, we plotted the model dynamics at distinct derivative operator orders $\alpha.$ Here the disease shows almost same nature at each values but just the amount of the given population varies at distinct values. The assembly of Figures 3 and 4 is dedicated to the model dynamics at the values of parameters $\beta_1, \; \beta_2, \; \beta_3$. From the given plots we can see that the impact of fractional order is measurable and every value demonstrates a unique nature of the given disease. The given parameter range is the best fit for the given model structure.

7.   Sensitivity analysis

Here sensitivity analysis is performed to show that which parameters are useful in reducing the infection spread of the disease. Forward sensitivity structure is taken as a important part of epidemic modelling, albeit its analysis become austere for tough biological dynamics. Epidemiologist and ecologist have given their much attention to the calculation of sensitivity analysis of $R_0$.

Definition 7.1. The normalized forward sensitivity index of the $R_0$ which based on differentiably on a variable $\Gamma$ is given by

There are three general techniques which are basically utilized to find the sensitivity indices: (ⅰ) By a Latin hypercube sampling technique; (ⅱ) by direct differentiation; (ⅲ) by linearising model and then simulating the obtained set of linear algebraic constraints. Here we adopt the direct differentiation technique because it gives analytical forms for the indices. These indices demonstrates the impact of different supports jointed with the infection spread and also provides us considerable knowledge about the comparative variations between $R_0$ and various parameters. Definitely, it provides the help in texturing the control techniques.

The parameter values $b_1$, $\beta_2$, $\beta_1$, $\lambda_h$, $\delta_h$, $\gamma_h$, $b_2$ and $\beta_3$ have a positive impact on the reproductive number $R_0$, which exemplify that the decay or growth of given parameters call by $10\%$ will decrease or increase the reproductive number by $10\%$, $1.70\%$, $8.2\%$, $4.27\%$, $2.3\%$, $1.7\%$, $1.7\%$ and $0.2\%$ respectively. But another side, the parameters index of $\mu_h$, $\gamma_v$ and $\delta_v$ shows that raising their values by 10% will reduce the value of reproductive number $R_0$ by 9.9%, 2.6% and 1.0% respectively. $\lambda_h$ has no impact on reproduction number. Graphical representations of the sensitivity index is explored in Figures 5–11 which demonstrate the sensitivity of different parameters.

8.   Conclusions

In this research, we have simulated a leptospirosis model dynamics by using a non-classical type derivative with the application of Banach contraction theorem and Picard successive approximation algorithm. Here the fractional order derivative called Caputo-Febrizo of order $\alpha$ has been used which exists a non-singular exponentially decay type kernel. The given analysis have been justified by applying well-known Sumudu transform alongwith proposed iterative scheme. Using the Banach theorem, the existence of unique equilibrium solution have been exemplified. With the help of Matlab, we have also presented numerical simulations to the approximate solutions which show the effectiveness of the theoretical results. From the given simulations, we concluded that the mentioned non-classical type model gives more emphatic structure as compare to the integer-order dynamics. Hence we can say that given variable order structure approaches to time responses with super-slow evolutions and super-fast transients directed to the steady-state, which impacts can not be smoothly textured by the ordinary order systems. In the future, the above model texture can be further rebuilt by more effective and advance models of dynamical systems. One may use Atangana Baleanu fractional derivative or Hilfer fractional derivative to compare the results with our model. One can also employ the fractional optimal control strategy to decrease the infected people and increased the susceptible people by choosing suitable control variables. The stochastic version of this model will also be useful to show the effect of white noise on the model.

Acknowledgments

This research was supported by King Mongkut's University of Technology Thonburi's Postdotoral Fellowship.

Conflict of interest

The authors declare no conflicts of interest.