Next Article in Journal
A New Look on the Profitability of Fixed and Indexed Mortgage Products
Next Article in Special Issue
An Optimal Control Problem Related to the RSS Model
Previous Article in Journal
Security Quantification for Discrete Event Systems Based on the Worth of States
Previous Article in Special Issue
Optimal Melanoma Treatment Protocols for a Bilinear Control Model
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Application of an Optimal Control Therapeutic Approach for the Memory-Regulated Infection Mechanism of Leprosy through Caputo–Fabrizio Fractional Derivative

1
School of Mathematics and Statistics, Beijing Technology and Business University, Beijing 100048, China
2
Centre for Mathematical Biology and Ecology, Department of Mathematics, Jadavpur University, Kolkata 700032, West Bengal, India
3
Department of Statistics, Visva-Bharati University, Santiniketan 731235, West Bengal, India
*
Author to whom correspondence should be addressed.
Mathematics 2023, 11(17), 3630; https://doi.org/10.3390/math11173630
Submission received: 27 July 2023 / Revised: 12 August 2023 / Accepted: 15 August 2023 / Published: 22 August 2023

Abstract

:
Leprosy (Hansen’s disease) is an infectious, neglected tropical skin disease caused by the bacterium Mycobacterium leprae (M. leprae). It is crucial to note that the dynamic behavior of any living microorganism such as M. leprae not only depends on the conditions of its current state (e.g., substrate concentration, medium condition, etc.) but also on those of its previous states. In this article, we have developed a three-dimensional mathematical model involving concentrations of healthy Schwann cells, infected Schwann cells, and M. leprae bacteria in order to predict the dynamic changes in the cells during the disease dissemination process; additionally, we investigated the effect of memory on system cell populations, especially on the M. leprae bacterial population, by analyzing the Caputo–Fabrizio fractionalized version of the model. Most importantly, we developed and investigated a fractionalized optimal-control-induced system comprising the combined drug dose therapy of Ofloxacin and Dapsone intended to achieve a more realistic treatment regime for leprosy. The main goal of our research article is to compare this fractional-order system with the corresponding integer-order model and also to distinguish the rich dynamics exhibited by the optimal-control-induced system based on different values of the fractional order ζ ( 0 , 1 ) . All of the analytical results are validated through proper numerical simulations and are compared with some real clinical data.

1. Introduction

Leprosy (also known as Hensen’s disease) is an age-old disease and is described in the literature of ancient civilizations. It is a chronic infectious disease caused by a type of bacteria called Mycobacterium leprae (M. leprae), which primarily affects the peripheral nerves, skin, and eyes, and finally leads to loss of ability to sense and touch accompanied by dry, flaky, reddish skin due to inflammation. In very advanced cases, apparent loss of various human organs, such as toes and fingers, and permanent blindness are observed [1]. The treatment of leprosy has certainly been advanced by the introduction of multidrug therapy (MDT) in 1981, following the recommendation of the World Health Organization (WHO) [2,3]. The number of new leprosy cases detected annually has remained quite stable over the last 15 years [4]. According to the WHO report of 2022, about 1500 people in the United States and 250,000 people around the world, especially India, Indonesia, and Brazil, are affected by this disease every year. This situation clearly indicates that, far from being eliminated as a public health problem, leprosy still causes considerable long-term morbidity in both the developing and developed worlds. After reviewing the various issues related to classification, treatment regimes, drug resistance, and existing elimination strategies, it has been observed that new anti-leprosy agents and treatment designs are required in the near future [5]. In order to investigate the pathogenesis of leprosy from a completely novel analytical and numerical point of view, considering only systems involving ordinary differential equations (ODE) with integer-order (IO) derivatives is not sufficient, and introducing fractionalized mathematical systems to introspect various aspects of memory’s effect on leprosy dynamics becomes mandatory in this scenario.
In 2016, Westerlund pointed out that all matter has a memory [6]. The dynamic behavior of living microorganisms such as M. leprae bacteria not only depends on the conditions of their current state but also on those of their previous states to better predict and interpret the pattern of the states at some point in the future [7,8]. It is to be noted that integer-order (IO) derivatives only take into account the local properties (at time t). In the real-world explanation, the IO differentiation explores the dynamics between two different points. In such a situation, a natural question may arise about the nonlocal behavior of the two points. To solve such limitations of local differentiation, researchers like Riemann and Liouville first introduced the concept of differentiation with nonlocal or fractional-order operators [9]. A fractional (fractional-order) derivative is a generalization of the integer-order (IO) derivative and integral. It originated first in the literature about the meaning of the 1 2 -order derivative, from L’Hospital’s to Leibnitz’s, in 1965 [10] and is a very promising tool for describing the memory phenomenon [11,12,13,14,15,16]. The kernel function of a fractional derivative is called the memory function [17,18]. Now, the leprosy model proposed by Ghosh et al. [19] is an IO model which in general is memory-free. Hence, this research work was unable to reflect any kind of memory effect through the model, which expresses the limitation of the system in investigating the overall impact of the drug resistance of M. leprae bacteria on the infection and disease dissemination process of leprosy.
The memory effect can be incorporated into an ODE-based system by introducing fractional-order ( ζ ( 0 , 1 ] ) derivatives as an index of memory [20]. The significance of varying the fractional order in ( 0 , 1 ] is that ζ tends to zero, which indicates that the fractionalized system has ideal memory and that the system is free from memory as ζ approaches 1. Although most of the early research studies on fractional-order systems were based on the Caputo and Riemann–Liouville fractional-order derivatives, it has been proven that these methods have certain drawbacks. For instance, the kernels of these methods have singularity at the endpoint of an interval of definition [21,22,23], which inevitably restricts their role and development in the modelling of different natural and human-made systems. Thus, to overcome these issues, several new definitions of fractional derivatives have been introduced in recent years [24,25,26,27,28]. The basic difference among these derivatives is the varying structure of the kernels, which should be chosen to satisfy the requirements of various systems. The main difference between the Caputo–Fabrizio (CF) and Caputo fractional derivatives is that the CF derivative is obtained using an exponential decay law while the Caputo derivative is completely based on the power law [23]. For these reasons, we have chosen the CF fractional derivative over the classical Caputo derivative for the construction of a fractionalized system in our research article.
Optimal control theory is yet another effective tool used in disease modelling [29,30]. It came into existence after the formulation of the most popular Pontryagin maximum principle [31] and it provides further insight into the dynamics of the diseases in addition to more appropriate control and preventive strategies [32]. Optimal control problems that involve fractional calculus are called fractional optimal control problems (FOCPs), and these are considered to be the generalized form of classic optimal control problems. There are several research papers in the literature that provide the theoretical basis and fundamentals of FOCPs, and most of these papers extensively investigate how to formulate FOCPs and derive optimal conditions for several states and control variables using analytical and numerical methods [33,34,35,36]. Nowadays, FOCPs have been applied to epidemiological and infectious- and noninfectious-disease-based cell dynamic mathematical models in order to develop faster and more accurate methods for controlling the disease. Hence, FOCPs are potentially more flexible tools for modelling epidemiological and biological systems which are closely related to memory effects.
Although there exist various types of therapy related to leprosy, none can fully cure the disease. Moreover, due to the memory function of the bacteria, new resistant M. leprae strains have been found and drug resistance capability has been observed to be built in most of the cases [37,38]. In a recent study across three countries, it was found that from new cases, 3 % were Dapsone-resistant and 2 % were Ofloxacin-resistant, and in samples from relapsed patients, 15 % and 8 % were seen to be Dapsone-resistant and Ofloxacin-resistant, respectively [39,40,41,42].
The dynamics of the mathematical model of leprosy proposed by Ghosh et al. [19] in 2021 was based on ODE systems. To achieve our goal, we have developed a three-compartmental fractional mathematical system in this research article by replacing the integer derivative with a Caputo–Fabrizio derivative of fractional order ζ . We have proved the existence and uniqueness of solutions of the fractional model with the help of the renowned Banach Fixed-Point theorem and investigated the stability by adopting Picard’s T-stability theory. Furthermore, we have studied the mechanisms of action and efficacy of the combined biologic therapy, introducing Dapsone and Ofloxacin by stating an optimal control problem to this fractional model. Our goal is to find the accurate optimal drug dosage regimen for the eradication of the disease. The optimal control problem for minimizing the cost of immune therapy and for the simultaneous optimization of the effect of this therapy on the infected Schwann cells and M. leprae bacteria are analyzed. The main objective of this research article is to compare from a quantitative point of view the fractional system and the corresponding optimal-control-induced system with the help of the Caputo-Fabrizio fractional derivative for different values of the fractional order ζ . All of our analytical results are validated through numerical simulations using MATLAB 2016a.

2. Preliminaries

Some crucial fundamental definitions from the theory of fractional calculus are presented in this section.
Definition 1. 
The Caputo fractional derivative operator of order ζ ( ζ 0 ) & n N { 0 } is defined by [43]
D t ζ ( u ( t ) ) = 1 Γ ( n ζ ) 0 t ( t ξ ) n ζ 1 d n d t n u ( ξ ) d ξ
where n 1 ζ < n .
Definition 2. 
Let v H ( a , b ) , b > a , 0 < ζ < 1 . Then, the time-fractional Caputo–Fabrizio fractional differential operator is defined as [44]
C F D t ζ ( v ( t ) ) = M ( ζ ) 1 ζ 0 t exp ζ ( t ξ ) 1 ζ v ( ξ ) d ξ , t 0 , 0 < ζ < 1
where M ( ζ ) is a normalization function which depends on ζ and satisfies the condition M ( 0 ) = M ( 1 ) = 1 .
Definition 3. 
The Caputo–Fabrizio (CF) fractional integral operator of order 0 < ζ < 1 is given by [45]
C F I t ζ ( v ( t ) ) = 2 ( 1 ζ ) ( 2 ζ ) M ( ζ ) v ( t ) + 2 ζ ( 2 ζ ) M ( ζ ) 0 t v ( ξ ) d ξ , t 0 .
Here, it is important to note that
C F D t ζ ( v ( t ) ) = 0 if v s . is a constant function .
Furthermore, it is imperative to observe that the previous definitions completely suggest that the fractional integral of a function of order 0 < ζ < 1 is actually represented by the average of the respective functions and their integral of order one. Furthermore, the equation
2 ( 1 ζ ) ( 2 ζ ) M ( ζ ) + 2 ζ ( 2 ζ ) M ( ζ ) = 1
holds true, which provides the following formula:
M ( ζ ) = 2 ( 2 ζ ) , 0 ζ < 1 .
Here, the specific form of the normalizing function M ( ζ ) given in (5) along with the boundary conditions is used throughout the study and, more specifically, for the purpose of numerical simulations.
Definition 4. 
The Laplace transform for the CF fractional operator of order 0 < ζ 1 for k N is given as follows [44]:
L C F D t k + ζ ( v ( t ) ) ( p ) = 1 1 ζ L v k + 1 ( t ) L exp ζ 1 ζ t = p k + 1 L ( v ( t ) ) p k v ( 0 ) p k 1 v ( 0 ) v k ( 0 ) p + ζ ( 1 p ) .
To be precise, we can say that
L C F D t ζ ( v ( t ) ) ( p ) = p L ( v ( t ) ) p + ζ ( 1 p ) , k = 0 L C F D t ζ + 1 ( v ( t ) ) ( p ) = p 2 L ( v ( t ) ) p v ( 0 ) v ( 0 ) p + ζ ( 1 p ) , k = 1 .

