Mathematical Analysis of Biodegradation Model under Nonlocal Operator in Caputo Sense

: To lower the concentration of organic pollutants in the efﬂuent stream, wastewater must be treated before being discharged into the environment. The question of whether wastewater treatment facilities can successfully reduce the concentration of micropollutants found in their inﬂuent streams is becoming increasingly pressing. The removal of micropollutants in treatment plants is investigated using a model that incorporates biodegradation and sorption as the key processes of micropollutant removal. This article provides the mathematical analysis of the wastewater model that describes the removal of micropollutant in treatment plants under a non-local operator in Caputo sense. The positivity of the solution is presented for the Caputo fractional model. The steady state’s solution of model and their stability is presented. The ﬁxed point theorems of Leray–Schauder and Banach are used to deduce results regarding the existence of the solution of the model. Ulam–Hyers (UH) types of stabilities are presented via functional analysis. The fractional Euler method is used to ﬁnd the numerical results of the proposed model. The numerical results are illustrated via graphs to show the effects of recycle ratio and the impact of fractional order on the evolution of the model.


Introduction and Motivation
Traditional wastewater treatment plants (WTP) are capable of removing bulk carboncontaining organic pollutants and also nutrients like nitrogen and phosphorus. They were not, however, constructed to eliminate organic pollutants identified in trace quantities, sometimes known as micropollutants [1]. Chemicals found in personal care items, such as active pharmaceutical ingredients and derivatives, as well as chemicals found in household goods, such as insecticides and surfactants, are examples of common micropollutants. Polycyclic aromatic hydrocarbons (PAHs), hormones, and metals, are examples of frequent micropollutants. Micropollutants have been proven to have an ecotoxic effect when released into aquatic environments. As a result, European legislation has been enacted that requires both businesses and governments to cut their emissions [2]. Pomies et al. [3] evaluated mathematical models for the removal of micropollutants in WTPs. In WTPs, micropollutants are eliminated via four major methods [3]. Biodegradation, cometabolism, mass transfer to the gas phase, and sorption onto particles are all examples of these processes. Not all pollutants are affected by these processes. By mass transfer into the gas phase, for instance, only volatile organic molecules are eliminated. For some pollutants, other processes may apply. Precipitation can, for example, eliminate heavy metals. One of the uses of mathematical models is determining the most important removal method for a specific pollutant. The activity of biomass in WTP can help eliminate micropollutants. This is also known as micropollutant bioconversion or biotransformation. Biodegradation and cometabolism are two mechanisms that might cause this change. In the context of the removal of micropollutants, there is no universally accepted definition of biodegradation [3]. The elimination of micropollutants via a mechanism unrelated to the microorganism's development is a popular view. In addition, it is commonly thought that biodegradation occurs solely in the aqueous phase. Delgadillo-Mirquez et al. [4] have looked at this hypothesis. The term "cometabolism" refers to the breakdown of micropollutants in the presence of another quickly degradable substrate. Micropollutants are not a carbon source for microbial development in this scenario [5]. A metabolic process degrades a large number of organic micropollutants found in wastewater treatment plants. Many additional models [6][7][8][9] have been developed to explore WTPs in different ways. Nelson et al. presented a mathematical model in [10] that shows alternative techniques to removing micropollutants in the activated sludge process. They outline the biological events and processes that their model includes. Three biological reactions can be found. The activated sludge mechanism is linked to two of these, Equations (1) and (2). The degradation of micropollutants by biomass is the third process, shown in Equation (3). There are two physical processes in the model: The sorption of micropollutants onto suspended particles to generate particulate micropollutants as shown in Equation (4), and the desorption of particulates as shown in Equation (5). The biomass (X B,H ) consumes the soluble substrate S s to create new biomass.
where Y H denotes the coefficient of heterotrophic yield. The following reaction represents the death of particulate biomass: where b H denotes the coefficient of heterotrophic decay. Biological removal of micropollutants is given by: Sorption is given by: where k sor denotes the kinetic constant for adsorption and TSS represents the total suspended solids. Desorption is given by: where k des denotes the kinetic constant for desorption. The Nelson model contains four contents explaining the activated sludge process and micropollutants. They formulated the dimensional model as: where M 2 represents the monod growth kinetics. M 2 and residence time τ are, respectively, defined as: The parameters k d , R bio and X TSS are, respectively, defined as: where c 2 represents the conversion of chemical demand to total suspended solids. Model (6) may now be converted to a dimensionless version by using the dimensionless variables as follows: The dimensionless model is described as follows: where The initial values are S * (0) = S * 0 > 0, X * (0) = X * 0 ≥ 0, C * s (0) = C * s 0 ≥ 0, and C * p (0) = C * p 0 ≥ 0. The description of all variables and parameters in model (6) and (7) are given below.
• S s , X B,H , C s and C p denote the soluble substrate, heterotroph (particulate biomass), soluble micropollutants, and particulate micropollutants, respectively. • S s,in , F, and V represent the substrate concentration in the feed, the flow rate through bioreactor, and the volume of the bioreactor, respectively. • C, X B,H,in and R denote the recycle concentration factor, heterotroph concentration, and the recycle ratio, respectively. • b H , C s,in , and X TSS denote the decay coefficient of particulate biomass, soluble micropollutants concentration, and the total suspended solids, respectively. • R bio , τ represent the biological removal rate and residence time.
• H, Q air , and C p,in represent the Henry coefficient of micropollutant, the flow rate of aeration, and the concentration of particulate micropollutants in the feed.
The memory and past history of a process can be determine by nonlocal operators in fractional calculus (FC). FC has lately been utilized to explain a wide range of difficult issues in the applied sciences [11,12]. Researchers have created mathematical theories to replicate the intricacies of nature using FC and to study memory and heredity features of a physical process. To construct mathematical models of real-world problems, differential operators (DOs) are used. DOs can be local or nonlocal (fractional). The three sorts of fractional operators that are most familiar are: The Caputo operator has a kernel with a power-law distribution [13]; a DO with exponential-decay law-called Caputo-Fabrizio operator [14]; and finally, an operator with Mittag-Leffler law-called the Atangana-Baleanu operator [15]. Researchers frequently use these operators to incorporate memory characteristics of theoretical models, which may be encountered in a wide range of disciplines of study [16][17][18].
Researchers have investigated the fractional differential equations from different points of view. Boundary value problems for Hilfer fractional differential inclusions with nonlocal integral boundary conditions have been studied by Wongcharoen et al. in [19]. Existence and Ulam-Hyers-Mittag-Leffler stability results of Ψ Hilfer nonlocal Cauchy problem is given in [20]. Hyers-Ulam stability and existence of solutions to the generalized Liouville-Caputo fractional differential equations have been investigated in [21]. Ahmed et al. investigated the solutions for the impulsive fractional pantograph differential equation via generalized anti-periodic boundary conditions in [22]. The numerical solution of fractional-order ordinary differential equations using the reformulated infinite state representation was proposed by Hinze and his coauthors [23]. Beside this, FC also has applications in biochemical models. Recently, Alqahtani et al. investigated the fractional order bioethanol production model [24]. Qureshi et al. used the fractional operator to study the blood ethanol model in [25]. In the literature [26], the authors analyzed the WWTPS model via fractional order. Some applications of fractional operators in a biochemical model are included in [27,28]. Inspired from the above studies, we study model (7) under a nonlocal operator in Caputo sense. We write model (7) in a Caputo operator as: subject to the initial conditions S * (0) = S * The fractional operator was not used to study model (7) in the literature. As a result, we examine the biodegradation model using the Caputo derivative. The aforementioned model's theoretical and computational features are investigated. Using well-known fixed point theorems like Banach's and Larey-nonlinear Schauder's alternative, the main purpose of this research is to investigate the existence and uniqueness of the model's solutions. The stability study is also looked at from several types of Ulam stabilities. Finally, we use the fractional Euler technique to get the approximate solutions to the biodegradation model's fractional model. The numerical results are graphically presented to study the evolution of the proposed model against different fractional orders. Now we provide some basic notions of fractional calculus, which are needed in the paper. FD and FI represent the fractional derivative and fractional integral, respectively. Definition 1 ([13]). For a continuous function H on [0, T], the FD in the case of Caputo is defined as: , an FI in Riemann-Liouvilli sense of order 0 < δ ≤ 1 is defined as: A and an open subset P of A , also 0 ∈ P. If U :P →B * r is a continuous and compact function. Then either U possesses a fixed point in P or ∃Y ∈ ∂P and ψ ∈ (0, 1) such that Y = ψU(Y ).

