1. Introduction
Mathematical modeling serves as an indispensable tool for comprehending the dynamic mechanisms underlying the spread of infectious diseases and devising effective strategies to mitigate economic losses. A pivotal achievement in mathematical epidemiology lies in the formulation, prediction, and identification of disease control measures through rigorous mathematical expressions. Traditionally, epidemiological models have compartmentalized a constant population into three categories: susceptible, infectious, and recovered. The evolution of these models has seen the inclusion of various compartment classes, adapting to the variation of available descriptive information. In the budding stages of research, the discourse predominantly revolved around integer-order differential equations. However, recent strides in the field have witnessed the emergence of fractional-order mathematical modeling, playing a substantial role in addressing a diverse array of challenges in engineering, epidemiology, ecology, and biology. Fractional calculus, as a generalization where the conventional integer order gives way to a fractional order, provides a more varied approach
[21] | S. B. Kurniawan, “Bioprocess of Trivalent Chromium using Bacillus subtilis and Pseudomonas putida”, Department of Environmental Eng., FTSLK ITS Surabaya, 2016. |
[22] | Sherene T., “Heavy Metal status of soils in Industrial Belts of Coimbatore District, Tamil Nadu”, Nature environment and pollution technology, Vol. 8, No. 3, pp. 613-618, 2009. |
[23] | S. S Askar, Dipankar Ghosh, “A fractional order SITR Mathematical model for forecasting of transmission of Covid 19 of India with lockdown effect”, Results in Physics 24, 104067, 2021. |
[24] | Subrata Paul and Animesh Mahata, “Dynamics of SIQR epidemic model with fractional order derivative”, Partial differential equation in applied mathematics, vol. 5, 100216, 2022. |
[25] | V. Lakshmikantham and A. S. Vatsala, “Theory of fractional differential inequalities and applications,” Communications in Applied Analysis, vol. 11, pp. 395–402, 2007. |
[21-25]
. Notably, as the order approaches one, the solution of the fractional order converges to that of the integer order system. The fractional order system introduces two critical properties—the memory property and hereditary property—which are not adequately captured by any integer order system. These properties contribute to a more accurate representation of real-world phenomena. The fractional order system aligns well with theoretical descriptions of various mechanisms, offering a sophisticated understanding of dynamics
[26] | V. Lakshmikantham and A. S. Vatsala, “Basic theory of fractional differential equations,” Nonlinear Analysis: Theory, Methods & Applications, vol. 69, no. 8, pp. 2677–2682, 2008. |
[28] | Wu, Jianyong, et al. "Sensitivity analysis of infectious disease models: methods, advances and their application." Journal of The Royal Society Interface 10.86(2013): 20121018. |
[29] | Y. Li, Y. Chen and I. Podlubny, “Stability of fractional-order nonlinear dynamic systems: Lyapunov direct method and generalized Mittag-Leffler stability,” Computers & Mathematics with Applications, vol. 59, no. 5, pp. 1810–1821, 2010. |
[26, 28-29]
. Over the last few decades, various types of fractional calculus, including Riemann–Liouville, Caputo, Atangana–Baleanu, Caputo–Fabrizio, Katugampola, Hadamard, etc., have been introduced extensively to scrutinize the dynamics of epidemic models. This ongoing exploration holds the promise of refining our capacity to model and comprehend intricate systems across diverse fields.
[1] | A. Kilbas, H. M. Srivastava and J. J. Trujillo, “Theory and Applications of Fractional Differential Equations”, Elsevier, Amsterdam, The Netherlands, 2006. |
[2] | Babakhani and V. Daftardar-Gejji, “Existence of positive solutions of nonlinear fractional differential equations”, Journal of mathematical analysis and applications, vol. 278, no. 2, pp. 434-442, 2003. |
[3] | C. F. Li, X, N. Luo, Y. Zhou and J. J, “Existence of positive solution of the boundary value problem for nonlinear fractional differential equations”, Computer and Mathematics with applications, vol. 59, pp. 1363-1375, 2010. |
[4] | D. Kumar and J. Singh, “Application of Homotopy analysis transform method to fractional biological population model”, Romanian reports physics, vol. 65, no. 1, pp. 63-75, 2013. |
[5] | Elhia, M., Balatif, O., Boujallal, L., & Rachik, M. Optimal control problem for a tuberculosis model with multiple infectious compartments and time delays. An International Journal of Optimization and Control: Theories & Applications (IJOCTA), 11(1), 75-91. (2021). |
[6] | Evirgen, Fırat, et al. "Real data-based optimal control strategies for assessing the impact of the Omicron variant on heart attacks." AIMS Bioengineering 10.3 (2023). |
[1-6]
Subtrapaul, et.,
[24] | Subrata Paul and Animesh Mahata, “Dynamics of SIQR epidemic model with fractional order derivative”, Partial differential equation in applied mathematics, vol. 5, 100216, 2022. |
[24]
proposed the dynamical behavior of the COVID-19 epidemic SIQR model in terms of fractional order. S. S. Askar, and Dipankar Ghosh,
[23] | S. S Askar, Dipankar Ghosh, “A fractional order SITR Mathematical model for forecasting of transmission of Covid 19 of India with lockdown effect”, Results in Physics 24, 104067, 2021. |
[23]
studied the h-stability of the Caputo fractional order differential equation by using some fractional comparison principle and the Lyapunov direct method. Mukherjee D, and Mondal R,
[15] | Mukharjee D and Mondal R, “Dynamical analysis of a Fractional Order Prey-Predator system with reserved area”, J Fract Calc Appl 11, 54-69, 2020. |
[15]
deal with the fractional order prey-predator model of the dynamical system for the reserved area and present the numerical analysis of hopf-bifurcation. Zizhen Zhang et al.,
[30] | Zhang, Zizhen, et al. "Dynamics of a fractional order mathematical model for COVID-19 epidemic." Advances in Difference Equations 2020. 1(2020): 1-16. |
[30]
formulated and analyzed the dynamics of the COVID-19 epidemic models with an isolated class of fractional order. Many researchers proposed many models in eco-epidemiology and analyzed the various stability properties, existence, and non-negativity of the solution
[7] | F. Haq, K. Shah, G. Rahman and M, Shahzad, “Existence and uniqueness of positive solutions to boundary value problem of fractional differential equations”, Sindh Univ. Res. Jour, vol. 48, no. 2, pp. 451-456, 2016. |
[8] | Evirgen, F., Ucar, E., ÖZDEMIR, N., Altun, E., & Abdeljawad, T. (2023). The impact of nonsingular memory on the mathematical model of Hepatitis C virus. Fractals, 2340065. |
[9] | F. Haq, K. Shah, G.-U. Rahman and M. Shahzad, “Numerical analysis of fractional order model of HIV-1 infection of CD4+ T cells”, Computational methods for differential Equations, vol. 5, no. 1, pp. 1-11, 2017. |
[10] | Fazal Haq and Muhammad Shahzad, “Numerical Analysis of Fractional Order Epidemic Model of Childhood Disease”, Discrete Dynamics in Nature and Society, vol. 2017, Article Id 4057089, pp. 7, 2017. |
[11] | Hamou, Abdelouahed Alla, et al. "Analysis and dynamics of a mathematical model to predict unreported cases of COVID-19 epidemic in Morocco." Computational and Applied Mathematics 41.6 (2022): 289. |
[12] | Ipung Fitri Purwanti and Tesya Paramita Putri, “Treatment of Chromium Contaminated Soil using Bioremediation”, AIP Conference Proceeding, 1903, 040008, 2017. |
[13] | Kashyap, A. J., Bhattacharjee, D., & Sarmah, H. K. (2021). A fractional model in exploring the role of fear in mass mortality of pelicans in the Salton Sea. An International Journal of Optimization and Control: Theories & Applications (IJOCTA), 11(3), 28–51. https://doi.org/10.11121/ijocta.2021.1123 |
[14] | Liao, Teh-lu. "Adaptive synchronization of two Lorenz systems." Chaos, Solitons & Fractals 9.9 (1998): 1555-1561. |
[15] | Mukharjee D and Mondal R, “Dynamical analysis of a Fractional Order Prey-Predator system with reserved area”, J Fract Calc Appl 11, 54-69, 2020. |
[7-15]
. Numerical approximation of the proposed model has been done by the Rung-Kutta method, Taylor series method, Adams Bashforth and Adams Moulton method, Adomain and Laplace Adomian decomposition method, Homotopy method and generalization of homotopy method
[14] | Liao, Teh-lu. "Adaptive synchronization of two Lorenz systems." Chaos, Solitons & Fractals 9.9 (1998): 1555-1561. |
[14]
. In this work, the Laplace Adomain decomposition method is used to study the approximate solution of the proposed model. The Laplace Adomian decomposition method was introduced by Adomian in 1980. This is a very effective method to find the explicit and numerical solution to the physical problem of ordinary and partial differential equations with nonlinear initial and boundary conditions. It can also be used to study many deterministic, stochastic diffusion models with and without delay. There is no perturbation and linearization is required and also no extra memory is required to solve this problem which is the main advantage of the method
[16] | Muhammad Rafiq and Javaid Ali,“Numerical Analysis of a bi-model covid-19 SITR model”, Alexandria Engineering Journal, vol. 61, pp. 227-235, 2022. |
[17] | Noor Badshah and Haji Akbar, “Stability analysis of fractional order SEIR model for malaria disease in Khyber Pakhtunkhwa”, Demonstratio Mathematica, vol. 54, pp. 326-334, 2021. |
[18] | Nuri Ozalp and Elif Demiric, “A fractional order SEIR model with vertical transmission”, Mathematics and computer modeling, 54, pp. 1-6, 2011. |
[19] | O. D. Makinde, “Adomain decomposition approach to a SIR epidemic model with constant vaccination strategy”, Applied Mathematics and Computation, vol. 184, no. 2, pp. 842-848, 2007. |
[20] | Ozdemir, Necati, Esmehan Uçar, and Derya Avcı. "Dynamic analysis of a fractional svir system modeling an infectious disease." Facta Universitatis, Series: Mathematics and Informatics (2022): 605-619. |
[16-20]
.
Sensitivity analysis (SA) is increasingly popular in developing more efficient models, requiring reliable statistical and mathematical methods to enhance modeling accuracy. SA serves various purposes, including developing decisions or recommendations, communicating, understanding, or quantifying systems, and refining models. During model development, SA verifies validity or accuracy, simplifies and calibrates weak processes or missing data, and identifies important parameters for further studies. These applications highlight SA's versatility and importance in refining models and optimizing decision-making processes.
The main aim of this paper is to study the mathematical model of Caputo fractional order soil pollution. Soil is the main requirement for living and non-living organisms. In the case of contamination of soil as such oil spill, excess use of heavy metals, and inorganic fertilizer, the soils are polluted. All polluted soil contains compounds that include metals, inorganic iron salts (phosphates, carbonates, sulfates, nitrates), and organic compounds (lipids, DNA, fatty acids, alcohol). These compounds are mainly entered through soil microbial activity and decomposition of organisms (animals and plants)
[10] | Fazal Haq and Muhammad Shahzad, “Numerical Analysis of Fractional Order Epidemic Model of Childhood Disease”, Discrete Dynamics in Nature and Society, vol. 2017, Article Id 4057089, pp. 7, 2017. |
[17] | Noor Badshah and Haji Akbar, “Stability analysis of fractional order SEIR model for malaria disease in Khyber Pakhtunkhwa”, Demonstratio Mathematica, vol. 54, pp. 326-334, 2021. |
[18] | Nuri Ozalp and Elif Demiric, “A fractional order SEIR model with vertical transmission”, Mathematics and computer modeling, 54, pp. 1-6, 2011. |
[10, 17, 18]
. The main soil pollution has an unfavorable two-way impact on food security, due to toxic contamination, which can reduce crop yields and be unsafe for consumption by humans and animals. Currently, the degradation of soil and land is affecting 3.2 million people which is 40% of the world’s population. It is necessary to treat the damaged soil and to maintain better soil management practices to limit soil pollution.
In this study, the work is carried out as follows: In Section 2, the mathematical model of the Caputo fractal fractional order SPTR model is presented. Section 3, analyses the existence, uniqueness, boundedness, and non-negative of the solution. In Section 4, we discuss the equilibrium and stability analysis. An adaptive control and sensitivity analysis are performed in Section 5. In Section 6, the numerical approximation of the solution of the model is done. In Section 7, the numerical simulation and discussion are given. The conclusion is given in Section 8.
2. Description of Caputo Fractional Order SPTR Model
The dynamic model for soil pollution comprises four distinct compartments: susceptible soil areas S(t), polluted soil areas (P), areas subject to necessary remedial treatments (T), and the recovered subclasses (R). Soil contamination typically arises from the presence of excessive heavy metals and plastics. Consequently, the application of both organic and inorganic remedies to the soil treatment (T) class becomes crucial, where certain areas may experience recovery while others might reintroduce contaminants to the susceptible soil due to treatment failures. In the contemporary landscape, there is a notable surge in the utilization of heavy metals in daily activities, leading to a recurrent cycle where previously recovered soil reverts to a susceptible state. In light of the above, the transmission dynamics of soil pollution under treatment measures are given in the compartmental representation illustrated in
Figure 1. This model encapsulates the intricate interplay between susceptible and treated soil areas, offering insights into the complex dynamics of soil pollution and the efficacy of remedial measures.
Figure 1. The schematic diagram for the soil pollution model.
The mathematical model of fractional order is presented as follows,
${}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\zeta \mathrm{},\mathrm{\wp}}S\left(t\right)=B-\frac{\mathit{\alpha SP}}{N}+\mathit{\beta T}+\mathit{\delta R}-\mathit{\mu S}$
${}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\zeta ,\wp}P\left(t\right)=\frac{\mathit{\alpha SP}}{N}-\mathit{\gamma P}-{\alpha}_{1}P-\mathit{\mu P}$
${}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\zeta \mathrm{},\wp}T\left(t\right)={\alpha}_{1}P-\mathit{\beta T}-{\gamma}_{1}T-\mathit{\mu T}$
${}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\zeta \mathrm{},\wp}R\left(t\right)=\mathit{\gamma P}+{\gamma}_{1}T-\mathit{\delta R}-\mathit{\mu R}$(1)
with the initial conditions,
$S\left(0\right)={S}_{0}\ge 0,\mathrm{}P\left(0\right)={P}_{0}\ge 0,\mathrm{}T\left(0\right)={T}_{0}\ge 0,\mathrm{}R\left(0\right)={R}_{0}\ge 0.$(2)
where the parameter and their descriptions are as follows, $N$ is the total soil area with $N=S+P+T+R$
B – recruitment rate of the soil,
$\alpha $ - the pollution rate of the soil,
$\beta $- recovered re-enter to the susceptible rate due to failure of the treatment,
$\gamma $ – natural recovered rate,
${\alpha}_{1}$ - remediation rate,
${\gamma}_{1}$- recovered rate due to treatment,
$\delta $- recovered re-enter to the susceptible class,
$\mu $- death rate of all the classes.
3. Existence and Boundedness
For the fractional order system, the existence and uniqueness of solutions in the region ${R}_{M}\times \left(0,t\right)$, where ${R}_{M}=\left\{\left(S,P,T,R\right)\in {R}_{+}^{4}:\mathrm{max}\left(\left|S\right|,\left|P\right|,\left|T\right|,\left|R\right|\le \eta \right)\right\}$ is considered, where $\eta $ is some positive constant.
Theorem 3.1 If
${\mathit{Y}}_{0}=\left({S}_{0},{P}_{0},{T}_{0},{R}_{0}\right)\in {R}_{M},$then for each
${Y}_{0}$ there exists a unique solution
$Y\left(t\right)\in {R}_{M}$ of the system (
1) with initial condition
${Y}_{0}$ for all
$t\ge 0$.
Proof: Let us consider the approach of
[9] | F. Haq, K. Shah, G.-U. Rahman and M. Shahzad, “Numerical analysis of fractional order model of HIV-1 infection of CD4+ T cells”, Computational methods for differential Equations, vol. 5, no. 1, pp. 1-11, 2017. |
[9]
and take
$R\left(Y\right)=({R}_{1}\left(Y\right),{R}_{2}\left(Y\right),{R}_{3}\left(Y\right),{R}_{4}(Y\left)\right)$,
${R}_{1}\left(Y\right)=B-\mathit{\alpha SP}+\mathit{\beta T}+\mathit{\delta R}-\mathit{\mu S}$
${R}_{2}\left(Y\right)=\mathit{\alpha SP}-\mathit{\gamma P}{-\alpha}_{1}P-\mathit{\mu P}$
${R}_{3}\left(Y\right)={\alpha}_{1}P-\mathit{\beta T}-{\gamma}_{1}T-\mathit{\mu T}$
${R}_{4}\left(Y\right)=\mathit{\gamma P}+{\gamma}_{1}T-\mathit{\delta R}-\mathit{\mu R}$(3)
The system of equations (
3) is reduced to,
${R}_{1}\left(Y\right)=B-\mathit{\alpha SP}+\mathit{\beta T}+\mathit{\delta R}-\mathit{\mu S}$
${R}_{2}\left(Y\right)=\mathit{\alpha SP}-\mathit{AP}$
${R}_{3}\left(Y\right)={\alpha}_{1}P-\mathit{BT}$
${\mathrm{}R}_{4}\left(Y\right)=\mathit{\gamma P}+{\gamma}_{1}T-\mathit{CR},\mathrm{}$(4)
where, $A=\gamma {+\alpha}_{1}+\mu $, $\mathit{B}=\beta +{\gamma}_{1}+\mu $, $C=\delta +\mu $. For any$,\stackrel{\u030c}{Y}\in {R}_{M}$,
$\u27e6R\left(Y\right)-R\left(\stackrel{\u030c}{Y}\right)\u27e7=\left|{R}_{1}\left(Y\right)-{R}_{1}\left(\stackrel{\u030c}{Y}\right)\right|+\left|{R}_{2}\left(Y\right)-{R}_{2}\left(\stackrel{\u030c}{Y}\right)\right|+\left|{R}_{3}\left(Y\right)-{R}_{3}\left(\stackrel{\u030c}{Y}\right)\right|+\left|{R}_{4}\left(Y\right)-{R}_{4}\left(\stackrel{\u030c}{Y}\right)\right|$
$=\left|B-\mathit{\alpha SP}+\mathit{\beta T}+\mathit{\delta R}-\mathit{\mu S}-B+\alpha \stackrel{\u030c}{S}\stackrel{\u030c}{P}+\beta \stackrel{\u030c}{T}+\delta \stackrel{\u030c}{R}-\mu \stackrel{\u030c}{S}\right|+\left|\mathrm{}\mathit{\alpha SP}-\mathit{AP}-\alpha \stackrel{\u030c}{S}\stackrel{\u030c}{P}+A\stackrel{\u030c}{P}\right|$
$+\left|{\alpha}_{1}P-\mathit{BT}-{\alpha}_{1}\stackrel{\u030c}{P}+B\stackrel{\u030c}{T}\right|+|\mathit{\gamma P}+{\gamma}_{1}T-\mathit{CR}-\gamma \stackrel{\u030c}{P}-{\gamma}_{1}\stackrel{\u030c}{T}+C\stackrel{\u030c}{R}|$
$\le \mathit{\alpha \eta}\left|S-\stackrel{\u030c}{S}\right|+\mathit{\alpha \eta}\left|P-\stackrel{\u030c}{P}\right|+\beta \left|T-\stackrel{\u030c}{T}\right|+\delta \left|R-\stackrel{\u030c}{R}\right|+\mu \left|S-\stackrel{\u030c}{S}\right|+\mathit{\alpha \eta}\left|S-\stackrel{\u030c}{S}\right|+\mathit{\alpha \eta}\left|P-\stackrel{\u030c}{P}\right|$
$+A\left|S-\stackrel{\u030c}{S}\right|+{\alpha}_{1}\left|P-\stackrel{\u030c}{P}\right|+B\left|T-\stackrel{\u030c}{T}\right|+\gamma \left|P-\stackrel{\u030c}{P}\right|+{\gamma}_{1}\left|T-\stackrel{\u030c}{T}\right|+C\left|R-\stackrel{\u030c}{R}\right|$
$\le \left(2\mathit{\alpha \eta}+\mu \right)\left|S-\stackrel{\u030c}{S}\right|+\left(2\left(\mathit{\alpha \eta}+{\alpha}_{1}+\gamma \right)+\mu \right)\left|P-\stackrel{\u030c}{P}\right|+\left(2\left(\beta +{\gamma}_{1}\right)+\mu \right)\left|T-\stackrel{\u030c}{T}\right|+\left(2\delta +\mu \right)\left|R-\stackrel{\u030c}{R}\right|$
$\le H\mathrm{}\Vert \left(S,P,T,R\right)-\left(\stackrel{\u030c}{S},\stackrel{\u030c}{P},\stackrel{\u030c}{T},\stackrel{\u030c}{R}\right)\Vert $
$\le H\Vert Y-\stackrel{\u030c}{Y}\Vert ,$
where $H=\mathrm{max}\{\left(2\mathit{\alpha \eta}+\mu \right),\left(2\left(\mathit{\alpha \eta}+{\alpha}_{1}+\gamma \right)+\mu \right),\mathit{}2\left(\beta +{\gamma}_{1}\right)+\mu ,\mathit{}2\delta +\mu \}$
Hence
$R\left(Y\right)$ satisfies the Lipshitz condition and so the existence and uniqueness of the solution of fractional order system (
1) with the initial conditions are confirmed.
Theorem 3.2 The system of the fractional order of (
1) in
${R}_{M}$ is non-negative and uniformly bounded.
Proof: The approach used by
[9] | F. Haq, K. Shah, G.-U. Rahman and M. Shahzad, “Numerical analysis of fractional order model of HIV-1 infection of CD4+ T cells”, Computational methods for differential Equations, vol. 5, no. 1, pp. 1-11, 2017. |
[9]
is followed. Define the function
$W\left(t\right)=S\left(t\right)+P\left(t\right)+T\left(t\right)+R\left(t\right)$ and
${}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\zeta \mathrm{},\wp}S\left(t\right)={}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\zeta \mathrm{},\mathrm{\wp}}S\left(t\right)+{}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\zeta \mathrm{},\mathrm{\wp}}P\left(t\right)+{}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\zeta \mathrm{},\mathrm{\wp}}T\left(t\right)+{}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\zeta \mathrm{},\mathrm{\wp}}R\left(t\right)$
$=B-\mathit{\alpha SP}+\mathit{\beta T}+\mathit{\delta R}-\mathit{\mu S}+\mathrm{}\mathit{\alpha SP}-\mathit{\gamma P}{-\alpha}_{1}P-\mathit{\mu P}$+${\alpha}_{1}P-\mathit{\beta T}-{\gamma}_{1}T-\mathit{\mu T}+\mathit{\gamma P}+{\gamma}_{1}T-\mathit{\delta R}-\mathit{\mu R}$
$=B-\left(S+P+T+R\right)\mu $
${}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\zeta \mathrm{},\mathrm{\wp}}W\left(t\right)+W\left(t\right)\mu \le B$
Applying the standard comparison theorem for the fractional order given in
[10] | Fazal Haq and Muhammad Shahzad, “Numerical Analysis of Fractional Order Epidemic Model of Childhood Disease”, Discrete Dynamics in Nature and Society, vol. 2017, Article Id 4057089, pp. 7, 2017. |
[10]
, we have,
$0\le W\left(t\right)\le W\left(0\right){E}_{\zeta}\left(-{\mu}^{\zeta}\right)=\frac{B}{\mu}{t}^{\zeta}{E}_{\zeta ,\mathrm{\wp}}\left(-{t}^{\zeta}\right),$
where ${E}_{\zeta}$ is the Mittag Leffler function. By taking $t\to \infty $, $0\le W\left(t\right)\le \frac{B}{\mu}$. Therefore, the solution is uniformly bounded in the region $R=\left\{\left(S,P,T,R\right)\in {R}_{+}^{4},W\left(t\right)\le \frac{B}{\mu}+\u03f5,\mathit{\u03f5}0\right\}.$
Now the non - non-negativity of the solution is studied for the fractional order system (
1). In (
1)
${}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\zeta \mathrm{},\mathrm{\wp}}S{\left(t\right)}_{S=0}=B,\mathrm{}{}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\zeta \mathrm{},\mathrm{\wp}}P{\left(t\right)}_{P=0}=0,\mathrm{}{}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\zeta \mathrm{},\mathrm{\wp}}T{\left(t\right)}_{T=0}=0,{}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\zeta \mathrm{},\mathrm{\wp}}R{\left(t\right)}_{R=0}=0$
Hence the solution of the system is non-negative.
4. Equilibrium and Stability Analysis
There are two equilibrium points exists for the model that are found by equating time derivatives in system (1) to zero as follows,
$B-\mathit{\alpha SP}+\mathit{\beta T}+\mathit{\delta R}-\mathit{\mu S}=0$
$\mathit{\alpha SP}-\mathit{\gamma P}{-\alpha}_{1}P-\mathit{\mu P}=0$
${\alpha}_{1}P-\mathit{\beta T}-{\gamma}_{1}T-\mathit{\mu T}=0$
$\mathit{\gamma P}+{\gamma}_{1}T-\mathit{\delta R}-\mathit{\mu R}=0$(5)
The two steady states of the model (
1) are given by,
Pollution free equilibrium point: ${P}_{0}\left(\frac{B}{\mu},0,0,0\right)$
Pollution extinct equilibrium point: ${P}_{*}({S}^{*},{P}^{*},{T}^{*},{R}^{*})$
where ${S}^{*}=\frac{\gamma +{\alpha}_{1}+\mu}{\alpha}$, ${T}^{*}=\frac{{\alpha}_{1}{P}^{*}}{\beta +{\gamma}_{1}+\mu}$, ${R}^{*}=\frac{1}{\delta +\mu}\left[\gamma +\frac{{\gamma}_{1}{\alpha}_{1}}{\beta +{\gamma}_{1}+\mu}\right]{P}^{*}$
${P}^{*}=\frac{\mu \left(\gamma +{\alpha}_{1}+\mu \right)}{\alpha}\left({R}_{0}-1\right){\left[\gamma +{\alpha}_{1}+\mu -\left(\frac{\beta {\alpha}_{1}}{\beta +{\gamma}_{1}+\mu}+\frac{1}{\delta}\left(\frac{1}{\delta +\mu}\left(\gamma +\frac{{\gamma}_{1}{\alpha}_{1}}{\beta +{\gamma}_{1}+\mu}\right)\right)\right)\right]}^{-1}$
${P}_{*}$ is positive, when (i) ${R}_{0}<1$ and $\gamma +{\alpha}_{1}+\mu <\left(\frac{\beta {\alpha}_{1}}{\beta +{\gamma}_{1}+\mu}+\frac{1}{\delta}\left(\frac{1}{\delta +\mu}\left(\gamma +\frac{{\gamma}_{1}{\alpha}_{1}}{\beta +{\gamma}_{1}+\mu}\right)\right)\right)$ (or)
(ii) ${R}_{0}>1$ and $\gamma +{\alpha}_{1}+\mu >\left(\frac{\beta {\alpha}_{1}}{\beta +{\gamma}_{1}+\mu}+\frac{1}{\delta}\left(\frac{1}{\delta +\mu}\left(\gamma +\frac{{\gamma}_{1}{\alpha}_{1}}{\beta +{\gamma}_{1}+\mu}\right)\right)\right)$
Basic Reproduction Number (BRN): [27] | Van den Driessche, Paul, and James Watmough. "Further notes on the basic reproduction number." Mathematical epidemiology (2008): 159-178. |
[27]
The next generation matrix (NGM) is used to find the BRN of the system (1)
Consider,${}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\mathit{\zeta},\wp}A\left(t\right)=F\left(Y\right)-V\left(Y\right)$, where$\mathit{F}\left(Y\right)=\left[\begin{array}{c}{f}_{1}\\ {f}_{2}\\ {f}_{3}\end{array}\right]=\left[\begin{array}{c}\mathit{\alpha SP}\\ 0\\ 0\end{array}\right]$
$V\left(Y\right)=\left[\begin{array}{c}{V}_{1}\\ {V}_{2}\\ {V}_{3}\end{array}\right]=\left[\begin{array}{c}-(\gamma +{\alpha}_{1}+\mu )\\ \mathit{\gamma P}+{\gamma}_{1}T-(\delta +\mu )R\\ {\alpha}_{1}P-(\beta +{\gamma}_{1}+\mu )T\end{array}\right]$
$F{\left(Y\right)}_{J}=\left[\begin{array}{c}\begin{array}{ccc}\frac{\partial {f}_{1}}{\partial P}& \frac{\partial {f}_{1}}{\partial T}& \frac{\partial {f}_{1}}{\partial R}\end{array}\\ \begin{array}{ccc}\frac{\partial {f}_{2}}{\partial P}& \frac{\partial {f}_{2}}{\partial T}& \frac{\partial {f}_{2}}{\partial R}\end{array}\\ \begin{array}{ccc}\frac{\partial {f}_{3}}{\partial P}& \frac{\partial {f}_{3}}{\partial T}& \frac{\partial {f}_{3}}{\partial R}\end{array}\end{array}\right]=\left[\begin{array}{c}\begin{array}{ccc}\mathit{\alpha S}& 0& 0\end{array}\\ \begin{array}{ccc}0& 0& 0\end{array}\\ \begin{array}{ccc}0& 0& 0\end{array}\end{array}\right]$
$V{\left(Y\right)}_{J}=\left[\begin{array}{c}\begin{array}{ccc}\frac{\partial {V}_{1}}{\partial P}& \frac{\partial {V}_{1}}{\partial T}& \frac{\partial {V}_{1}}{\partial R}\end{array}\\ \begin{array}{ccc}\frac{\partial {V}_{2}}{\partial P}& \frac{\partial {V}_{2}}{\partial T}& \frac{\partial {V}_{2}}{\partial R}\end{array}\\ \begin{array}{ccc}\frac{\partial {V}_{3}}{\partial P}& \frac{\partial {V}_{3}}{\partial T}& \frac{\partial {V}_{3}}{\partial R}\end{array}\end{array}\right]=\left[\begin{array}{ccc}-(\gamma +{\alpha}_{1}+\mu )& 0& 0\\ {\alpha}_{1}& -(\beta +{\gamma}_{1}+\mu )& 0\\ \gamma & {\gamma}_{1}& -(\delta +\mu )\end{array}\right]$
${V}^{-1}=\frac{1}{\mathit{ADC}}\left[\begin{array}{ccc}\mathit{DC}& 0& {\alpha}_{1}\gamma +\mathit{B\gamma}\\ {\alpha}_{1}C& \mathit{AC}& -\mathit{A\gamma}\\ {\alpha}_{1}{\gamma}_{1}+\mathit{D\gamma}& 0& \mathit{AD}\end{array}\right]$,
where $A=-(\gamma +{\alpha}_{1}+\mu ),\mathit{D}=-(\beta +{\gamma}_{1}+\mu ),\mathit{C}=-(\delta +\mu ).$It follows that spectral radius of the matrix is $\rho \left(F{V}^{-1}\right)=\mathrm{max}\left({\lambda}_{1,2,3}\right)$ at pollution free equilibrium point ${P}_{0}$. Therefore, the basic reproduction number is,
${\mathit{R}}_{0}=\frac{\mathit{\alpha B}}{\mu (\gamma +{\alpha}_{1}+\mu )}.$(6)
Theorem: 4.1 If all the eigen values have negative real part then the pollution free equilibrium point
${P}_{0}({S}_{0},0,0,0)$ of the system (
1) is locally asymptotically stable.
Proof: From the Matignon condition, we observe that the pollution free equilibrium point is locally asymptotically stable if and only if all the eigenvalues
${\lambda}_{i}$ of the Jacobian
$\left(J\right({P}_{0}\left)\right)$ satisfy
$|\mathrm{arg}({\lambda}_{i}\left)\right|>\frac{\mathit{a\pi}}{2}.$ The Jacobian matrix of the system (
1) is
$J\left[P\right]=\left[\begin{array}{cccc}\frac{\partial {F}_{1}}{\partial S}& \frac{\partial {F}_{1}}{\partial P}& \frac{\partial {F}_{1}}{\partial T}& \frac{\partial {F}_{1}}{\partial R}\\ \frac{\partial {F}_{2}}{\partial S}& \frac{\partial {F}_{2}}{\partial P}& \frac{\partial {F}_{2}}{\partial T}& \frac{\partial {F}_{2}}{\partial R}\\ \frac{\partial {F}_{3}}{\partial S}& \frac{\partial {F}_{3}}{\partial P}& \frac{\partial {F}_{3}}{\partial T}& \frac{\partial {F}_{3}}{\partial R}\\ \frac{\partial {F}_{4}}{\partial S}& \frac{\partial {F}_{4}}{\partial P}& \frac{\partial {F}_{4}}{\partial T}& \frac{\partial {F}_{4}}{\partial R}\end{array}\right]=\mathrm{}\left[\begin{array}{cccc}-\mathit{\alpha P}-\mu & \mathit{\alpha S}& \beta & \delta \\ \alpha & -(\gamma +{\alpha}_{1}+\mu )& 0& 0\\ 0& {\alpha}_{1}& -(\beta +{\gamma}_{1}+\mu )& 0\\ 0& \gamma & {\gamma}_{1}& -(\delta +\mu )\end{array}\right]$(7)
$J\left[{P}_{0}\right]=\left[\begin{array}{cccc}-\mu & \frac{\mathit{\alpha B}}{\mu}& \beta & \delta \\ 0& -(\gamma +{\alpha}_{1}+\mu )& 0& 0\\ 0& {\alpha}_{1}& -(\beta +\gamma +\mu )& 0\\ 0& \gamma & {\gamma}_{1}& -(\delta +\mu )\end{array}\right]$.
Here are the eigenvalues, ${\lambda}_{1}=-\mu ,\mathit{}{\lambda}_{2}=-\left(\delta +\mu \right),{\mathit{\lambda}}_{3,4}=\left[\begin{array}{cc}-\left(\beta +\gamma +\mu \right)& 0\\ {\gamma}_{1}& -(\delta +\mu )\end{array}\right]$. Therefore, the eigenvalues ${\lambda}_{3}=-\left(\beta +\gamma +\mu \right),{\mathit{\lambda}}_{4}=-(\delta +\mu )$. It can be observed that all the eigenvalues of the Jacobian matrix are negative. Hence by using the Matignon condition, we observe that the pollution-free equilibrium is locally asymptotically stable.
Theorem 4.2: The pollution extinction equilibrium point${\mathit{P}}_{*}\left({S}^{*},{P}^{*},{T}^{*},{R}^{*}\right)$ is globally asymptotically stable if and only if ${R}_{0}>1$ and $\frac{\mu}{S}+\frac{{\beta}^{2}T{T}^{*}}{4\alpha}+\frac{{\delta}^{2}R{R}^{*}}{4(\gamma +{\gamma}_{1})}\le \frac{{P}^{*}}{S}.$
Proof: Consider the following nonlinear Lyapunov function,
$W\left(t\right)=\left(S-{S}^{*}-{S}^{*}\mathrm{ln}\left(\frac{S}{{S}^{*}}\right)\right)+\left(P-{P}^{*}-{P}^{*}\mathrm{ln}\left(\frac{P}{{P}^{*}}\right)\right)+\left(T-{T}^{*}-{T}^{*}\mathrm{ln}\left(\frac{T}{{T}^{*}}\right)\right)+\left(R-{R}^{*}-{R}^{*}\mathrm{ln}\left(\frac{R}{{R}^{*}}\right)\right).$(8)
Take the Lyapunov fractional derivative on both sides,
$\frac{\mathit{dW}\left(t\right)}{\mathit{dt}}=\left(1-\frac{{S}^{*}}{S}\right)\stackrel{\u0307}{S}+\left(1-\frac{{P}^{*}}{P}\right)\stackrel{\u0307}{P}+\left(1-\frac{{T}^{*}}{T}\right)\stackrel{\u0307}{T}+\left(1-\frac{{R}^{*}}{R}\right)\stackrel{\u0307}{R}\mathrm{}$(9)
$=\left(1-\frac{{S}^{*}}{S}\right)\left[B-\mathit{\alpha SP}+\mathit{\beta T}+\mathit{\delta R}-\mathit{\mu S}\right]+\left(1-\frac{{P}^{*}}{P}\right)\left[\mathit{\alpha SP}-\left(\gamma +{\alpha}_{1}+\mu \right)P\right]$
$+\left(1-\frac{{T}^{*}}{T}\right)\left[{\alpha}_{1}P-\left(\beta +{\gamma}_{1}+\mu \right)T\right]+\left(1-\frac{{R}^{*}}{R}\right)[\mathit{\gamma P}+{\gamma}_{1}T-(\delta +\mu \left)R\right]$
$=\left(1-\frac{{S}^{*}}{S}\right)\left[\alpha {S}^{*}{P}^{*}+\beta {T}^{*}+\delta {R}^{*}-\mu {S}^{*}-\mathit{\alpha SP}+\mathit{\beta T}+\mathit{\delta R}-\mathit{\mu S}\right]+\left(1-\frac{{P}^{*}}{P}\right)\left[\mathit{\alpha SP}-\alpha {S}^{*}P\right]$
$+\left(1-\frac{{T}^{*}}{T}\right)\left[{\alpha}_{1}P-\frac{{\alpha}_{1}{P}^{*}}{{T}^{*}}T\right]+\left(1-\frac{{R}^{*}}{R}\right)\left[\mathit{\gamma P}+{\gamma}_{1}T-\left(\frac{\gamma {P}^{*}}{{R}^{*}}+\frac{{\gamma}_{1}{T}^{*}}{{R}^{*}}\mathrm{}\right)R\right]$
After some simplification we obtain,
$=-{\left[\sqrt{\alpha}\left(\frac{{(T}^{*}-T)}{\sqrt{T{T}^{*}}}-\frac{\beta \sqrt{T{T}^{*}}}{2\sqrt{\alpha}}\left(S-{S}^{*}\right)\right)\right]}^{2}-{\left[\frac{\sqrt{\gamma +{\gamma}_{1}}}{R{R}^{*}}\left({R}^{*}-R\right)-\frac{\delta \sqrt{R{R}^{*}}}{2\sqrt{\gamma +{\gamma}_{1}}}\left({S}^{*}-S\right)\right]}^{2}$
$+\left[-\frac{{P}^{*}}{S}+\mathrm{}\frac{\mu}{S}+\frac{{\beta}^{2}T{T}^{*}}{4\alpha}+\frac{{\delta}^{2}R{R}^{*}}{4\left(\gamma +{\gamma}_{1}\right)}\right]{\left(S-{S}^{*}\right)}^{2}$
Finally, we conclude that, ${}_{t}{}^{\mathit{CFF}}{\mathbb{D}}^{\mathit{\zeta},\wp}W\left(t\right)\le 0$, if and only if $\frac{\mu}{S}+\frac{{\beta}^{2}T{T}^{*}}{4\alpha}+\frac{{\delta}^{2}R{R}^{*}}{4(\gamma +{\gamma}_{1})}\le \frac{{P}^{*}}{S}.$
The first derivative of a Lyapunov function is instrumental in evaluating the global stability of endemic equilibrium locations. It provides essential insights that can be further enriched by second derivative analysis. While the second derivative informs us about curvature through its sign, the first derivative offers valuable information on disease propagation. Examining the first derivative gives us insights into how the system evolves and spreads over time. Complemented by the second derivative examination, this analysis contributes to a comprehensive understanding of the dynamics and stability of endemic equilibria in disease models.
Taking the second derivative of (
9),
$\frac{d\stackrel{\u0307}{W}\left(t\right)}{\mathit{dt}}=\frac{d}{\mathit{dt}}\left(\left(1-\frac{{S}^{*}}{S}\right)\stackrel{\u0307}{S}\right)+\frac{d}{\mathit{dt}}\left(\left(1-\frac{{P}^{*}}{P}\right)\stackrel{\u0307}{P}\right)+\frac{d}{\mathit{dt}}\left(\left(1-\frac{{T}^{*}}{T}\right)\stackrel{\u0307}{T}\right)+\frac{d}{\mathit{dt}}\left(\left(1-\frac{{R}^{*}}{R}\right)\stackrel{\u0307}{R}\right)$
$={\left(\frac{\stackrel{\u0307}{S}}{S}\right)}^{2}{S}^{*}+{\left(\frac{\stackrel{\u0307}{P}}{P}\right)}^{2}{P}^{*}+{\left(\frac{\stackrel{\u0307}{T}}{T}\right)}^{2}{T}^{*}+{\left(\frac{\stackrel{\u0307}{R}}{R}\right)}^{2}{R}^{*}+\left(1-\frac{{S}^{*}}{S}\right)\stackrel{\u0308}{S}+\left(1-\frac{{P}^{*}}{P}\right)\stackrel{\u0308}{P}+\left(1-\frac{{T}^{*}}{T}\right)\stackrel{\u0308}{T}+\left(1-\frac{{R}^{*}}{R}\right)\stackrel{\u0308}{R}$
where,
$\stackrel{\u0308}{S}=-\alpha \stackrel{\u0307}{S}P-\alpha \stackrel{\u0307}{P}S+\beta \stackrel{\u0307}{T}+\delta \stackrel{\u0307}{R}-\mu \stackrel{\u0307}{S}$
$\stackrel{\u0308}{P}=\alpha \stackrel{\u0307}{S}P+\alpha \stackrel{\u0307}{P}S-\gamma \stackrel{\u0307}{P}-{\alpha}_{1}\stackrel{\u0307}{P}-\mu -\gamma \stackrel{\u0307}{P}$
$\stackrel{\u0308}{T}={\alpha}_{1}\stackrel{\u0307}{P}-\beta \stackrel{\u0307}{T}-{\gamma}_{1}\stackrel{\u0307}{T}-\mu \stackrel{\u0307}{T}$
$\stackrel{\u0308}{R}=\gamma \stackrel{\u0307}{P}+{\gamma}_{1}\stackrel{\u0307}{T}-\mathit{\delta R}-\mathit{\mu R}$
$\frac{d\stackrel{\u0307}{W}\left(t\right)}{\mathit{dt}}={\left(\frac{\stackrel{\u0307}{S}}{S}\right)}^{2}{S}^{*}+{\left(\frac{\stackrel{\u0307}{P}}{P}\right)}^{2}{P}^{*}+{\left(\frac{\stackrel{\u0307}{T}}{T}\right)}^{2}{T}^{*}+{\left(\frac{\stackrel{\u0307}{R}}{R}\right)}^{2}{R}^{*}+\left(1-\frac{{S}^{*}}{S}\right)\left[-\alpha \stackrel{\u0307}{S}P-\alpha \stackrel{\u0307}{P}S+\beta \stackrel{\u0307}{T}+\delta \stackrel{\u0307}{R}-\mu \stackrel{\u0307}{S}\right]$
$+\left(1-\frac{{P}^{*}}{P}\right)\left[\alpha \stackrel{\u0307}{S}P+\alpha \stackrel{\u0307}{P}S-\gamma \stackrel{\u0307}{P}-{\alpha}_{1}\stackrel{\u0307}{P}-\mu -\gamma \stackrel{\u0307}{P}\right]+\left(1-\frac{{T}^{*}}{T}\right)\left[{\alpha}_{1}\stackrel{\u0307}{P}-\beta \stackrel{\u0307}{T}-{\gamma}_{1}\stackrel{\u0307}{T}-\mu \stackrel{\u0307}{T}\right]$
$+\mathrm{}\left(1-\frac{{R}^{*}}{R}\right)\left[\gamma \stackrel{\u0307}{P}+{\gamma}_{1}\stackrel{\u0307}{T}-\mathit{\delta R}-\mathit{\mu R}\right]$
$={\left(\frac{\stackrel{\u0307}{S}}{S}\right)}^{2}{S}^{*}+{\left(\frac{\stackrel{\u0307}{P}}{P}\right)}^{2}{P}^{*}+{\left(\frac{\stackrel{\u0307}{T}}{T}\right)}^{2}{T}^{*}+{\left(\frac{\stackrel{\u0307}{R}}{R}\right)}^{2}{R}^{*}-\mu \left[\stackrel{\u0307}{S}+\stackrel{\u0307}{P}+\stackrel{\u0307}{T}+\stackrel{\u0307}{R}\right]-\left(\frac{{S}^{*}}{S}\right)\left[-\alpha \stackrel{\u0307}{S}P-\alpha \stackrel{\u0307}{P}S+\beta \stackrel{\u0307}{T}+\delta \stackrel{\u0307}{R}-\mu \stackrel{\u0307}{S}\right]$
$-\left(\frac{{P}^{*}}{P}\right)\left[\alpha \stackrel{\u0307}{S}P+\alpha \stackrel{\u0307}{P}S-\gamma \stackrel{\u0307}{P}-{\alpha}_{1}\stackrel{\u0307}{P}-\mu -\gamma \stackrel{\u0307}{P}\right]-\left(\frac{{T}^{*}}{T}\right)\left[{\alpha}_{1}\stackrel{\u0307}{P}-\beta \stackrel{\u0307}{T}-{\gamma}_{1}\stackrel{\u0307}{T}-\mu \stackrel{\u0307}{T}\right]-\left(\frac{{R}^{*}}{R}\right)\left[\gamma \stackrel{\u0307}{P}+{\gamma}_{1}\stackrel{\u0307}{T}-\mathit{\delta R}-\mathit{\mu R}\right]$
Now, substituting the values of the derivative $\stackrel{\u0307}{S},\stackrel{\u0307}{P},\stackrel{\u0307}{T},\stackrel{\u0307}{R}$ and rearranging the terms of positive and negative groups, so that we have,
$\frac{{d}^{2}W\left(t\right)}{\mathit{dt}}={\mathcal{H}}_{1}-{\mathcal{H}}_{2}$
Therefore, we see that,
i. If ${\mathcal{H}}_{1}{\mathcal{H}}_{2}$ then $\frac{{d}^{2}W\left(t\right)}{\mathit{dt}}>0,$
ii. If ${\mathcal{H}}_{1}<{\mathcal{H}}_{2}$ then $\frac{{d}^{2}W\left(t\right)}{\mathit{dt}}<0,$
iii. If ${\mathcal{H}}_{1}={\mathcal{H}}_{2}$ then $\frac{{d}^{2}W\left(t\right)}{\mathit{dt}}=0,$
This interpretation explains the model's global behavior with respect to Lyapunov second order derivative.
Lemma 4.1: The given model (1) has no periodic orbits.
Proof: The proof has carried with a Dulac plus Poincare Bendixson theorem as follows:
$h\left(S,P\right)=\frac{1}{\mathit{SP}}$ , where all the parameters are positive.
$\nabla .\mathit{hf}=\frac{\partial}{\partial S}\left(h.{f}_{1}\right)+\frac{\partial}{\partial P}\left(h.{f}_{2}\right)$
$=\frac{\partial}{\partial S}(h.\left(B-\mathit{\alpha SP}+\mathit{\beta T}+\mathit{\delta R}-\mathit{\mu S}\right))+\frac{\partial}{\partial P}\left(h.(\mathit{\alpha SP}-\mathit{\gamma P}{-\alpha}_{1}P-\mathit{\mu P}\right))$
$=-\left[\frac{B}{{S}^{2}P}+\frac{\mathit{\beta T}}{{S}^{2}P}+\frac{\mathit{\delta R}}{{S}^{2}P}\right]<0$(10)
Hence, from the Dulac principle, we conclude that there is no periodic solution for the given region R.
5. Numerical Approximation
The numerical approximate solution of a system (
1) is obtained using Laplace adomain decomposition method,
Applying Laplace transform on both sides of equation (
1),
Substitute the initial condition
(2) in the above equations,
${s}^{\alpha}L\left[S\left(t\right)\right]-{s}^{\alpha -1}S\left(0\right)=L\left(B-\mathit{\alpha SP}+\mathit{\beta T}+\mathit{\delta R}-\mathit{\mu S}\right)$
${s}^{\alpha}L\left[P\left(t\right)\right]-{s}^{\alpha -1\mathrm{}}P\left(0\right)=L(\mathit{\alpha SP}-(\gamma +{\alpha}_{1}+\mu \left)P\right)$
${s}^{\alpha}L\left(T\left(t\right)\right)-{s}^{\alpha -1}T\left(0\right)=L\left({\alpha}_{1}P-\left(\beta +{\gamma}_{1}+\mu \right)T\right)$
${s}^{\alpha}L\left(R\left(t\right)\right)-{s}^{\alpha -1}R\left(0\right)=L(\mathit{\gamma P}+{\gamma}_{1}T-(\delta +\mu \left)R\right)$
Now,
$L\left[S\left(t\right)\right]=\frac{S\left(0\right)}{s}+\frac{1}{{s}^{\alpha}}L\left(B-\mathit{\alpha SP}+\mathit{\beta T}+\mathit{\delta R}-\mathit{\mu S}\right)$
$L\left[P\left(t\right)\right]=\frac{P\left(0\right)}{s}+\frac{1}{{s}^{\alpha}}L(\mathit{\alpha SP}-(\gamma +{\alpha}_{1}+\mu \left)P\right)$
$L\left(T\left(t\right)\right)=\frac{T\left(0\right)}{s}+\frac{1}{{s}^{\alpha}}L\left({\alpha}_{1}P-\left(\beta +{\gamma}_{1}+\mu \right)T\right)$
$L\left(R\left(t\right)\right)=\frac{R\left(0\right)}{s}+\frac{1}{{s}^{\alpha}}L(\mathit{\gamma P}+{\gamma}_{1}T-(\delta +\mu \left)R\right)$(11)
Consider the infinite series,
$S\left(t\right)=\sum _{k=0}^{\mathrm{\infty}}{S}_{k}\left(t\right),P\left(t\right)=\sum _{k=0}^{\mathrm{\infty}}{P}_{k}\left(t\right),\mathrm{}T\left(t\right)=\sum _{k=0}^{\mathrm{\infty}}{T}_{k}\left(t\right),R\left(t\right)=\sum _{k=0}^{\mathrm{\infty}}{R}_{k}\left(t\right)$(12)
Also consider the decomposition of nonlinear terms as follows,
$S\left(t\right)P\left(t\right)=\sum _{k=0}^{\infty}{A}_{k}\left(t\right)$(13)
where the decomposition polynomial ${A}_{k}$ is defined as follows,
${A}_{k}=\frac{1}{\mathrm{\u1d26}(k+1)}\frac{{d}^{k}}{d{\lambda}^{k}}{\left[\sum _{i=0}^{k}{\lambda}^{i}{S}_{i}\left(t\right)\sum _{i=0}^{k}{\lambda}^{i}{P}_{i}\left(t\right)\right]}_{\lambda =0\mathrm{}}$
The first three polynomials are,
${A}_{0}={S}_{0}\left(t\right){P}_{0}\left(t\right),$
${A}_{1}={S}_{0}\left(t\right){P}_{1}\left(t\right)+{S}_{1}\left(t\right){P}_{0}\left(t\right)$
${A}_{2}=2{S}_{0}\left(t\right){P}_{2}\left(t\right)+2{S}_{1}\left(t\right){P}_{1}\left(t\right)+2{S}_{2}\left(t\right){P}_{0}\left(t\right)$
Substituting above equations in equation (
11),
$L\left[\sum _{k=0}^{\mathrm{\infty}}{S}_{k}\left(t\right)\right]=\frac{S\left(0\right)}{s}+\frac{1}{{s}^{\alpha}}L\left(B-\alpha \sum _{k=0}^{\mathrm{\infty}}{A}_{k}\left(t\right)+\beta \sum _{k=0}^{\mathrm{\infty}}{T}_{k}\left(t\right)+\delta \sum _{k=0}^{\mathrm{\infty}}{R}_{k}\left(t\right)-\mu \sum _{k=0}^{\mathrm{\infty}}{S}_{k}\left(t\right)\right)$
$L\left[\sum _{k=0}^{\mathrm{\infty}}{P}_{k}\left(t\right)\right]=\frac{P\left(0\right)}{s}+\frac{1}{{s}^{\alpha}}L(\alpha \sum _{k=0}^{\mathrm{\infty}}{A}_{k}\left(t\right)-(\gamma +{\alpha}_{1}+\mu \left)\sum _{k=0}^{\mathrm{\infty}}{P}_{k}\left(t\right)\right)$
$L\left(\sum _{k=0}^{\mathrm{\infty}}{T}_{k}\left(t\right)\right)=\frac{T\left(0\right)}{s}+\frac{1}{{s}^{\alpha}}L\left({\alpha}_{1}\sum _{k=0}^{\mathrm{\infty}}{P}_{k}\left(t\right)-\left(\beta +{\gamma}_{1}+\mu \right)\sum _{k=0}^{\mathrm{\infty}}{T}_{k}\left(t\right)\right)$
$L\left(\sum _{k=0}^{\mathrm{\infty}}{R}_{k}\left(t\right)\right)=\frac{R\left(0\right)}{s}+\frac{1}{{s}^{\alpha}}L(\gamma \sum _{k=0}^{\mathrm{\infty}}{P}_{k}\left(t\right)+{\gamma}_{1}\sum _{k=0}^{\mathrm{\infty}}{T}_{k}\left(t\right)-(\delta +\mu \left)\sum _{k=0}^{\mathrm{\infty}}{R}_{k}\left(t\right)\right)$(14)
Comparing both sides of equation of (
14),
$L\left({S}_{0}\right)=\frac{{N}_{1}}{s},$
$L\left({S}_{1}\right)=\frac{1}{{s}^{\alpha}}B-\frac{\alpha}{{s}^{\alpha}}L\left({A}_{0}\right)+\frac{\beta}{{s}^{\alpha}}L\left({T}_{0}\right)+\frac{\delta}{{s}^{\alpha}}L\left({R}_{0}\right)-\frac{\mu}{{s}^{\alpha}}L\left({S}_{0}\right)$
$L\left({S}_{2}\right)=\frac{1}{{s}^{\alpha}}B-\frac{\alpha}{{s}^{\alpha}}L\left({A}_{1}\right)+\frac{\beta}{{s}^{\alpha}}L\left({T}_{1}\right)+\frac{\delta}{{s}^{\alpha}}L\left({R}_{1}\right)-\frac{\mu}{{s}^{\alpha}}L\left({S}_{1}\right)$
In general,$L\left({S}_{k+1}\right)=\frac{1}{{s}^{\alpha}}B-\frac{\alpha}{{s}^{\alpha}}L\left({A}_{k}\right)+\frac{\beta}{{s}^{\alpha}}L\left({T}_{k}\right)+\frac{\delta}{{s}^{\alpha}}L\left({R}_{k}\right)-\frac{\mu}{{s}^{\alpha}}L\left({S}_{k}\right)$(15)
$L\left({P}_{0}\right)=\frac{{N}_{2}}{s},$
$L\left({P}_{1}\right)=\frac{\alpha}{{s}^{\alpha}}L\left({A}_{0}\right)-\frac{\gamma +{\alpha}_{1}+\mu}{{s}^{\alpha}}L\left({P}_{0}\right))$
$L\left({P}_{2}\right)=\frac{\alpha}{{s}^{\alpha}}L\left({A}_{1}\right)-\frac{\gamma +{\alpha}_{1}+\mu}{{s}^{\alpha}}L\left({P}_{1}\right))$
In general,$L\left({P}_{k+1}\right)=\frac{\alpha}{{s}^{\alpha}}L\left({A}_{k}\right)-\frac{\gamma +{\alpha}_{1}+\mu}{{s}^{\alpha}}L\left({P}_{k}\right))$(16)
$L\left({T}_{0}\right)=\frac{{N}_{3}}{s},$
$L\left({T}_{1}\right)=\frac{1}{{s}^{\alpha}}{\alpha}_{1}L\left({P}_{0}\right)-\frac{\beta +{\gamma}_{1}+\mu}{{s}^{\alpha}}L\left({T}_{0}\right)$
$L\left({T}_{2}\right)=\frac{1}{{s}^{\alpha}}{\alpha}_{1}L\left({P}_{1}\right)-\frac{\beta +{\gamma}_{1}+\mu}{{s}^{\alpha}}L\left({T}_{1}\right)$
In general,$L\left({T}_{k+1}\right)=\frac{1}{{s}^{\alpha}}{\alpha}_{1}L\left({P}_{k}\right)-\frac{\beta +{\gamma}_{1}+\mu}{{s}^{\alpha}}L\left({T}_{k}\right))$(17)
$L\left({R}_{0}\right)=\frac{{N}_{4}}{s},$
$L\left({R}_{1}\right)=\frac{1}{{s}^{\alpha}}\gamma \mathrm{}L\left({P}_{0}\right)+\frac{1}{{s}^{\alpha}}{\gamma}_{1}\mathrm{}L\left({T}_{0}\right)-\frac{1}{{s}^{\alpha}}(\delta +\mu )L\left({R}_{0}\right)$
$L\left({R}_{2}\right)=\frac{1}{{s}^{\alpha}}\gamma \mathrm{}L\left({P}_{1}\right)+\frac{1}{{s}^{\alpha}}{\gamma}_{1}\mathrm{}L\left({T}_{1}\right)-\frac{1}{{s}^{\alpha}}(\delta +\mu )L\left({R}_{1}\right)$
In general,$\mathrm{}L\left({R}_{k+1}\right)=\frac{1}{{s}^{\alpha}}\gamma \mathrm{}L\left({P}_{k}\right)+\frac{1}{{s}^{\alpha}}{\gamma}_{1}\mathrm{}L\left({T}_{k}\right)-\frac{1}{{s}^{\alpha}}(\delta +\mu )L\left({R}_{k}\right)$(18)
Taking Laplace inverse transform (LIT) on both sides of the above equations,
${S}_{1}=\frac{B{t}^{\alpha}}{\Gamma (\alpha +1)}-\frac{{t}^{\alpha}}{\Gamma \left(\alpha +1\right)}[\alpha {N}_{1}{N}_{2}+\beta {N}_{3}+\delta {N}_{4}-\mu {N}_{1}]$
${S}_{2}=\frac{B{t}^{\alpha}}{\Gamma (\alpha +1)}-\frac{{t}^{\alpha}}{\Gamma \left(\alpha +1\right)}[\alpha {A}_{1}+\beta {T}_{1}+\delta {R}_{1}-\mu {S}_{1}]$
${P}_{1}=\frac{{t}^{\alpha}}{\Gamma \left(\alpha +1\right)}\left[\alpha {N}_{1}{N}_{2}-\left(\gamma +{\alpha}_{1}+\mu \right){N}_{2}\right]$
${P}_{2}=\frac{{t}^{\alpha}}{\Gamma \left(\alpha +1\right)}\left[\alpha {A}_{1}-\left(\gamma +{\alpha}_{1}+\mu \right){N}_{1}\right]$
${T}_{1}=\frac{{t}^{\alpha}}{\Gamma \left(\alpha +1\right)}\left[{\alpha}_{1}{N}_{2}-\left(\beta +{\gamma}_{1}+\mu \right){N}_{3}\right]$
${T}_{2}=\frac{{t}^{\alpha}}{\Gamma \left(\alpha +1\right)}\left[{\alpha}_{1}{P}_{1}-\left(\beta +{\gamma}_{1}+\mu \right){T}_{1}\right]$
${R}_{1}=\frac{{t}^{\alpha}}{\Gamma \left(\alpha +1\right)}\left[\gamma {N}_{2}+{\gamma}_{1}{N}_{3}-\left(\delta +\mu \right){N}_{4}\right]$
${R}_{2}=\frac{{t}^{\alpha}}{\Gamma \left(\alpha +1\right)}\left[\gamma {P}_{1}+{\gamma}_{1}{T}_{1}-\left(\delta +\mu \right){R}_{1}\right].$
Then we have,$S\left(t\right)={N}_{1}+{S}_{1}+{S}_{2}+{S}_{3}+\cdots $
$P\left(t\right)={N}_{2}+{P}_{1}+{P}_{2}+{P}_{3}+\cdots $
$T\left(t\right)={N}_{3}+{T}_{1}+{T}_{2}+{T}_{3}+\cdots $
$R\left(t\right)={N}_{4}+{R}_{1}+{R}_{2}+{R}_{3}+\cdots $
Substituting the initial and parameter values, we get the required numerical approximate solution for the equation (
1).
6. Sensitivity Analysis
6.1. Local Sensitivity Analysis
The sensitivity analysis of each parameter of the basic reproduction number with respect to the model statements is performed in order to control the pollution spread through the soil. Sensitivity indices of each parameter in the BRN are calculated to implement the sensitivity analysis
[27] | Van den Driessche, Paul, and James Watmough. "Further notes on the basic reproduction number." Mathematical epidemiology (2008): 159-178. |
[27]
. The normalized forward sensitivity indices are obtained for different behavior of each parameter in prevalence of pollution transmission. Let M be a normalized forward index variables with respect to the dependent variable
$\rho $ and is defined as follows,
${\mathrm{\aleph}}_{\rho}^{M}=\frac{\partial M}{\partial \rho}\times \frac{\rho}{M}.$
By using explicit form of basic reproduction number, we can find an analytical expression for various parameters of the normalized forward sensitivity indices. The normalized sensitivity analysis at the baseline for each parameter are predicted as follows,
${\mathrm{\aleph}}_{B}^{{R}_{0}}=\frac{\partial {R}_{0}}{\partial B}\times \frac{B}{{R}_{0}}=1$
${\mathrm{\aleph}}_{\alpha}^{{R}_{0}}=\frac{\partial {R}_{0}}{\partial \alpha}\times \frac{\alpha}{{R}_{0}}=1$
${\mathrm{\aleph}}_{\mu}^{{R}_{0}}=\frac{\partial {R}_{0}}{\partial \mu}\times \frac{\mu}{{R}_{0}}=-\frac{(1+\gamma +{\alpha}_{1}+\mu )}{\gamma +{\alpha}_{1}+\mu}$
${\mathrm{\aleph}}_{{\alpha}_{1}}^{{R}_{0}}=\frac{\partial {R}_{0}}{\partial {\alpha}_{1}}\times \frac{{\alpha}_{1}}{{R}_{0}}=-\frac{{\alpha}_{1}}{(\gamma +{\alpha}_{1}+\mu )}$
${\mathrm{\aleph}}_{\gamma}^{{R}_{0}}=\frac{\partial {R}_{0}}{\partial \gamma}\times \frac{\gamma}{{R}_{0}}=-\frac{\gamma}{(\gamma +{\alpha}_{1}+\mu )}$
Sensitivity analysis is used to identify the important aspects of different parameters contribution in the results of the basic reproduction number which are based on their uncertainty estimation. Partial rank correlation is the acceleration method to find the statistical influence of the monotonicity relationship between any parameters and the basic reproduction number. The performance sensitivity index is presented in the above calculations which exhibit that positive signs are presented to be very proportionally and highly sensitive to the parameter values of
${R}_{0}$ at the same time those with the negative sign predict less sensitivity to decrease
${R}_{0}$ and other parameters are neutrally sensitive. The positive estimated parameters B and
$\alpha $ have been increased (decrease) by 10 % and with the same percentage
${R}_{0}$ is also increased (decrease). The remaining negative parameters decrease (increase) by 10% and the basic reproduction number
${R}_{0}$ by increase (decrease) are estimated in
Table 1.
The irrespective effect of the control parameter in basic reproduction number
${R}_{0}$ for the most sensitive parametric representation is illustrated through different contour plots which are shown in
Figures 8 and 9. In addition, the
Figure 7 contour plot shows that increasing the necessary expected remediation rate for the pollution in the soil substance will ultimately reduce the invested amount of secondary pollution number. At the same time, the increased amount of pollution transmission has increased the secondary pollution number. Also, the reduction of pollution transferred into susceptible will decline the basic reproduction number
${R}_{0}.$ Figure 8 shows that when the amount gradual increase of recovery rate along with pollution transmission rate acknowledges the reduction of
${R}_{0}.$ In addition, the decline rate of transmission in the soil saturation has decreased the BRN.
Table 1. The evaluated parameters of the model (1) at baseline values for the normalized forward sensitivity indices of ${R}_{0}$.
Parameter | Sensitivity indices | Critical values |
B | ${\aleph}_{B}^{{R}_{0}}$ | 1 |
$\alpha $ | ${\aleph}_{\alpha}^{{R}_{0}}$ | 1 |
${\alpha}_{1}$ | ${\aleph}_{{\alpha}_{1}}^{{R}_{0}}$ | -0.7741025239 |
$\gamma $ | ${\aleph}_{\gamma}^{{R}_{0}}$ | -0.1557454146 |
$\mu $ | ${\mathrm{\aleph}}_{\mu}^{{R}_{0}}$ | -1.7824 |
6.2. Global Sensitivity Analysis
The selection of three sampling sites within Coimbatore city, particularly around the 1194 electroplating facilities primarily situated on black soils, underscored the interconnectedness of industrial and agricultural landscapes in the region. The strategic selection of sampling sites near industries and within 1km distance from primary waste and effluent disposal points facilitated a comprehensive assessment of potential soil contamination, highlighting the urgent need for remedial actions to mitigate the adverse effects of heavy metal pollution on both environmental and public health fronts.
[22] | Sherene T., “Heavy Metal status of soils in Industrial Belts of Coimbatore District, Tamil Nadu”, Nature environment and pollution technology, Vol. 8, No. 3, pp. 613-618, 2009. |
[22]
Table 2. Analyzing Sensitivity indices of Cr using the Morris method and Normalized Forward sensitivity indices (${R}_{0}=0.12857<1$).
Parameter | Values 22] | Morris Method | NFSI | Ranges |
B | 0.0234 | 0 | 0.010523 | 0.1-0.3 |
$\alpha $ | 0.2 | 0.0026402 | 0.086314 | 0.1-0.5 |
${\alpha}_{1}$ | 0.9 | 0.35434 | -0.16239 | 0.05-0.15 |
$\beta $ | 0.0233 | 0.25161 | 0.007436 | 0.01-0.03 |
$\gamma $ | 0.9 | 0.0076827 | 0.50965 | 0.05-0.15 |
${\gamma}_{1}$ | 0.35 | 0.23044 | 0.52066 | 0.05-0.15 |
$\delta $ | 0.05 | 0.14256 | -0.49234 | 0.1-0.3 |
$\mu $ | 0.02 | 0.010727 | -1.0433 | 0.01-0.03 |
Table 3. Analyzing Sensitivity indices of Ni using the Morris method and Normalized Forward sensitivity indices (${R}_{0}=0.45<1$).
Parameter | Values 22] | Morris Method | NFSI | Ranges |
B | 0.0234 | 0 | 0.01298 | 0.1-0.3 |
$\alpha $ | 0.7 | 0.0026259 | 0.033217 | 0.2-0.7 |
${\alpha}_{1}$ | 0.9 | 0.35209 | -0.13826 | 0.05-0.15 |
$\beta $ | 0.0233 | 0.25952 | 0.04273 | 0.01-0.03 |
$\gamma $ | 0.9 | 0.0072494 | 0.52432 | 0.05-0.15 |
${\gamma}_{1}$ | 0.35 | 0.22444 | 0.52874 | 0.05-0.15 |
$\delta $ | 0.05 | 0.14363 | -0.48483 | 0.1-0.3 |
$\mu $ | 0.02 | 0.010451 | -1.0026 | 0.01-0.03 |
Figure 2. Parameter behaviors in terms Cd transmission using the sensitivity assessment.
Table 4. Analyzing Sensitivity indices of Cd using Morris method and Normalized Forward sensitivity indices (${R}_{0}=0.57857<1$).
Parameter | Values 22] | Morris Method | NFSI | Ranges |
B | 0.0234 | 0 | 0.013388 | 0.1-0.3 |
$\alpha $ | 0.9 | 0.002596 | 0.018777 | 0.5-0.9 |
${\alpha}_{1}$ | 0.9 | 0.34795 | -0.13806 | 0.05-0.15 |
$\beta $ | 0.0233 | 0.26974 | 0.053694 | 0.01-0.03 |
$\gamma $ | 0.9 | 0.0070178 | 0.52282 | 0.05-0.15 |
${\gamma}_{1}$ | 0.35 | 0.22026 | 0.53043 | 0.05-0.15 |
$\delta $ | 0.05 | 0.14165 | -0.48194 | 0.1-0.3 |
$\mu $ | 0.02 | 0.010788 | -1.0433 | 0.01-0.03 |
Figure 3. Parameter behaviors in terms Ni transmission using the sensitivity assessment.
Figure 4. Parameter behaviors in terms Cr transmission using the sensitivity assessment.