3. The Basic Integer-Order Model and the Caputo–Fabrizio Fractionalized Mathematical Model Formulation

In recent years, fractional-order derivatives have gained huge importance in the field of modeling real-world biological phenomena. The fractional-order derivative is in fact a much generalized version of the integer-order derivative. In this research article, we now introduce the basic three-dimensional nonlinear ODE-based mathematical model developed in [19] that describes the fundamental disease dynamics of leprosy.
d S u d t = ν 1 S u 1 S u S u m a x β 1 S u B l , d S i d t = β 1 S u B l μ S i , d B l d t = ν 2 B l 1 B l B l m a x β 2 S u B l + σ S i ,
with initial values S u ( 0 ) = S u 0 0 , S i ( 0 ) = S i 0 0 and B l ( 0 ) = B 0 0 at t = 0 . Here, S u ( t ) , S i ( t ) and B l ( t ) denote the concentrations of healthy Schwann cells, infected Schwann cells and M. leprae bacteria at any time t. ν 1 and ν 2 describe the intrinsic growth rates of the S u ( t ) and B l ( t ) populations, where S u m a x and B l m a x are the carrying capacity of the same. β 1 is the rate at which healthy cells are infected by M. leprae and μ denotes the natural mortality rate of S i cells. The bacterial clearance rate results from the infection and the proliferation rates of newly produced free M. leprae bacteria, which are indicated by the parameters β 2 and σ , respectively. Modifying the above system in terms of the CF (Caputo–Fabrizio) fractional differential system of equations, we obtain
C F D t ζ S u ( t ) = ν 1 S u 1 S u S u m a x β 1 S u B l , C F D t ζ S i ( t ) = β 1 S u B l μ S i , C F D t ζ B l ( t ) = ν 2 B l 1 B l B l m a x β 2 S u B l + σ S i
with initial values S u ( 0 ) = S u 0 0 , S i ( 0 ) = S i 0 0 and B l ( 0 ) = B 0 0 at t = 0 .

3.1. The Iterative Scheme

We now consider system (7). The term S u B l is a nonlinear term and, hence, applying the Laplace transformation operator ( L ) on both sides of the system (7), we obtain that
p L S u ( t ) S u ( 0 ) p + ζ ( 1 p ) = L ν 1 S u 1 S u S u m a x β 1 S u B l ,
p L S i ( t ) S i ( 0 ) p + ζ ( 1 p ) = L β 1 S u B l μ S i ,
p L B l ( t ) B l ( 0 ) p + ζ ( 1 p ) = L ν 2 B l 1 B l B l m a x β 2 S u B l + σ S i .
The set in Equation (8) can now be rewritten in the following form:
L ( S u ( t ) ) = S u ( 0 ) p + p + ζ ( 1 p ) p L ν 1 S u 1 S u S u m a x β 1 S u B l , L ( S i ( t ) ) = S i ( 0 ) p + p + ζ ( 1 p ) p L ( β 1 S u B l μ S i ) , L ( B l ( t ) ) = B l ( 0 ) p + p + ζ ( 1 p ) p L ν 2 B l 1 B l B l m a x β 2 S u B l + σ S i .
Using the inverse Laplace, we obtain
S u ( t ) = S u ( 0 ) + L 1 p + ζ ( 1 p ) p L ν 1 S u 1 S u S u m a x β 1 S u B l , S i ( t ) = S i ( 0 ) + L 1 p + ζ ( 1 p ) p L β 1 S u B l μ S i , B l ( t ) = B l ( 0 ) + L 1 p + ζ ( 1 p ) p L ν 2 B l 1 B l B l m a x β 2 S u B l + σ S i .
We now present the series solutions generated by this method as follows:
S u = n = 0 S u n , S i = n = 0 S i n , B l = n = 0 B l n .
Furthermore, the series solution representation of the only existing nonlinear term S u B l is given as
S u B l = n = 0 G n where G n = k = 0 n S u k k = 0 n B l k k = 0 n 1 S u k k = 0 n 1 B l k .
We now use the initial conditions to achieve the following recursive formulas:
S u n + 1 = S u n ( 0 ) + L 1 p + ζ ( 1 p ) p L ν 1 S u n 1 S u n S u m a x β 1 S u n B l n , S i n + 1 = S i n ( 0 ) + L 1 p + ζ ( 1 p ) p L β 1 S u n B l n μ S i n , B l n + 1 = B l n ( 0 ) + L 1 p + ζ ( 1 p ) p L ν 2 B l n 1 B l n B l m a x β 2 S u n B l n + σ S i n .
The approximate solution is assumed to be obtained as a limit when n , i.e., S u ( t ) = lim n S u n ( t ) , S i ( t ) = lim n S i n ( t ) and B l ( t ) = lim n B l n ( t ) .

3.2. Stability Analysis

In this section, first, we present the detailed definition of the T-stability of Picard’s iteration [46].
Definition 5. 
Suppose T is a self-map on a complete metric space ( Y , d ) . Consider an iteration y n + 1 = g ( T , y n ) . Furthermore, let us assume that P ( T ) is the fixed-point set of T with P ( T ) ϕ and let { y n } converge to some point y P ( T ) . Let { z n } Y and define { u n } = d ( Z n + 1 , g ( T , z n ) ) . Now, if u n 0 implies that z n y , then the iteration method y n + 1 = g ( T , y n ) is said to be T-stable. Furthermore, note that the convergence of { z n } guarantees that { z n } must be bounded above. If all these conditions hold true for y n + 1 = g ( T , y n ) , then Picard’s iteration method is called T-stable.
Let ( X , . ) be a Banach space. As every Banach space is a complete metric space with the metric induced by the associated norm, Definition 5 holds true for ( X , . ) also.
Theorem 1. 
Let T be a self-map on the space ( X , . ) , which satisfies the following:
T x T y     Λ x T x   +   ϱ x y for all x , y X
where Λ 0 and ϱ [ 0 , 1 ) . Suppose T has a fixed point. Then, T is Picard’s T-stable.
Now, let us define a self-map T as
T ( S u n ( t ) ) = S u n + 1 = S u n ( 0 ) + L 1 p + ζ ( 1 p ) p L ν 1 S u n 1 S u n S u m a x β 1 S u n B l n , T ( S i n ( t ) ) = S i n + 1 = S i n ( 0 ) + L 1 p + ζ ( 1 p ) p L β 1 S u n B l n μ S i n , T ( B l n ( t ) ) = B l n + 1 = B l n ( 0 ) + L 1 p + ζ ( 1 p ) p L ν 2 B l n 1 B l n B l m a x β 2 S u n B l n + σ S i n
For all m , n N , let us first construct the following differences:
T ( S u n ( t ) ) T ( S u m ( t ) ) = S u n ( t ) S u m ( t ) + L 1 p + ζ ( 1 p ) p L ν 1 S u n 1 S u n S u m a x β 1 S u n B l n L 1 p + ζ ( 1 p ) p L ν 1 S u m 1 S u m S u m a x β 1 S u m B l m , T ( S i n ( t ) ) T ( S i m ( t ) ) = S i n ( t ) S i m ( t ) + L 1 p + ζ ( 1 p ) p L β 1 S u n B l n μ S i n L 1 p + ζ ( 1 p ) p L β 1 S u m B l m μ S i m , T ( B l n ( t ) ) T ( B l m ( t ) ) = B l n ( t ) B l m ( t ) + L 1 p + ζ ( 1 p ) p L ν 2 B l n 1 B l n B l m a x β 2 S u n B l n + σ S i n L 1 p + ζ ( 1 p ) p L ν 2 B l m 1 B l m B l m a x β 2 S u m B l m + σ S i m
where p + ζ ( 1 p ) p is a Lagrange multiplier in fractional form. As all the solutions S u n , S i n , B l n are Cauchy sequences in the Banach space ( X , . ) , it is true that S u n S u m 0 , S i n S i m 0 and B l n B l m 0 as n , m . Due to this similar behavior of the solutions, i.e., comparative influence of the solutions in [47], we have
S u n ( t ) S u m ( t )     S i n ( t ) S i m ( t ) , S u n ( t ) S u m ( t )     B l n ( t ) B l m ( t ) .
Now, applying the norm on both sides of the first equation of (16), we obtain
T ( S u n ( t ) ) T ( S u m ( t ) ) =   S u n ( t ) S u m ( t ) +   L 1 p + ζ ( 1 p ) p L ν 1 S u n 1 S u n S u m a x β 1 S u n B l n   L 1 p + ζ ( 1 p ) p L ν 1 S u m 1 S u m S u m a x β 1 S u m B l m =   S u n ( t ) S u m ( t ) + L 1 p + ζ ( 1 p ) p L ν 1 ( S u n ( t ) S u m ( t ) +   ν 1 S u m a x S u n ( S u n S u m ) + ν 1 S u m a x S u m ( S u n S u m ) +   β 1 B l n ( S u n S u m ) + β 1 S u m ( B l n B l m ) ] ] .
Using triangle inequality, we obtain
T ( S u n ( t ) ) T ( S u m ( t ) )   S u n ( t ) S u m ( t ) + L 1 p + ζ ( 1 p ) p L ν 1 ( S u n ( t ) S u m ( t ) + ν 1 S u m a x S u n ( S u n S u m ) + ν 1 S u m a x S u m ( S u n S u m ) + β 1 B l n ( S u n S u m ) + β 1 S u m ( B l n B l m ) ] ] .
Then, using the relations in (17), we obtain
T ( S u n ( t ) ) T ( S u m ( t ) )   S u n ( t ) S u m ( t ) + L 1 p + ζ ( 1 p ) p L ν 1 ( S u n ( t ) S u m ( t ) + ν 1 S u m a x S u n ( S u n S u m ) + ν 1 S u m a x S u m ( S u n S u m ) + β 1 B l n ( S u n S u m ) + β 1 S u m ( S u n S u m ) ] ]   S u n ( t ) S u m ( t ) [ 1 + ν 1 E 1 ( ζ ) 2 M 1 ν 1 S u m a x E 2 ( ζ )   β 1 ( M 1 + M 3 ) E 3 ( ζ ) ]
where E 1 ( ζ ) , E 2 ( ζ ) and E 3 ( ζ ) are functions of L 1 p + ζ ( 1 p ) p L ( . ) and S u n   <   M 1 , S i n   <   M 2 and B l n   <   M 3 . Proceeding similarly, we obtain from the second and third equations of (16)
T ( S i n ( t ) ) T ( S i m ( t ) )   S i n ( t ) S i m ( t ) 1 + β 1 ( M 1 + M 3 ) E 3 ( ζ ) μ E 4 ( ζ )
and
T ( B l n ( t ) ) T ( B l m ( t ) )   B l n ( t ) B l m ( t ) [ 1 + ν 2 E 5 ( ζ ) 2 M 3 ν 2 B l m a x E 6 ( ζ ) β 2 ( M 1 + M 3 ) E 3 ( ζ ) + σ E 7 ( ζ ) ]
where E 4 ( ζ ) , E 5 ( ζ ) , E 6 ( ζ ) and E 7 ( ζ ) are functions of L 1 p + ζ ( 1 p ) p L ( . ) and
1 + ν 1 E 1 ( ζ ) 2 M 1 ν 1 S u m a x E 2 ( ζ ) β 1 ( M 1 + M 3 ) E 3 ( ζ ) < 1 , 1 + β 1 ( M 1 + M 3 ) E 3 ( ζ ) μ E 4 ( ζ ) < 1 , 1 + ν 2 E 5 ( ζ ) 2 M 3 ν 2 B l m a x E 6 ( ζ ) β 2 ( M 1 + M 3 ) E 3 ( ζ ) + σ E 7 ( ζ ) < 1 .
So, we can conclude that the self-map T defined in (15) has a fixed point. In view of (21) and also choosing ϱ = ( 0 , 0 , 0 ) and
Λ = 1 + ν 1 E 1 ( ζ ) 2 M 1 ν 1 S u m a x E 2 ( ζ ) β 1 ( M 1 + M 3 ) E 3 ( ζ ) , 1 + β 1 ( M 1 + M 3 ) E 3 ( ζ ) μ E 4 ( ζ ) , 1 + ν 2 E 5 ( ζ ) 2 M 3 ν 2 B l m a x E 6 ( ζ ) β 2 ( M 1 + M 3 ) E 3 ( ζ ) + σ E 7 ( ζ ) ,
we can see that all the conditions of Theorem 1 are satisfied. Thus, the self-mapping T is Picard’s T-stable. Here, it is important to note that Λ is a constant, not a function.
Summarizing the previous discussions, we now present the following theorem.
Theorem 2. 
Consider system (7) with the set of equations in system (14). Let T be a self-map as defined by (15). If the conditions (21) are satisfied by T, then T has a fixed point and, hence, T is Picard’s T-stable.