Mathematical Analysis
In this section, we look at some of the model's key features, such as non-negative solutions, the existence of a unique solution for the suggested model, and several forms of Ulam's stability results. Using fractional Euler's technique, we establish the numerical solutions for the investigated model.

Theoretical Analysis
In this section, we examine the theoretical aspects of the proposed model. Lemma 1 (Non-negativity of the solution). If the initial conditions are non-negative then the solutions of (8) are non-negative, for all t > 0.
Proof. As all of the parameters and initial conditions are positive. So from the proposed model we get As a consequence, the suggested model's solutions are all non-negative.

Steady States and Its Stability
• The proposed model has two steady states. The first one is the washout branch (W 0 ), when X * = 0. The washout branch is given by The second one is the no-washout branch (W * ), when X * > 0. The no-washout branch is given by: Now, we present the local stability (LS) of the steady states. Consider the Jacobian matrix as: • For the LS of the washout solution branch, consider the eigenvalues of J at the washout branch as: , the WSS is LS if τ * is enough low: .
• For an LS of no-washout solution branch, from the second equation of model (8), we have For the no-washout branch X * = 0, so For the above equation, the characteristic polynomail of J is as follows: where a = 1 The no-washout solution is of interest when all solution are positive. So, in these cases, the coefficients of J are positive. From the above equations, we see that a > 0, a 3 > 0, Thus b > 0, so no-washout solution branch is LS.

Applications of Fixed Point Theory to the Existence of Solution
We establish that there is at least one unique solution to the given model using the Leray-Schauder theorem and the Banach fixed point theorems in this section. Let us write the model's right-hand side (8) as follows: Model (8) becomes: Application of FI gives: Let Equation (12) gets the form: Let A = C [0, T], R 4 denotes a Banach space associated with a norm: where |Y (t * )| = S * (t * ) + X * (t * ) + C * s (t * ) + C * p (t * ) , where S * , X * , C * s , C * p ∈ A . If (13) admits a fixed point then it must have a solution. Let us define an operator U : A → A . In the view of (13), we define U as: Converting the examined model to a fixed point problem to establish fixed point theory, i.e, Y = UY. The following theorem guarantees that the model has at least one solution.
Theorem 5 (Existence of solution). Assume that: (C 1 ) For all t * ∈ [0, T] and Y ∈ A ∃ α Z > 0 and β Z > 0 such that: and Z is a continuous function. Then, Equation (13) has at least one solution.
Proof. Consider the operator which is defined by (15). First, we prove that U maps bounded balls into bounded balls in A . Consider a bounded ball as: Using mathematics and assumption (C 1 ), we get Thus, U(B * r ) ⊂ B * r . Now, it is shown that U maps bounded balls into equi-continuous sets of A . For Y ∈ B * r and z 1 , Clearly, |(UY )(z 2 ) − (UY )(z 1 )| → 0 as z 2 → z 1 . Thanks to the Arzela-Ascoli theorem, U : A → A is completely continuous. Finally, we show that all solutions to Y = ψUY, where ψ ∈ (0, 1), are bounded. Then, for t * ∈ [0, T], one gets After some manipulation and taking the maximum value of t * ∈ [0, T], one has From assumption (C 2 ), ∃ζ > 0 such that Y = ζ. Now, let P = {Y ∈ A : Y < ζ}. Note that U : P → B * r is completely continuous. From P, Y ∈ ∂P with Y = ψUY, for some ψ ∈ (0, 1). Hence, by the fixed point theorem of Leray-Schauder, U possesses a fixed point Y ∈ P. As an outcome, the considered model must have at least one solution.
Proof. We use the contraction theorem of Banach, to prove that the model under discussion has a unique solution. Let B * r = {Y ∈ A : Y ≤ r}. As proved above, UB * r ⊂ B * r . Lastly, we show that U : B * r → B * r is a contraction. For any Y 1 , Y 2 ∈ A , one gets Since¯h Z T ζ Γ(ζ+1) < 1, it means that U is a contraction. By the Banach contraction theorem, U possesses a unique fixed point. As a result, there is at most one solution to the problem under consideration.