3.3. Existence of the Solutions

Using fixed-point theory, we now show the existence of the solutions of system (7) in this subsection. For this, let us first observe that
S u ( t ) S u 0 ( t ) = 2 ( 1 ζ ) M ( ζ ) ( 2 ζ ) ν 1 S u ( t ) 1 S u ( t ) S u m a x β 1 S u ( t ) B l ( t ) + 2 ζ M ( ζ ) ( 2 ζ ) 0 t ν 1 S u ( y ) 1 S u ( y ) S u m a x β 1 S u ( y ) B l ( y ) d y , S i ( t ) S i 0 ( t ) = 2 ( 1 ζ ) M ( ζ ) ( 2 ζ ) β 1 S u B l μ S i + 2 ζ M ( ζ ) ( 2 ζ ) 0 t β 1 S u ( y ) B l ( y ) μ S i ( y ) d y , B l ( t ) B l 0 ( t ) = 2 ( 1 ζ ) M ( ζ ) ( 2 ζ ) ν 2 B l ( t ) 1 B l ( t ) B l m a x β 2 S u ( t ) B l ( t ) + σ S i + 2 ζ M ( ζ ) ( 2 ζ ) 0 t ν 2 B l ( y ) 1 B l ( y ) B l m a x β 2 S u ( y ) B l ( y ) + σ S i ( y ) d y .
Let T 1 be an operator on H to itself, i.e., T 1 : H H . Here, T 1 is chosen as an operator for the entire system. Applying it, we obtain that
T 1 ( S u ( t ) ) = 2 ( 1 ζ ) M ( ζ ) ( 2 ζ ) K 1 ( t , S u ( t ) ) + 2 ζ M ( ζ ) ( 2 ζ ) 0 t K 1 ( y , S u ( y ) ) d y , T 1 ( S i ( t ) ) = 2 ( 1 ζ ) M ( ζ ) ( 2 ζ ) K 2 ( t , S i ( t ) ) + 2 ζ M ( ζ ) ( 2 ζ ) 0 t K 2 ( y , S i ( y ) ) d y , T 1 ( B l ( t ) ) = 2 ( 1 ζ ) M ( ζ ) ( 2 ζ ) K 3 ( t , B l ( t ) ) + 2 ζ M ( ζ ) ( 2 ζ ) 0 t K 3 ( y , B l ( y ) ) d y
where
K 1 ( t , S u ( t ) ) = ν 1 S u ( t ) 1 S u ( t ) S u m a x β 1 S u ( t ) B l ( t ) , K 2 ( t , S i ( t ) ) = β 1 S u ( t ) B l ( t ) μ S i ( t ) , K 3 ( t , B l ( t ) ) = ν 2 B l ( t ) 1 B l ( t ) B l m a x β 2 S u ( t ) B l ( t ) + σ S i ( t ) .
Let P H be bounded. We aim to show that T 1 ( P ) ¯ is compact to ensure the existence and boundedness of the solutions of system (7), where T 1 is defined as in (22). We can see that there exist positive reals κ 1 , κ 2 and κ 3 such that S u   <   κ 1 , S i   <   κ 2 and B l   <   κ 3 . From the first equation of (22), we can write
T 1 ( S u ( t ) ) =   2 ( 1 ζ ) M ( ζ ) ( 2 ζ ) K 1 ( t , S u ( t ) ) + 2 ζ M ( ζ ) ( 2 ζ ) 0 t K 1 ( y , S u ( y ) ) d y 2 ( 1 ζ ) M ( ζ ) ( 2 ζ ) K 1 ( t , S u ( t ) ) + 2 ζ M ( ζ ) ( 2 ζ ) 0 t K 1 ( y , S u ( y ) ) d y 2 ( 1 ζ ) M ( ζ ) ( 2 ζ ) + a 1 2 ζ M ( ζ ) ( 2 ζ ) K 1 ( t , S u ( t ) ) R 1 2 ( 1 ζ ) M ( ζ ) ( 2 ζ ) + a 1 2 ζ M ( ζ ) ( 2 ζ )
which implies
T 1 ( S u ( t ) )   2 R 1 M ( ζ ) ( 2 ζ ) ( 1 + ζ a 1 ζ )
and also proceeding similarly, we can obtain
T 1 ( S i ( t ) )   2 R 2 M ( ζ ) ( 2 ζ ) ( 1 + ζ a 2 ζ ) , T 1 ( B l ( t ) )   2 R 3 M ( ζ ) ( 2 ζ ) ( 1 + ζ a 3 ζ )
where
R 1 = max t [ 0 , 1 ] S u [ 0 , κ 1 ] K 1 ( t , S u ( t ) ) ,
R 2 = max t [ 0 , 1 ] S i [ 0 , κ 2 ] K 2 ( t , S i ( t ) ) ,
R 3 = max t [ 0 , 1 ] B l [ 0 , κ 3 ] K 3 ( t , B l ( t ) ) .
Hence, we have proved that T 1 ( P ) is bounded. Let, t 2 > t 1 and S u , S i , B l P . So, for a given ϵ > 0 , there exists η satisfying that ( t 2 t 1 )   <   η , and we can write the following:
K 1 ( t 2 , S u ( t 2 ) ) K 1 ( t 1 , S u ( t 1 ) )   ν 1 S u ( t 2 ) S u ( t 1 ) +   ν 1 S u m a x S u ( t 2 ) + S u ( t 1 ) S u ( t 2 ) S u ( t 1 ) +   β 1 B l S u ( t 2 ) S u ( t 1 )   ν 1 S u ( t 2 ) S u ( t 1 ) + 2 κ 1 ν 1 S u m a x S u ( t 2 ) S u ( t 1 ) +   β 1 κ 3 S u ( t 2 ) S u ( t 1 )   ν 1 + 2 κ 1 ν 1 S u m a x + β 1 κ 3 S u ( t 2 ) S u ( t 1 ) .
Assuming that if the function S u ( t ) is Lipschitz-continuous, i.e., for some real number L 1 0 and for all t 1 , t 2 , the inequality S u ( t 2 ) S u ( t 1 )   L 1 t 2 t 1 holds, we can rewrite (23) as
K 1 ( t 2 , S u ( t 2 ) ) K 1 ( t 1 , S u ( t 1 ) ) G 1 t 2 t 1
where G 1 = L 1 ν 1 + 2 κ 1 ν 1 S u m a x + β 1 κ 3 . Similarly, we have
K 2 ( t 2 , S i ( t 2 ) ) K 2 ( t 1 , S i ( t 1 ) ) G 2 t 2 t 1 ,
K 3 ( t 2 , B l ( t 2 ) ) K 3 ( t 1 , B l ( t 1 ) ) G 3 t 2 t 1
if S i ( t ) and B l ( t ) are also Lipschitz-continuous, i.e., for some real numbers L 2 , L 3 0 , the conditions
S i ( t 2 ) S i ( t 1 ) L 2 t 2 t 1 , B l ( t 2 ) B l ( t 1 ) L 3 t 2 t 1 hold , respectively , for all t 1 , t 2 .
Furthermore,
T 1 ( S u ( t 2 ) ) T 1 ( S u ( t 1 ) ) 2 ( 1 ζ ) M ( ζ ) ( 2 ζ ) K 1 ( t 2 , S u ( t 2 ) ) K 1 ( t 1 , S u ( t 1 ) ) + 2 ζ M ( ζ ) ( 2 ζ ) R 1 K 1 ( t 2 , S u ( t 2 ) ) K 1 ( t 1 , S u ( t 1 ) ) 2 ( 1 ζ ) M ( ζ ) ( 2 ζ ) G 1 t 2 t 1 + 2 ζ M ( ζ ) ( 2 ζ ) R 1 G 1 t 2 t 1 ( using inequality ( 24 ) ) .
Finally, choosing
η = ϵ 2 ( 1 ζ ) M ( ζ ) ( 2 ζ ) G 1 + 2 ζ M ( ζ ) ( 2 ζ ) R 1 G 1 ,
we can see that T 1 ( S u ( t 2 ) ) T 1 ( S u ( t 1 ) ) ϵ holds.
Similarly proceeding and using inequalities (25) and (26), we can also easily show that T 1 ( S i ( t 2 ) ) T 1 ( S i ( t 1 ) ) ϵ and T 1 ( B l ( t 2 ) ) T 1 ( B l ( t 1 ) ) ϵ hold, which guarantees the equicontinuity of T 1 . Hence, according to the well-known Arzela–Ascoli theorem, we can say that T 1 ( P ) ¯ is compact. Next, we present the following theorem by summarising the previous discussions on the existence of the solutions of system (7), and then we move forward to show the uniqueness of the solutions of system (7).
Theorem 3. 
Let P H be bounded and T 1 be defined as in (22). Then, there exist κ 1 , κ 2 and κ 3 such that if the functions S u ( t ) , S i ( t ) and B l ( t ) are Lipschitz-continuous, i.e., if for some real numbers, L 1 , L 2 and L 3 0 , the following conditions hold
S u ( t 2 ) S u ( t 1 )   L 1 t 2 t 1 , S i ( t 2 ) S i ( t 1 )   L 2 t 2 t 1 , B l ( t 2 ) B l ( t 1 )   L 3 t 2 t 1 ,
for all t 1 , t 2 , then T 1 ( P ) ¯ is compact. Thus, all the solutions of system (7) exist and are bounded.

3.4. Uniqueness of the Solutions

To prove the uniqueness of the solutions of system (7), let us consider the mapping T 1 again which was defined previously. Now,
T 1 ( S u ( t ) ) T 1 ( S u ˜ ( t ) ) = 2 ( 1 ζ ) M ( ζ ) ( 2 ζ ) ( K 1 ( t , S u ( t ) ) K 1 ( t , S u ˜ ( t ) ) ) + 2 ζ M ( ζ ) ( 2 ζ ) 0 t ( K 1 ( y , S u ( y ) ) K 1 ( y , S u ˜ ( y ) ) ) 2 D 1 M ( ζ ) ( 2 ζ ) S u ( t ) S u ˜ ( t ) .
Similarly, we can obtain
T 1 ( S i ( t ) ) T 1 ( S i ˜ ( t ) ) 2 D 2 M ( ζ ) ( 2 ζ ) S i ( t ) S i ˜ ( t ) , T 1 ( B l ( t ) ) T 1 ( B l ˜ ( t ) ) 2 D 3 M ( ζ ) ( 2 ζ ) B l ( t ) B l ˜ ( t )
where D 1 , D 2 , D 3 R . Hence, the operator T 1 is a contraction if the following conditions hold:
2 D 1 M ( ζ ) ( 2 ζ ) S u ( t ) S u ˜ ( t )   <   1 , 2 D 2 M ( ζ ) ( 2 ζ ) S i ( t ) S i ˜ ( t )   <   1 , 2 D 3 M ( ζ ) ( 2 ζ ) B l ( t ) B l ˜ ( t )   <   1
which ensures the existence of unique positive solutions of system (7) using fixed-point theorem.

4. Equilibria and Stability

Our CF fractionalized mathematical model (7) has two equilibria, namely the disease-free equilibrium E 0 and the endemic equilibrium E * . Here, E 0 is given as E 0 = ( S u m a x , 0 , 0 ) . The value of the basic reproduction number R 0 is given as R 0 = β 1 σ S u m a x μ ( β 2 S u m a x ν 2 ) [19]. R 0 is actually interpreted as the secondary number of new infections in a completely susceptible healthy Schwann cell population and, based on the above, we now present the following theorem on the stability of E 0 for our system (7) as follows:
Theorem 4. 
The disease-free equilibrium E 0 of system (7) is locally asymptotically stable if R 0 < 1 .
To obtain the coordinates of the endemic equilibrium E * , we now set the right-hand sides of system (7) to zero. Hence, we obtain the values of S u * , S i * and B l * , which are already described in the article by Ghosh et al. [19]. In this context, we now present the following theorem, which describes the required criterion about the stability of E * [48].
Theorem 5. 
If the matrix ( I ( 1 ζ ) J ) is invertible, then the endemic equilibrium E * of the CF fractionalized system (7) is locally asymptotically stable if all the roots of the characteristic equation d e t ( x ( I ( 1 ζ ) J ) ζ J ) = 0 of system (7) evaluated at E * are negative real or have negative real parts where J denotes the Jacobian matrix of system (7) at E * = ( S u * , S i * , B l * ) .

5. Optimal-Control-Induced Caputo–Fabrizio Fractional Mathematical Model