Ulam-Hyers Stability
We will explore the theory for several forms of Ulam-Hyers (UH) stability in this part. We begin by defining several Ulam stabilizers for the model under discussion. Let σ > 0 and a continuous mapping L Z : [0, T] × R 4 → R + . Consider where t * ∈ [0, T] and σ = max(σ j ) T ,for j = 1, 2, 3, 4.

Definition 7.
The model under consideration is UH stable, if for every σ > 0, and for every solution M ∈ A of (12) ∃J Z > 0 and a solution Y ∈ A with where t * ∈ [0, T] and J Z = max(J Z j ) T for j = 1, 2, 3, 4.

Definition 9.
The suggested model is UHR stable if ∃ U L Z > 0 such that for each σ > 0 and for every M ∈ A ∃ a solution Y ∈ A such that where t * ∈ [0, T].

Definition 10.
The suggested model is generalized UHR stable if ∃ U L Z > 0 such that for each σ > 0 and for every M ∈ A ∃ a solution Y ∈ A such that where t * ∈ [0, T].

Remark 1.
Let M ∈ A be a solution of (20) iff ∃G ∈ A possesses the properties as follows:

Remark 2.
Let M ∈ A be a solution of (21) iff ∃B ∈ A with the following properties: . First, we prove fundamental results that are required for discussing Ulam's stabilities in the provided model. We will make another assumption that will be useful in the discussion. We assume that: (A 4 ) For any t * ∈ [0, T], ∃ and increasing function L Z ∈ A and Υ L Z > 0 such that Lemma 2. Let δ ∈ (0, 1] and M ∈ A be a solution of (20); then . (28) Proof. Since M ∈ A is a solution of (20). So by the second point of Remark 1, one gets Using fractional integral, the solution of (29) is given by Using the first part of Remark 1, we have That is the end of the proof.

Lemma 3.
Let δ ∈ (0, 1] and M ∈ A be a solution of (21); then . (30) Proof. Since M ∈ A is a solution of (20). In the context of Remark 2's second part, one may write: In the context of Remark 2's first part, one may write: The proof is completed.
We are now in a position to demonstrate the model's UH and UHR stability.
Theorem 11 (Generalized UH stability). Assume that for all Y ∈ A , the function Z is continuous. Let the assumption (A 1 ) and¯h Z T δ Γ(δ+1) < 1 be satisfied. Then model (12) is UH stable and, consequently, generalized UH stable.
Proof. Assume that σ > 0 and M ∈ A is any solution of the relation (20). Let Y ∈ A be a unique solution of the system (12). Using Equation (13) and Lemma 1, we get After some calculation, one may get . As a result, the requirement of UH stability is achieved. Therefore, model (12) is UH stable. Next, set L Z (σ) = σJ Z with L Z (0) = 0. Hence, model (12) is generalized UH stable.
The UHR and generalized UHR stability of the considered system are shown in the following theorems.
Proof. Assume that σ > 0 and M ∈ A is any solution of the relation (22). Let Y ∈ A be a unique solution of the system (12). Using Equation (13) and Lemma 2, we get Doing simple calculation, one may get . By letting: We reach our goal: As a result, the suggested model is UHR stable. Next, set σ = 1 in the relation which is given above, along with L Z (0) = 0, the considered model is generalized UHR stable.

Computational Analysis
In this section, for the considered model, the numerical results are deduced via Euler's method. Next, we find the numerical result for the first equation of the system (8). Consider the first equation as: To find the solutions of the Equation (32), let [0, d] be the set of points where the solution lies. The solution of equation (32), is E p (t), which cannot be evaluated directly. Instead, we produce a collection of points t * k , t *

k+1
, which we utilize for our iterative method. For this, we divide [0, d] into j sub-intervals t * k , t * k+1 of each width h = d j , utilizing the nodes t * k = kh, where k = 0, 1, . . . , j. Assume that S * (t * ), D δ t * S * (t * ) and D 2δ t * S * (t * ) are continuous on [0, d]. S * (t * ) can be expanded by generalized Taylor's formula about t * = t * 0 = 0, ∀ t * , ∃ C 1 so that: When D δ t * [S * (t * )](t * 0 ) = G 1 (t * 0 , S * (t * 0 )) and h = t * 1 are put into (33), the above equation becomes: Taking h too small results neglecting second order (h 2δ ), Similarly, a collection of points is achieved for the approximate solution, where t * k+1 = t * k + h; therefore, Thus, Equation (35) is a numerical result for the first equation of the considered model. Similarly, for the last two compartments of the considered system, we have Hence, Equations (36)-(38) represent the solution of the proposed model.

Numerical Simulations
In this part, we illustrate the numerical results graphically for various scenario. The parameter values used for the graphical analysis are as follows: µ m ax = 6, S s,in = 0.2, Y H = 0.67, b H = 0.22, K S = 0.02, C s,in = 100 and c 2 = 0.9. For the other parameters, we suggest two types of compound used for micropollutant biotransformation parameters. One is highly biodegradable and the other is slowly biodegradable compounds. The values of k * sor , k * bio , and k * d are given in Table 1. The graphical illustration is discussed for the two different compounds; one is IBP and other is SMX. In Figure 1, the effect of the recycle ratio is presented for IBP and δ = 0.8. After an abrupt increase, it is noticed that the soluble substrate is decreasing when the recycle ratio is increasing. However, particulate biomass is increasing when the recycle ratio increases. The particulate biomass attains its maximum value at the higher value of R * . It is also observed that, as the recycle ratio increases, then soluble micropollutants decrease and particulate micropollutants increase. Thus, the recycle ratio R * has an important role in the activated sludge processes. In Figure 1, the effect of the recycle ratio is presented for SMX and δ = 0.8. After an abrupt increase, it is noticed that the soluble substrate is decreased when the recycle ratio is increasing. However, particulate biomass is increasing when the recycle ratio increases. The particulate biomass attains its maximum value at the higher value of R * = 0.98. It is seen that, as the R * soluble micropollutant increases rapidly and then become stable with the passage of time, particulate micropollutants decrease. We observe different dynamics of the soluble micropullaunt and particulate micropollutants in Figures 1 and 2 because we used two different compounds for biodegradation. The IBP is highly biodegradable while SMX has low biodegradability. In Figures 3 and 4, we show the effect of fractional order on the evolution of the model for IBP and SMX, respectively. As δ increases, the rate of decay or growth becomes slower, and vice versa. Thus, the fractional order δ has an important role in the rate of decay or growth process of the model.

Conclusions
In this paper, we developed a fractional-order mathematical model for the activated sludge process that incorporated the three major processes for micropollutant removal: Biological removal (biodegradation and/or cometabolism), volatilization, and sorption. We used a basic biotransformation model with no cometabolism and a rate of biodegradation that was a linear function of micropollutant concentration. We showed that the fractional order model has a non-negative solution. We calculated two steady states of solution and derived the local stability of both steady state solutions. We deduced the results that are concerned with existence of at least one solution and at most one solution of the fractional order biodegradation model through fixed point theorems of Leray-Schauder and Stephen Banach. We constructed different forms of UH stability via functional analysis. We deduced numerical results of the considered model through the Euler numerical method. We illustrated the achieved results through graphs to study the dynamics of the different classes of the considered model. The elimination of soluble micropollutants is improved when they are highly biodegradable and have a high sorption capacity, according to the results. By either raising the value of the biological degradation kinetic constant or concentrating the substrate in the feed, the limiting concentration can be reduced. Concentrating the feed improves the elimination of biodegradable pollutants. The numerical simulations shown that the recycle ratio has a great impact on the soluble and particulate micropollutants. The effect of fractional order on the evolution of the chemical process of the model has been observed. In future, one could study the considered model in different approaches: • Most of the dynamical systems have been studied for bifurcation theory [29][30][31]. One can study the static bifurcation for the considered model.

Conflicts of Interest:
The authors declare no conflict of interest.