Optimal control is a very useful tool for controlling the progression of a disease in the human body. Furthermore, this tool has gained major importance lately for the investigation of efficient and cost-effective drug treatment policies for various infectious diseases and for other different important biological problems based on fractional mathematical models [36,43,49,50]. In our research article, we have analyzed our formulated CF fractionalized model (7) by incorporating two control functions, u 1 ( t ) and u 2 ( t ) ; one is an effect of the drug Ofloxacin and another is of Dapsone on various cell densities, respectively. Here, Ofloxacin blocks the occurrence of new infections, and by preventing the formation of folic acid, Dapsone specifically inhibits replication of M. leprae bacteria. The optimal-control-induced CF fractional mathematical model is presented as
C F D t ζ S u ( t ) = ν 1 S u 1 S u S u m a x β 1 ( 1 u 1 ( t ) ) S u B l , C F D t ζ S i ( t ) = β 1 ( 1 u 1 ( t ) ) S u B l μ S i , C F D t ζ B l ( t ) = ν 2 ( 1 u 2 ( t ) ) B l 1 B l B l m a x β 2 ( 1 u 1 ( t ) ) S u B l + σ S i
with initial values S u ( 0 ) = S u 0 0 , S i ( 0 ) = S i 0 0 and B l ( 0 ) = B 0 0 at t = 0 .
Our main aim is to decrease the number of infected cells and the bacterial density as well as to increase the healthy cell concentrations. Let us now consider the state system given by (27) with the class of admissible controls defined as
U = { ( u 1 ( · ) , u 2 ( · ) ) : u 1 , u 2 are Lebesgue measurable functions on [ 0 , 1 ] and 0 u i ( t ) 1 for i = 1 , 2 } .
So, the objective function for the CF fractionalized optimal control system (27), (28) is given as
J ( u 1 ( · ) , u 2 ( · ) ) = 0 t f 1 2 C 1 u 1 2 ( t ) + 1 2 C 2 u 2 2 ( t ) + S i 2 ( t ) + B l 2 ( t ) d t
where C 1 and C 2 measure the cost associated with the control functions u 1 ( t ) and u 2 ( t ) , respectively. Then, we find the optimal controls u 1 and u 2 to minimize the cost function
J ( u 1 , u 2 ) = 0 t f ψ ( S u ( t ) , S i ( t ) , B l ( t ) , u 1 ( t ) , u 2 ( t ) , t ) d t
subject to the constraints 0 C F D t ζ ( S u ( t ) ) = α 1 , 0 C F D t ζ ( S i ( t ) ) = α 2 and 0 C F D t ζ ( B l ( t ) ) = α 3 , where α j = α j ( S u , S i , B l , u 1 , u 2 , t ) and j = 1 , 2 , 3 , and the given initial conditions are S u ( 0 ) = S u 0 and S i ( 0 ) = S i 0 , B l ( 0 ) = B l 0 .
Now, we first present a formulation of a generalized fractional optimal control problem (FOCP) and deduce the necessary conditions for its optimality. For this, let us consider a generalized FOCP as
J ( v ) = 0 t f L ( t , x , v ) d t d t
subject to the constraints 0 C F D t ζ ( x ( t ) ) = g ( t , x , v ) with initial condition x ( 0 ) = x 0 . Here, x ( t ) and v ( t ) are state and control vectors, respectively, and L and g are differentiable functions with 0 < ζ 1 .
Theorem 6. 
We define a Hamiltonian as follows:
H ( t , x , v , ϕ ) = L ( t , x , v ) + ϕ g ( t , x , v )
where ϕ C 1 [ 0 , t f ] is a function. If ϕ , x , v satisfy the following equations:
0 C F D t ζ ( x ( t ) ) = H ( t , x ( t ) , v ( t ) , ϕ ( t ) ) ϕ , t C F D t f ζ ( ϕ ( t ) ) = H ( t , x ( t ) , v ( t ) , ϕ ( t ) ) x , H ( t , x ( t ) , v ( t ) , ϕ ( t ) ) v = 0 , ϕ ( t f ) = 0
then ( x , v ) is the minimizer of (31).
Proof. 
Let us first substitute Equation (32) in (31); hence, we obtain
J ( v ) = 0 t f H ( t , x , v , ϕ ) ϕ * g ( t , x , v ) d t .
The necessary condition for the optimality of the FOCP is
δ J ( v ) = 0 .
To obtain the optimal control laws, we choose the variation of equation (33) as
δ J ( v ) = δ x H x + δ v s . H v + δ ϕ H ϕ δ ϕ 0 C F D t ζ ( x ( t ) ) ϕ 0 C F D t ζ ( δ x ( t ) ) d t
where δ x , δ v , δ ϕ are variations of x, v and ϕ , respectively.
Again,
0 t f ϕ ( t ) 0 C F D t ζ ( δ x ( t ) ) d t = 0 t f δ x 0 C F D t ζ ( ϕ ( t ) ) d t t C F I t f 1 ζ ( ϕ ( t ) δ x .
Substituting (36) in (35), we obtain
δ J ( v ) = 0 t f δ ( x ) H x 0 C F D t ζ ( ϕ ( t ) ) + δ v s . H v + δ ϕ H ϕ 0 C F D t ζ ( x ( t ) ) d t + t C F I t f 1 ζ ( ϕ ( t ) ) δ x | t = t f .
Now, we know δ J ( v ) = 0 . Hence, considering Equation (37), the coefficients of δ x , δ v , δ ϕ must be equal to zero, which leads us to the following equations:
0 C F D t ζ ( x ( t ) ) = H ( t , x ( t ) , v ( t ) , ϕ ( t ) ) ϕ , t C F D t f ζ ( ϕ ( t ) ) = H ( t , x ( t ) , v ( t ) , ϕ ( t ) ) x , H ( t , x ( t ) , v ( t ) , ϕ ( t ) ) v = 0 , t C F I t f 1 ζ ( ϕ ( t ) ) | t = t f = ϕ ( t f ) = 0 .
In addition, the following necessary conditions must also hold for the optimality of the FOCP defined in (31), which are noted here in the form of the following lemma.
Lemma 1. 
The following conditions hold true for the generalized FOCP described in (31):
t C F D t f ζ ( ϕ ( t ) ) = H ( t , x ( t ) , v ( t ) , ϕ ( t ) ) x
and 0 C F D t ζ ( ϕ ( t f t ) ) = H ( t f t , x ( t f t ) , v ( t f t ) , ϕ ( t f t ) ) x
where 0 < ζ 1 .
Proof. 
The definition of the CF fractional derivative (2) is given as
0 C F D t ζ ( f ( t ) ) = 1 1 ζ 0 t exp ζ ( t x ) 1 ζ f ( x ) d x , t 0 , 0 < ζ < 1 .
From the above-mentioned definition of the CF fractional derivative, it follows that
t f t C F D t f ζ ( ϕ ( t f t ) ) = 1 1 ζ t f t t f exp ζ ( t f t x ) 1 ζ ϕ ( x ) d x .
Now, assuming t f x = y , from Equation (40), we obtain that
t f t C F D t f ζ ( ϕ ( t f t ) ) = 1 1 ζ t 0 exp ζ ( y t ) 1 ζ ϕ ( t f y ) ( d y ) = 1 1 ζ 0 t exp ζ ( y t ) 1 ζ ϕ ( t f y ) d y .
Hence, the optimality conditions are achieved as
0 C F D t ζ ( x ( t ) ) = H ( t , x ( t ) , v ( t ) , ϕ ( t ) ) ϕ , 0 C F D t ζ ( ϕ ( t f t ) ) = H ( t f t , x ( t f t ) , v ( t f t ) , ϕ ( t f t ) ) x and H v = 0 where H : = H ( t , x ( t ) , v ( t ) , ϕ ( t ) ) .
We now shift our focus from the generalized point of view, specifically to our CF fractionalized optimal-control-induced (CFOC) system (27), (29), (30). Let us first consider the following modified cost function:
Z ^ = 0 T H ( S u , S i , B l , u 1 , u 2 , t ) j = 1 3 θ j α j ( S u , S i , B l , u 1 , u 2 , t ) d t .
Hence, the Hamiltonian is defined as
H ( S u , S i , B l , u 1 , u 2 , t ) = ψ ( S u , S i , B l , u 1 , u 2 , t ) + j = 1 3 θ j α j ( S u , S i , B l , u 1 , u 2 , t ) .
Then, utilizing (42), the necessary and sufficient conditions for the CF fractional optimal control (CFOC) problem defined in (27), (29), (30) are given as
0 C F D ζ θ 1 = H S u , 0 C F D ζ θ 2 = H S i , 0 C F D ζ θ 3 = H B l , H u i = 0 , i = 1 , 2 , and 0 C F D t ζ ( S u ( t ) ) = H θ 1 , 0 C F D t ζ ( S i ( t ) ) = H θ 2 , . 0 C F D t ζ ( B l ( t ) ) = H θ 3 .
Moreover, θ 1 , θ 2 and θ 3 are Lagrange’s multipliers, which express the necessary and sufficient conditions in terms of the Hamiltonian for the fractional optimal control problem defined above.
Now, consider system (27). Let us consider the Hamiltonian defined in (44). Rewriting it in the following form, we obtain
H ( S u , S i , B l , u 1 , u 2 , θ ) = 1 2 C 1 u 1 2 + 1 2 C 2 u 2 2 + S i 2 + B l 2 + θ 1 ν 1 S u 1 S u S u m a x β 1 ( 1 u 1 ( t ) ) S u B l + θ 2 β 1 ( 1 u 1 ( t ) ) S u B l μ S i + θ 3 ν 2 ( 1 u 2 ( t ) ) B l 1 B l B l m a x β 2 ( 1 u 1 ( t ) ) S u B l + σ S i .
Theorem 7. 
If u 1 * , u 2 * are optimal controls of the given CFOC system defined by (27), (29), (30), and if S u * , S i * , B l * are the corresponding optimal paths, then there exist co-state variables θ 1 * , θ 2 * , θ 3 * such that besides the given control system being satisfied, the following conditions, i.e., the co-state equations, hold true also. The co-state equations are given as
0 C F D t ζ θ 1 * = θ 1 * ν 1 2 S u ν 1 S u m a x β 1 ( 1 u 1 ) B l θ 2 * β 1 ( 1 u 1 ) B l + θ 3 * β 2 ( 1 u 1 ) B l , 0 C F D t ζ θ 2 * = θ 2 * μ θ 3 * σ 2 S i , 0 C F D t ζ θ 3 * = θ 1 * β 1 ( 1 u 1 ) S u θ 2 * β 1 ( 1 u 1 ) S u θ 3 * ν 2 ( 1 u 2 ) 1 2 B l B l m a x β 2 ( 1 u 1 ) S u 2 B l
with transversality conditions θ 1 * ( t f ) = 0 , θ 2 * ( t f ) = 0 , θ 3 * ( t f ) = 0 , and the optimality conditions are given by
H ( S u * ( t ) , S i * ( t ) , B l * ( t ) , u 1 * ( t ) , u 2 * ( t ) ) = m i n 0 u i 1 H ( S u * , S i * , B l * , u i * ) , u 1 * ( t ) = min 1 , max 0 , S u B l ( β 1 θ 2 * β 2 θ 3 * β 1 θ 1 * ) C 1 , u 2 * ( t ) = min 1 , max 0 , θ 3 * ν 2 B l 1 B l B l m a x C 2 .
Proof. 
The adjoint system (46) is obtained from H as follows
d θ 1 d t = H S u , d θ 2 d t = H S i , d θ 3 d t = H B l
with transversality conditions given as θ 1 ( t f ) = θ 2 ( t f ) = θ 3 ( t f ) = 0 . Furthermore, the characterization of the CF fractionalized optimal controls u 1 * ( t ) and u 2 * ( t ) are achieved by solving the following equations
H u 1 = 0 , H u 2 = 0
on the interior of the control set and utilizing the properties of the control space U . □
The analytical sections of our study come to an end here and next we proceed to the numerical outcomes for validation of the analytical portions of our proposed systems.

6. Numerical Simulation

In this section, we perform numerical simulations for both the Caputo–Fabrizo fractional system denoted by (7) and also the CF fractionalized optimal control system (27). All the numerical results are compared with the analytical and theoretical outcomes previously achieved. We chose the initial values according to the cardinal rule of scientific hypothesis. Some values of the parameters were chosen from the research articles [19,51,52], and the other values were estimated. The values of the parameters which we have used here are described in the following table denoted by Table 1. All of our numerical findings here were obtained using MATLAB 2016a. Throughout the study, the interval of consideration was chosen as [ 0 , 2.5 × 10 3 ] .
Here, we want to mention that Table 1 actually refers to the values of the system parameters for system (7) and system (27) for the fractional order ζ = 1 . During the simulations of the figures for ζ ( 0 , 1 ) , we adopted the technique proposed by Atangana et al. [53]. To avoid the dimension mismatching of ( t i m e ) 1 and ( t i m e ) ( ζ ) between the left- and right-hand sides of the systems, the dimensions of the system parameters were modified accordingly and the corresponding values were utilized for the numerical simulations.
In Figure 1 and Figure 2, the behavior of the cell densities and the phase portrait diagrams of the healthy Schwann cells, infected Schwann cells and M. leprae bacteria for system (7) in the 3-D phase space are exhibited, respectively, for ζ = 1 , i.e., for the classical integer-order system. Figure 1 depicts oscillatory periodic solutions and stable limit cycles. On the contrary, in Figure 2, stable solutions are observed whenever the value of ν 2 is increased from 0.03 to 0.05 , which describes that the intrinsic growth rate of the bacterial population is a very crucial parameter for demonstrating the dynamical shift in system (7).
We now move on to understand the behavior of system (7) in the previous memory states, i.e., for the noninteger cases or fractional-order cases for the value of ζ ( 0 , 1 ) . Considering ζ = 0.8 , Figure 3 shows that system (7) produces oscillatory solutions and stable limit cycles for ν 2 = 0.03 , while S u cell, S i cell and B l bacteria reach the stable concentrations of approximately 50 mm 3 , 120 mm 3 and 1100 mm 3 , respectively, at the endemic steady state E * described in Figure 4.
Next, a comparison of system (7)’s behavior is investigated for two different values of ζ . Keeping the whole parameter set fixed, we varied the fractional order from ζ = 1 to ζ = 0.6 . The numerical outcomes in Figure 5 exhibit unstable behavior with sustained oscillations of the cell population densities of system (7) for ζ = 1 , but moving towards the previous memory state for ζ = 0.6 , Figure 6 shows that after a little initial fluctuation, the system state populations become asymptotically stable. In addition, in Figure 7, the stability regions of the interior equilibrium E * ( S u * , S i * , B l * ) for the CF fractionalized system (7) for three different values of ζ , i.e., for ζ = 1 , 0.8 , 0.6 , are clearly demonstrated in Figure 7a, Figure 7b and Figure 7c, respectively.
In Figure 8 and Figure 9, the effect of the optimal control treatment policy has been demonstrated on the Caputo–Fabrizio fractionalized optimal control (CFOC) system for the fractional orders ζ = 0.9 and ζ = 0.6 , respectively. In both cases, the concentrations of healthy Schwann cells ( S u ) are observed to be increased, and also the densities of S i and B l are decreased in the body of a leprosy-infected person. The bacterial concentration B l decreases to 510 mm 3 in Figure 8c, while it decreases to a stable concentration of 420 mm 3 depicted in Figure 9c. This indicates that the CFOC system (27) acts better when more previous memory states are considered for ζ = 0.6 . Now, looking into the optimal control profiles of u 1 * and u 2 * , we can see that the drug therapy Dapsone denoted by u 2 * ( t ) needs to be increased after 520 days in case of ζ = 0.9 , while it should be increased after nearly 1100 days up to the range 0.8–0.9 for ζ = 0.6 described in Figure 8d and Figure 9d, respectively. This happens due to the memory effect of M. leprae-induced infection and as the previous memory stages of the bacteria are extremely correlated with the drug resistance scenarios during the treatment [7]. Indeed, the M. leprae bacteria are highly Dapsone-resistant [39,54]. To tackle the dissemination of leprosy into the human body and to effectively inhibit bacterial drug resistance, applying the optimal control treatment approach for the previous memory state of ζ = 0.6 is more realistic in nature and appropriate than the present state or the states very adjacent to ζ = 0.9 .

7. Discussion and Conclusions

Leprosy is a chronic infectious bacterial disease caused by the bacterium M. leprae. Despite several attempts by clinical and experimental scientists in the past, complete elimination of leprosy has still not been attained. The drug resistance characteristics of the bacterium, its associated relapse problems, and designing a perfect drug dose regimen in a cost-effective way remain a great challenge [37,55,56] for researchers worldwide. Applications of new innovative tools and techniques are essential in this situation. Keeping this in mind, we have formulated and analyzed a three-dimensional CF fractional-order-based mathematical system and, most importantly, investigated the impacts of memory effects on the CF fractionalized optimal control system by incorporating a combined drug therapy of Ofloxacin and Dapsone. In this section, we will discuss the crucial interpretations obtained from the analytical and numerical outcomes for both system (7) and system (27).
Firstly, we formulated an iterative scheme using the Laplace and inverse Laplace transformations in Section 3.1. Following that, we established the stability of the solutions using Picard’s T-stability iterative criterion in Theorem 2 in Section 3.2. For demonstrating the existence and uniqueness of the solutions (Theorem 3 in Section 3.3 and Section 3.4) of system (7), we used the well-known Banach fixed-point theorem and Arzela–Ascoli theorem. The formula for the basic reproduction number R 0 was derived and the local asymptotic stability of E 0 for the CF fractionalized system (7) was investigated with respect to the threshold value of R 0 = 1 . Besides this, we also described the stability criterion of the endemic equilibrium E * of system (7) in Theorem 5 in Section 4. Furthermore, the CFOC system (27) was investigated by suitably defining the control set U in Equation (28) and the objective function (29) in Section 5. A generalized FOCP was formed in (31) and optimality conditions were proven in detail for this FOCP denoted by (42). Then, the formulas in (42) were applied to achieve the necessary and sufficient optimality conditions for system (27), and also the values of the optimal control pair u 1 * and u 2 * and the co-state or the adjoint equations with the corresponding transversality conditions were described elaborately in Theorem 7.
In 2021, Ghosh et al. [19] described three strategies and, among them, Strategy-III was found to be the most effective one. However, serious matters of concern in determining a realistic and accurate therapy for leprosy are the high cost and extreme adverse effects of the combined drug therapy. Here, in this research article, from Figure 8d and Figure 9d, we can observe that after introducing the memory effect, the amount of the drugs needed initially is much less for system (27) than for Strategy-III in both the memory-free model and integer-ordered system, which is much better in terms of cost-effectiveness and safe therapeutics.
Furthermore, comparing Figure 8 and Figure 9, we can notice that if the value of ζ is decreased from 0.9 to 0.6 , oscillatory solutions appear for the specific range of the parameter set, but in both cases, under the optimal treatment policy, the densities of the cell populations approach a stable concentration. The dynamics of the optimal control profile of Ofloxacin remain almost similar for both ζ = 0.9 , 0.6 , but the control profile of Dapsone provides notable differences. More specifically, for ζ = 0.6 , the drug dose of Dapsone, i.e., u 2 * ( t ) , needs to be increased after 1100 days to tackle the highly Dapsone-drug-resistant M. leprae [39,54] and the associated infection procedure. Thus, to build a perfect regimen for combined therapy, acquiring enough knowledge from the previous memory states about drug-resistance issues [38,57,58] is essential and, hence, CFOC system (27) with ζ = 0.6 is very fruitful in this context. The more we reduce the value of ζ and approach the previous memory states, the more accurate the precision will be. Still, future works on leprosy in this aspect should also focus on investigating the memory stages for the value of ζ < 0.6 , and before implementation, the outcomes should be validated properly by clinical and experimental researchers.
The memory effect for a biological problem actually involves the entire history before the present instant, and the implication of fractional calculus can be viewed as a suitable and exclusive tool for modeling these types of phenomena with hysteresis. Our present study reveals that the whole M. leprae-induced infection in leprosy not only depends on the present state but also on the previous memory stages of the infection process. With the different circumstances and effects that the whole infection mechanism experiences in various memory stages, the CF fractionalized system responds differently to the current conditions. The knowledge accumulated by the metabolic activities in preceding stages confers fitness to the bacteria in an evolutionary way, which validates that the history-dependent drug-resistant behavior is essentially a clear manifestation of the memory effect in leprosy. Thus, to overcome the drug resistance scenario, high-cost-related problems and side-effects of the combined therapy more accurately and realistically, the three-dimensional Caputo–Fabrizio fractionalized optimal control mathematical model (27) for the fractional order ζ = 0.6 presented in our research article should certainly be considered by mathematical and clinical scientists for all future studies on leprosy.

Author Contributions

Conceptualization, S.G.; methodology, S.G. and P.K.R.; software, S.R.; validation, S.G. and S.R.; formal analysis, S.G. and H.B.; investigation, S.G. and X.C.; writing—original draft preparation, S.G.; writing—review and editing, S.G.; visualization, S.G., X.C. and H.B.; supervision, P.K.R. and X.C.; funding acquisition, X.C. and P.K.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research work was supported by the Major Project of National Statistical Science Foundation of China (2021LD01), School of Mathematics and Statistics, Beijing Technology and Business University, Beijing, P.R. China. Furthermore, the research was funded by the W.B. Government State Fellowship, Government of West Bengal, India, and DST-FIST, Department of Mathematics, Jadavpur University, Kolkata 700032, India.

Data Availability Statement

Not applicable.

Acknowledgments

We sincerely thank the anonymous reviewers for their valuable suggestions which strengthened our manuscript greatly in various aspects.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Schreuder, P.A.; Noto, S.; Richardus, J.H. Epidemiologic trends of leprosy for the 21st century. Clin. Dermatol. 2016, 34, 24–31. [Google Scholar] [CrossRef] [PubMed]
  2. World Health Organization. Global leprosy situation, 2010. Wkly. Epidemiol. Rec. 2010, 85, 337–348. [Google Scholar]
  3. World Health Organization. Global leprosy situation, 2012. Wkly. Epidemiol. Rec. 2012, 87, 317–328. [Google Scholar]
  4. Ghosh, S.; Saha, S.; Roy, P.K. Critical observation of WHO recommended multidrug therapy on the disease leprosy through mathematical study. J. Theor. Biol. 2023, 567, 111496. [Google Scholar] [CrossRef] [PubMed]
  5. Richardus, J.H.; Habbema, J.D.F. The impact of leprosy control on the transmission of M. leprae: Is elimination being attained? Lepr. Rev. 2007, 78, 330–337. [Google Scholar] [CrossRef] [PubMed]
  6. Westerlund, S. Dead matter has memory! Phys. Scr. 1991, 43, 174. [Google Scholar] [CrossRef]
  7. Wolf, D.M.; Fontaine-Bodin, L.; Bischofs, I.; Price, G.; Keasling, J.; Arkin, A.P. Memory in microbes: Quantifying history-dependent behavior in a bacterium. PLoS ONE 2008, 3, e1700. [Google Scholar] [CrossRef]
  8. Yang, C.Y.; Bialecka-Fornal, M.; Weatherwax, C.; Larkin, J.W.; Prindle, A.; Liu, J.; Garcia-Ojalvo, J.; Süel, G.M. Encoding membrane-potential-based memory within a microbial community. Cell Syst. 2020, 10, 417–423. [Google Scholar] [CrossRef]
  9. Samko, S.G.; Kilbas, A.A.; Marichev, O.I. Fractional Integrals and Derivatives; Gordon and Breach Science Publishers: Yverdon-les-Bains, Switzerland, 1993; Volume 1. [Google Scholar]
  10. Miller, K.S.; Ross, B. An Introduction to the Fractional Calculus and Fractional Differential Equations; John Wiley & Sons: New York, NY, USA, 1993. [Google Scholar]
  11. Ghosh, M.K.; Sardar, T.; Cao, X.; Roy, P.K. Mathematical study of a memory induced biochemical system. IEEE/CAA J. Autom. Sin. 2018, 5, 1142–1149. [Google Scholar] [CrossRef]
  12. Rossikhin, Y.A.; Shitikova, M.V. Application of fractional calculus for dynamic problems of solid mechanics: Novel trends and recent results. Appl. Mech. Rev. 2010, 63, 010801. [Google Scholar] [CrossRef]
  13. Stiassnie, M. On the application of fractional calculus for the formulation of viscoelastic models. Appl. Math. Model. 1979, 3, 300–302. [Google Scholar] [CrossRef]
  14. Lundstrom, B.N.; Higgs, M.H.; Spain, W.J.; Fairhall, A.L. Fractional differentiation by neocortical pyramidal neurons. Nat. Neurosci. 2008, 11, 1335–1342. [Google Scholar] [CrossRef] [PubMed]
  15. Caputo, M. Linear models of dissipation whose Q is almost frequency independent—II. Geophys. J. Int. 1967, 13, 529–539. [Google Scholar] [CrossRef]
  16. Bagley, R.L.; Torvik, P.J. Fractional calculus-a different approach to the analysis of viscoelastically damped structures. AIAA J. 1983, 21, 741–748. [Google Scholar] [CrossRef]
  17. Mathai, A.M.; Haubold, H.J. Special Functions for Applied Scientists; Springer: New York, NY, USA, 2008; Volume 4. [Google Scholar]
  18. Kneller, G.R.; Hinsen, K. Fractional Brownian dynamics in proteins. J. Chem. Phys. 2004, 121, 10278–10283. [Google Scholar] [CrossRef] [PubMed]
  19. Ghosh, S.; Chatterjee, A.N.; Roy, P.K.; Grigorenko, N.; Khailov, E.; Grigorieva, E. Mathematical Modeling and Control of the Cell Dynamics in Leprosy. Comput. Math. Model. 2021, 32, 52–74. [Google Scholar] [CrossRef]
  20. Du, M.; Wang, Z.; Hu, H. Measuring memory with the order of fractional derivative. Sci. Rep. 2013, 3, 3431. [Google Scholar] [CrossRef]
  21. Chinnathambi, R.; Rihan, F.A.; Alsakaji, H.J. A fractional-order model with time delay for tuberculosis with endogenous reactivation and exogenous reinfections. Math. Methods Appl. Sci. 2021, 44, 8011–8025. [Google Scholar] [CrossRef]
  22. Al-Mdallal, Q.M.; Omer, A.S.A. Fractional-order Legendre-collocation method for solving fractional initial value problems. Appl. Math. Comput. 2018, 321, 74–84. [Google Scholar] [CrossRef]
  23. Al-Mdallal, Q.M.; Hajji, M.A. A convergent algorithm for solving higher-order nonlinear fractional boundary value problems. Fract. Calc. Appl. Anal. 2015, 18, 1423–1440. [Google Scholar] [CrossRef]
  24. Singh, J.; Kumar, D.; Baleanu, D. New aspects of fractional Biswas–Milovic model with Mittag–Leffler law. Math. Model. Nat. Phenom. 2019, 14, 303. [Google Scholar] [CrossRef]
  25. Owolabi, K.M.; Atangana, A. Analysis and application of new fractional Adams–Bashforth scheme with Caputo–Fabrizio derivative. Chaos Solitons Fractals 2017, 105, 111–119. [Google Scholar] [CrossRef]
  26. Kumar, D.; Singh, J.; Al Qurashi, M.; Baleanu, D. Analysis of logistic equation pertaining to a new fractional derivative with non-singular kernel. Adv. Mech. Eng. 2017, 9, 1–8. [Google Scholar] [CrossRef]
  27. Tateishi, A.A.; Ribeiro, H.V.; Lenzi, E.K. The role of fractional time-derivative operators on anomalous diffusion. Front. Phys. 2017, 5, 52. [Google Scholar] [CrossRef]
  28. Atangana, A.; Alkahtani, B.S.T. Analysis of the Keller–Segel model with a fractional derivative without singular kernel. Entropy 2015, 17, 4439–4453. [Google Scholar] [CrossRef]
  29. Yang, L.; Su, Y.; Yang, X.; Wang, Z. The Optimal Control Strategy of Virus Transmission Based on Caputo-Fabrizio Order. Front. Phys. 2021, 9, 731972. [Google Scholar] [CrossRef]
  30. Owolabi, K.M.; Patidar, K.C.; Shikongo, A. Numerical solution for a problem arising in angiogenic signalling. AIMS Math. 2019, 4, 43–63. [Google Scholar] [CrossRef]
  31. Pontryagin, L.S. Mathematical Theory of Optimal Processes; CRC Press: Boca Raton, FL, USA, 1987. [Google Scholar]
  32. Okyere, E.; Oduro, F.T.; Amponsah, S.K.; Dontwi, I.K. Fractional order optimal control model for malaria infection. arXiv 2016, arXiv:1607.01612. [Google Scholar]
  33. Agrawal, O.P. A general formulation and solution scheme for fractional optimal control problems. Nonlinear Dyn. 2004, 38, 323–337. [Google Scholar] [CrossRef]
  34. Agrawal, O.P. A formulation and numerical scheme for fractional optimal control problems. J. Vib. Control 2008, 14, 1291–1299. [Google Scholar] [CrossRef]
  35. Agrawal, O.P.; Baleanu, D. A Hamiltonian formulation and a direct numerical scheme for fractional optimal control problems. J. Vib. Control 2007, 13, 1269–1281. [Google Scholar] [CrossRef]
  36. Baba, B.A.; Bilgehan, B. Optimal control of a fractional order model for the COVID-19 pandemic. Chaos Solitons Fractals 2021, 144, 110678. [Google Scholar] [CrossRef] [PubMed]
  37. Cambau, E.; Saunderson, P.; Matsuoka, M.; Cole, S.T.; Kai, M.; Suffys, P.; Rosa, P.S.; Williams, D.; Gupta, U.D.; Lavania, M.; et al. Antimicrobial resistance in leprosy: Results of the first prospective open survey conducted by a WHO surveillance network for the period 2009–15. Clin. Microbiol. Infect. 2018, 24, 1305–1310. [Google Scholar] [CrossRef] [PubMed]
  38. Benjak, A.; Avanzi, C.; Singh, P.; Loiseau, C.; Girma, S.; Busso, P.; Fontes, A.N.B.; Miyamoto, Y.; Namisato, M.; Bobosha, K.; et al. Phylogenomics and antimicrobial resistance of the leprosy bacillus Mycobacterium leprae. Nat. Commun. 2018, 9, 352. [Google Scholar] [CrossRef]
  39. Williams, D.L.; Araujo, S.; Stryjewska, B.M.; Scollard, D. Dapsone resistance in leprosy patients originally from American Samoa, United States, 2010–2012. Emerg. Infect. Dis. 2018, 24, 1584–1585. [Google Scholar] [CrossRef] [PubMed]
  40. Matsuoka, M.; Budiawan, T.; Aye, K.S.; Kyaw, K.; Tan, E.V.; Cruz, E.D.; Gelber, R.; Saunderson, P.; Balagon, V.; Pannikar, V. The frequency of drug resistance mutations in Mycobacterium leprae isolates in untreated and relapsed leprosy patients from Myanmar, Indonesia and the Philippines. Lepr. Rev. 2007, 78, 343–352. [Google Scholar] [CrossRef]
  41. Aubry, A.; Rosa, P.S.; Chauffour, A.; Fletcher, M.L.; Cambau, E.; Avanzi, C. Drug resistance in leprosy: An update following 70 years of chemotherapy. Infect. Dis. Now 2022, 52, 243–251. [Google Scholar] [CrossRef]
  42. Ahuja, M.; Pathak, V.K. Ofloxacin resistance in multibacillary new leprosy cases from Purulia, West Bengal: A threat to effective secondary line treatment for rifampicin-resistant leprosy cases. J. Glob. Antimicrob. Resist. 2022, 30, 282–285. [Google Scholar] [CrossRef]
  43. Vellappi, M.; Kumar, P.; Govindaraj, V.; Albalawi, W. An optimal control problem for mosaic disease via Caputo fractional derivative. Alex. Eng. J. 2022, 61, 8027–8037. [Google Scholar] [CrossRef]
  44. Caputo, M.; Fabrizio, M. A new definition of fractional derivative without singular kernel. Prog. Fract. Differ. Appl. 2015, 1, 73–85. [Google Scholar]
  45. Losada, J.; Nieto, J.J. Properties of a new fractional derivative without singular kernel. Prog. Fract. Differ. Appl. 2015, 1, 87–92. [Google Scholar]
  46. Qing, Y.; Rhoades, B.E. T-stability of Picard iteration in metric spaces. Fixed Point Theory Appl. 2008, 2008, 418971. [Google Scholar] [CrossRef]
  47. Gao, F.; Li, X.; Li, W.; Zhou, X. Stability analysis of a fractional-order novel hepatitis B virus model with immune delay based on Caputo-Fabrizio derivative. Chaos Solitons Fractals 2021, 142, 110436. [Google Scholar] [CrossRef]
  48. Li, H.; Cheng, J.; Li, H.B.; Zhong, S.M. Stability analysis of a fractional-order linear system described by the Caputo–Fabrizio derivative. Mathematics 2019, 7, 200. [Google Scholar] [CrossRef]
  49. Kamocki, R. Pontryagin maximum principle for fractional ordinary optimal control problems. Math. Methods Appl. Sci. 2014, 37, 1668–1686. [Google Scholar] [CrossRef]
  50. Chatterjee, A.N.; Al Basir, F.; Almuqrin, M.A.; Mondal, J.; Khan, I. SARS-CoV-2 infection with lytic and non-lytic immune responses: A fractional order optimal control theoretical study. Results Phys. 2021, 26, 104260. [Google Scholar] [CrossRef]
  51. Ghosh, S.; Rana, S.; Roy, P.K. Leprosy: Considering the Effects on Density-Dependent Growth of Mycobacterium leprae. Differ. Equations Dyn. Syst. 2022, 1–15. [Google Scholar] [CrossRef]
  52. Masaki, T.; Qu, J.; Cholewa-Waclaw, J.; Burr, K.; Raaum, R.; Rambukkana, A. Reprogramming adult Schwann cells to stem cell-like cells by leprosy bacilli promotes dissemination of infection. Cell 2013, 152, 51–67. [Google Scholar] [CrossRef]
  53. Khan, M.A.; Atangana, A. Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative. Alex. Eng. J. 2020, 59, 2379–2389. [Google Scholar] [CrossRef]
  54. Wu, Z.; Wang, C.; Wang, Z.; Shi, Y.; Jiang, H.; Wang, H. Risk factors for Dapsone Resistance in Leprosy Patients: A systematic meta-analysis. J. Glob. Antimicrob. Resist. 2022, 30, 459–467. [Google Scholar] [CrossRef]
  55. Maymone, M.B.; Venkatesh, S.; Laughter, M.; Abdat, R.; Hugh, J.; Dacso, M.M.; Rao, P.N.; Stryjewska, B.M.; Dunnick, C.A.; Dellavalle, R.P. Leprosy: Treatment and management of complications. J. Am. Acad. Dermatol. 2020, 83, 17–30. [Google Scholar] [CrossRef] [PubMed]
  56. Rosa, P.S.; D’Espindula, H.R.; Melo, A.C.; Fontes, A.N.; Finardi, A.J.; Belone, A.F.; Sartori, B.G.; Pires, C.A.; Soares, C.T.; Marques, F.B.; et al. Emergence and transmission of drug-/multidrug-resistant Mycobacterium leprae in a former leprosy colony in the brazilian amazon. Clin. Infect. Dis. 2020, 70, 2054–2061. [Google Scholar] [CrossRef] [PubMed]
  57. Matsuoka, M. Drug resistance in leprosy. Jpn. J. Infect. Dis. 2010, 63, 1–7. [Google Scholar] [CrossRef] [PubMed]
  58. Williams, D.L.; Gillis, T.P. Drug-resistant leprosy: Monitoring and current status. Lepr. Rev. 2012, 83, 269–281. [Google Scholar] [CrossRef]
Figure 1. Behavior of the solutions of the system populations and 3-D phase portrait diagrams of the CF fractionalized system (7) depicting the oscillatory dynamics and appearance of limit cycles for ζ = 1 . Values of ν 2 = 0.03 , S u m a x = 1000 were used to simulate the subfigures in this figure, where all the other parameter values were chosen from Table 1. (a) Behavior of the trajectories of system (7); (b) 3-D phase diagram for system (7) in S u S i B l space.
Figure 1. Behavior of the solutions of the system populations and 3-D phase portrait diagrams of the CF fractionalized system (7) depicting the oscillatory dynamics and appearance of limit cycles for ζ = 1 . Values of ν 2 = 0.03 , S u m a x = 1000 were used to simulate the subfigures in this figure, where all the other parameter values were chosen from Table 1. (a) Behavior of the trajectories of system (7); (b) 3-D phase diagram for system (7) in S u S i B l space.
Mathematics 11 03630 g001
Figure 2. Behavior of the solutions of the system populations and 3-D phase portrait diagrams of the CF fractionalized system (7) depicting the stable behavior for ζ = 1 . Values of ν 2 = 0.05 and S u m a x = 1000 were used to simulate the subfigures in this figure where other parameter values were chosen from Table 1. (a) Behavior of the trajectories of system (7); (b) 3-D phase diagram for system (7) in S u S i B l space.
Figure 2. Behavior of the solutions of the system populations and 3-D phase portrait diagrams of the CF fractionalized system (7) depicting the stable behavior for ζ = 1 . Values of ν 2 = 0.05 and S u m a x = 1000 were used to simulate the subfigures in this figure where other parameter values were chosen from Table 1. (a) Behavior of the trajectories of system (7); (b) 3-D phase diagram for system (7) in S u S i B l space.
Mathematics 11 03630 g002
Figure 3. Time series and phase portrait diagram of the CF fractionalized system (7) depicting the unstable oscillatory behavior of the system state populations and appearance of stable limit cycles for ζ = 0.8 . Values of ν 2 = 0.03 and S u m a x = 1200 were used to simulate the subfigures in this figure where other parameter values were chosen from Table 1. (a) Behavior of the trajectories of system (7); (b) 3-D phase diagram for system (7) in S u S i B l space.
Figure 3. Time series and phase portrait diagram of the CF fractionalized system (7) depicting the unstable oscillatory behavior of the system state populations and appearance of stable limit cycles for ζ = 0.8 . Values of ν 2 = 0.03 and S u m a x = 1200 were used to simulate the subfigures in this figure where other parameter values were chosen from Table 1. (a) Behavior of the trajectories of system (7); (b) 3-D phase diagram for system (7) in S u S i B l space.
Mathematics 11 03630 g003
Figure 4. Time series and phase portrait diagram of the CF fractionalized system (7) depicting the asymptotically stable behavior of the system state populations for ζ = 0.8 . Here, values of ν 2 = 0.05 and S u m a x = 1200 were used to simulate the subfigures in this figure. Values of all the other parameters were chosen from Table 1. (a) Behavior of the trajectories of system (7); (b) 3-D phase diagram for system (7) in S u S i B l space.
Figure 4. Time series and phase portrait diagram of the CF fractionalized system (7) depicting the asymptotically stable behavior of the system state populations for ζ = 0.8 . Here, values of ν 2 = 0.05 and S u m a x = 1200 were used to simulate the subfigures in this figure. Values of all the other parameters were chosen from Table 1. (a) Behavior of the trajectories of system (7); (b) 3-D phase diagram for system (7) in S u S i B l space.
Mathematics 11 03630 g004
Figure 5. Times series and phase portrait diagrams of the CF fractionalized system (7) depicting the sustained oscillatory unstable behavior of the system state populations and appearance of limit cycles for ζ = 1 . Here, values of ν 2 = 0.035 and S u m a x = 1200 were used for simulating the subfigures in this figure. All the other parameter values were chosen from Table 1. (a) Behavior of the trajectories of system (7); (b) 3-D phase diagram for system (7) in S u S i B l space; (c) 2-D phase diagram of system (7) in the S u B l plane; (d) 2-D phase diagram for system (7) in the S u S i plane; (e) 2-D phase diagram for system (7) in the S i B l plane.
Figure 5. Times series and phase portrait diagrams of the CF fractionalized system (7) depicting the sustained oscillatory unstable behavior of the system state populations and appearance of limit cycles for ζ = 1 . Here, values of ν 2 = 0.035 and S u m a x = 1200 were used for simulating the subfigures in this figure. All the other parameter values were chosen from Table 1. (a) Behavior of the trajectories of system (7); (b) 3-D phase diagram for system (7) in S u S i B l space; (c) 2-D phase diagram of system (7) in the S u B l plane; (d) 2-D phase diagram for system (7) in the S u S i plane; (e) 2-D phase diagram for system (7) in the S i B l plane.
Mathematics 11 03630 g005
Figure 6. Times series and phase portrait diagrams of the CF fractionalized system (7) depicting the asymptotically stable dynamics of the system state populations for ζ = 0.6 . Here, major observations are noted by studying the comparative behavior of system (7) between the CF-fractional order ζ = 1 (unstable in the form of stable periodic solutions and limit cycles shown in Figure 5) and ζ = 0.6 (stable behavior demonstrated in this figure, i.e., in Figure 6). Here, values of ν 2 = 0.035 and S u m a x = 1200 were used for simulating the subfigures in this figure. All the other parameter values were chosen from Table 1. (a) Behavior of the trajectories of system (7); (b) 3-D phase diagram for system (7) in S u S i B l space; (c) 2-D phase diagram of system (7) in the S u B l plane; (d) 2-D phase diagram for system (7) in the S u S i plane; (e) 2-D phase diagram for system (7) in the S i B l plane.
Figure 6. Times series and phase portrait diagrams of the CF fractionalized system (7) depicting the asymptotically stable dynamics of the system state populations for ζ = 0.6 . Here, major observations are noted by studying the comparative behavior of system (7) between the CF-fractional order ζ = 1 (unstable in the form of stable periodic solutions and limit cycles shown in Figure 5) and ζ = 0.6 (stable behavior demonstrated in this figure, i.e., in Figure 6). Here, values of ν 2 = 0.035 and S u m a x = 1200 were used for simulating the subfigures in this figure. All the other parameter values were chosen from Table 1. (a) Behavior of the trajectories of system (7); (b) 3-D phase diagram for system (7) in S u S i B l space; (c) 2-D phase diagram of system (7) in the S u B l plane; (d) 2-D phase diagram for system (7) in the S u S i plane; (e) 2-D phase diagram for system (7) in the S i B l plane.
Mathematics 11 03630 g006
Figure 7. Representation of the stability regions for the Caputo–Fabrizio (CF) fractionalized system (7) for different values of CF fractional order ζ . Values of the parameters were chosen as ν 1 = 0.4 , B l m a x = 530 , μ = 0.1 , β 1 = 0.00032 , σ = 0.3 , β 2 = 0.00024 and all the values of other parameters were taken from Table 1. (a) Stability region for CF fractional order ζ = 1 ; (b) stability region for CF fractional order ζ = 0.8 ; (c) stability region for CF fractional order ζ = 0.6 .
Figure 7. Representation of the stability regions for the Caputo–Fabrizio (CF) fractionalized system (7) for different values of CF fractional order ζ . Values of the parameters were chosen as ν 1 = 0.4 , B l m a x = 530 , μ = 0.1 , β 1 = 0.00032 , σ = 0.3 , β 2 = 0.00024 and all the values of other parameters were taken from Table 1. (a) Stability region for CF fractional order ζ = 1 ; (b) stability region for CF fractional order ζ = 0.8 ; (c) stability region for CF fractional order ζ = 0.6 .
Mathematics 11 03630 g007
Figure 8. Trajectories of the cell populations of the optimal-control-induced Caputo–Fabrizio fractional (CFOC) system (27). Scenarios for both with and without control are exhibited here in subfigures (a), (b) and (c), respectively, denoted by green and red color. In subfigure (d), optimal control profiles of u 1 * and u 2 * are shown for CFOC system (27). Values of the parameters were chosen as ζ = 0.9 , ν 2 = 0.03 , S u m a x = 1200 for simulating the subfigures of this figure and other parameter values were selected from Table 1. (a) Behavior of the densities of healthy Schwann cells S u ; (b) behavior of the densities of infected Schwann cells S i ; (c) behavior of the concentrations of M. leprae bacteria B l ; (d) dynamics of the optimal control profiles u 1 * and u 2 * for system (27).
Figure 8. Trajectories of the cell populations of the optimal-control-induced Caputo–Fabrizio fractional (CFOC) system (27). Scenarios for both with and without control are exhibited here in subfigures (a), (b) and (c), respectively, denoted by green and red color. In subfigure (d), optimal control profiles of u 1 * and u 2 * are shown for CFOC system (27). Values of the parameters were chosen as ζ = 0.9 , ν 2 = 0.03 , S u m a x = 1200 for simulating the subfigures of this figure and other parameter values were selected from Table 1. (a) Behavior of the densities of healthy Schwann cells S u ; (b) behavior of the densities of infected Schwann cells S i ; (c) behavior of the concentrations of M. leprae bacteria B l ; (d) dynamics of the optimal control profiles u 1 * and u 2 * for system (27).
Mathematics 11 03630 g008
Figure 9. Trajectories of the cell populations of the optimal-control-induced Caputo–Fabrizio fractional (CFOC) system (27). Scenarios for both with and without control are exhibited here in subfigures (a), (b) and (c), respectively, denoted by green and red colors. In subfigure (d), optimal control profiles of u 1 * and u 2 * are shown for CFOC system (27). Values of the parameters were chosen as ζ = 0.6 , ν 2 = 0.03 , S u m a x = 1200 for simulating the subfigures of this figure and other parameter values were selected from Table 1. (a) Behavior of the densities of healthy Schwann cells S u ; (b) behavior of the densities of infected Schwann cells S i ; (c) behavior of the concentrations of M. leprae bacteria B l ; (d) dynamics of the optimal control profiles u 1 * and u 2 * for system (27).
Figure 9. Trajectories of the cell populations of the optimal-control-induced Caputo–Fabrizio fractional (CFOC) system (27). Scenarios for both with and without control are exhibited here in subfigures (a), (b) and (c), respectively, denoted by green and red colors. In subfigure (d), optimal control profiles of u 1 * and u 2 * are shown for CFOC system (27). Values of the parameters were chosen as ζ = 0.6 , ν 2 = 0.03 , S u m a x = 1200 for simulating the subfigures of this figure and other parameter values were selected from Table 1. (a) Behavior of the densities of healthy Schwann cells S u ; (b) behavior of the densities of infected Schwann cells S i ; (c) behavior of the concentrations of M. leprae bacteria B l ; (d) dynamics of the optimal control profiles u 1 * and u 2 * for system (27).
Mathematics 11 03630 g009aMathematics 11 03630 g009b
Table 1. List of parameter values used in numerical simulation for systems (7) and (27) for ζ = 1 .
Table 1. List of parameter values used in numerical simulation for systems (7) and (27) for ζ = 1 .
ParameterParameter DefinitionAssigned ValueReferences
ν 1 growth rate of S u 0.4 [19,51]
ν 2 growth rate of B l 0.01–0.05 [19,52]
S u m a x carrying capacity of S u 600–1400 [19,51,52]
B l m a x carrying capacity of B l 400–550 [19,51,52]
μ natural death rate of S i 0.1 Estimated
β 1 infection rate0.0003–0.0046 [19]
σ proliferation rate of B l 0.1–0.35 [19,51]
β 2 clearance rate of B l 0.00015–0.0003 [19,51]
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Cao, X.; Ghosh, S.; Rana, S.; Bose, H.; Roy, P.K. Application of an Optimal Control Therapeutic Approach for the Memory-Regulated Infection Mechanism of Leprosy through Caputo–Fabrizio Fractional Derivative. Mathematics 2023, 11, 3630. https://doi.org/10.3390/math11173630

AMA Style

Cao X, Ghosh S, Rana S, Bose H, Roy PK. Application of an Optimal Control Therapeutic Approach for the Memory-Regulated Infection Mechanism of Leprosy through Caputo–Fabrizio Fractional Derivative. Mathematics. 2023; 11(17):3630. https://doi.org/10.3390/math11173630

Chicago/Turabian Style

Cao, Xianbing, Salil Ghosh, Sourav Rana, Homagnic Bose, and Priti Kumar Roy. 2023. "Application of an Optimal Control Therapeutic Approach for the Memory-Regulated Infection Mechanism of Leprosy through Caputo–Fabrizio Fractional Derivative" Mathematics 11, no. 17: 3630. https://doi.org/10.3390/math11173630

APA Style

Cao, X., Ghosh, S., Rana, S., Bose, H., & Roy, P. K. (2023). Application of an Optimal Control Therapeutic Approach for the Memory-Regulated Infection Mechanism of Leprosy through Caputo–Fabrizio Fractional Derivative. Mathematics, 11(17), 3630. https://doi.org/10.3390/math11173630

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop