mathematics

Article
Mathematical Analysis of Oxygen Uptake Rate in Continuous
Process under Caputo Derivative
Rubayyi T. Alqahtani 1 , Abdullahi Yusuf 2,3, * and Ravi P. Agarwal 4

1 Department of Mathematics and Statistics, College of Science Al Imam Mohammad Ibn Saud Islamic
University (IMSIU), Riyadh 11566, Saudi Arabia; rr-gahtani@hotmail.com
2 Department of Computer Engineering, Biruni University, 34010 Istanbul, Turkey
3 Department of Mathematics, Federal University Dutse, Jigawa 7156, Nigeria
4 Department of Mathematics, Texas A M University-Kingsville, Kingsville, TX 78363, USA;
Ravi.Agarwal@tamuk.edu
* Correspondence: ayusuf@biruni.edu.tr

Abstract: In this paper, the wastewater treatment model is investigated by means of one of the most
robust fractional derivatives, namely, the Caputo fractional derivative. The growth rate is assumed
to obey the Contois model, which is often used to model the growth of biomass in wastewaters. The
characteristics of the model under consideration are derived and evaluated, such as equilibrium,
stability analysis, and steady-state solutions. Further, important characteristics of the fractional
wastewater model allow us to understand the dynamics of the model in detail. To this end, we
discuss several important analyses of the fractional variant of the model under consideration. To
observe the efficiency of the non-local fractional differential operator of Caputo over its counter-
classical version, we perform numerical simulations.





 Keywords: wastewater model; stability analysis; Caputo fractional operator; existense and unique-
ness results; numerical dynamics
Citation: Alqahtani, R.T.; Yusuf, A.;
Agarwal, R.P. Mathematical Analysis
of Oxygen Uptake Rate in
Continuous Process under Caputo
Derivative. Mathematics 2021, 9, 675. 1. Introduction
https://doi.org/10.3390/math9060675 Mathematical models as routes of interpreting real-life situations are included in the
literature. These models display a vital update in the quantification and assessment of
Academic Editor: Mariano Torrisi real-life problems and preventive measures [1–3]. Mathematical modeling has, in many
ways, proven to be a very versatile and effective way of studying the dynamics of many
Received: 25 February 2021
situations. To build and test control measures that are persuasive, statistical analysis and
Accepted: 12 March 2021
numerical simulations may be used. It is also widely accepted that by using mathematical
Published: 22 March 2021
models, the prediction to control and solve real-life problems can be reached. A relevant
argument here is that it is important to obtain acceptable criteria for the particular situation
Publisher’s Note: MDPI stays neutral
under consideration and to use these criteria to determine the effects of potential control
with regard to jurisdictional claims in
measures. The dilemma, then, is how, from an optimal point of view, to implement
published maps and institutional affil-
such measures. Several noteworthy attempts have recently been made to establish results
iations.
through mathematical modelling [4–8].
One of the key areas to see the application of mathematical models is the field of
biochemical engineering, where oxygen is a key substrate in aerobic bioprocesses that is
used for growth, maintenance, and other metabolic routes, including product synthesis [9].
Copyright: © 2021 by the authors.
Oxygen is continuously supplied via the gas phase due to its low solubility in broths,
Licensee MDPI, Basel, Switzerland.
which are typically aqueous solutions, and thus, knowledge of the oxygen transfer rate
This article is an open access article
is needed for the design and scale-up of the bioreactor, which can be addressed by a
distributed under the terms and
mathematical modelling approach. Both the transition of oxygen mass from the phase of
conditions of the Creative Commons
liquid to gas and its microorganism intake are of decisive importance, as the transfer of
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
oxygen is the control stage for microbial growth in many processes and may influence
4.0/).
the evolution of bioprocesses. Early in the history of biochemical engineering, this fact

Mathematics 2021, 9, 675. https://doi.org/10.3390/math9060675 https://www.mdpi.com/journal/mathematics

Mathematics 2021, 9, 675 2 of 19

was established and reviewed [10–14]. The oxygen uptake rate is one of the essential
physiological features of culture growth, and has been used to optimize the fermentation
method. The real oxygen uptake rate is usually estimated from the experimentally defined
oxygen uptake rate [15,16].
We consider a wastewater treatment model with the Contois model for the growth
rate. As the main feeds into a tank, the aerator is named [17–19], and the usage theory
can be established. At this point, by biological oxidation of the substrate, one bacteria
population degrades the pollutant. By absorbing oxygen, the reaction produces an aerobic
environment. The mixture is transmitted to a settler tank in the second stage. Here, because
of gravity, the solid pieces settle and collect at the bottom. Part of the sludge inside the
aerator is recycled to promote oxidation. The topic of waste water treatment has also been
studied from many points of view, such as complex observation [20,21] and control [22,23],
although the anaerobic component is considered in the majority of situations [24–26].
However, this model has never been studied using the contemporary fractional op-
erators with the Contois model. Additionally, it has been multifariously proved that
natural phenomena can be more effectively enunciated by fractional models than by
differential integer-order equations. The fractional calculus has taken on the value and pop-
ularity of modeling realistic cases, in particular those with memory effects, because of this
advantage [27–30]. Fractional operators have been used in so many fields where [31–35] is
used. We were inspired to study and evaluate a new fractional version of Caputo’s operator
model given this value. This is the first time that the Caputo operator has been hired for
the model being considered, to the best of our knowledge.

2. Formation of the Model

We present the model to be investigated in this research here. The model with dimen-
sional form is given by:

dS VXµ(S, X )
V = F ( S0 − S ) − ,
dt Yx/s
(1)
dX
V = − FβX + VX µ(S, X ) − Kd VX.
dt
The oxygen accumulation in the liquid phase can be expressed as

VdC Xµ(S, X )
= − FβC + VK I a(C0 − C ) − V , (2)
dt Yx/O2

where
- µ(S, X ) is the Contois model for the specific growth rate, which is expressed as
 
µ(S, X ) = µm Ks XS+S ;
- The rate of oxygen transfer from the gas to the liquid phase is provided by

OTR = K Ia (C0 − C ); (3)

- Yx/O2 is the overall yield of cells from oxygen; and

- K Ia is the oxygen mass transfer coefficient.
A detailed definition of the model can be found in Nomenclature.

The Non-Dimensional Model

Non-dimensionalization is very important in mathematical modeling. Here, we
aim to non-dimensionalize (1) and (2) to study various mathematical properties of the
model. Equations (1) and (2) can be written in non-dimensional form by the following
equations using non-dimensional variable C ∗ = SC0 for oxygen concentration, S∗ = SS0 for
Mathematics 2021, 9, 675 3 of 19

Ks X
the concentration of substrate, X ∗ = S0 as the microorganism concentration, and t∗ = µm t
as time
dS∗ 1 S∗ X ∗
= (1 − S ∗ ) − ∗ ∗ ,
dt τ α (S + X ∗ )
dX ∗ − βX ∗ S∗ X ∗
= ∗ + ∗ − Kd∗ X ∗ , (4)
dt τ (S + X ∗ )
dC ∗ − βC ∗ S∗ X ∗
= ∗ + K ∗Ia (C0∗ − C ∗ ) − .
dt τ γ(S∗ + X ∗ )

In Equation (4), τ ∗ is the main experimental control parameter. It is assumed for the
investigation of the mathematical model that the substrate concentration entering into the
bioreactor is positive S0 > 0 and the rate of change of death is positive Kd∗ > 0. In addition,
it is assumed that 0 < Kd∗ < 1 and
A1: All the parameters are nonnegative.
A2: 0 ≤ µ(S, X ) ≤ 1, µ(0, X ) = 0.
In Equation (4),
Kd
• Kd∗ = µm represents the dimensionless death rate,
Vµ
• τ ∗ = F m represents the dimensionless residence time,
• K ∗Ia = KµIma represents the dimensionless oxygen mass transfer coefficient, and
• C0∗ = CS00 , α∗ = Ks Yx/s , and γ = Ks Yx/O2 are dimensionless parameters.
In the rest of the analysis of Model 4, we drop the star symbols from Equation (4).

3. Analysis of the Model in Classical Sense

In this section, our aim is to analyze the properties of the model in a classical sense.

3.1. Equilibria and Stability Analysis

We consider the number of equilibrium solutions of the model (4). It is obvious that
the model has a branch of the washout given by
 
C0 KIa τ
E0 =(S, X, C ) = 1, 0, . (5)
KIa τ + β

One can obtain the model’s non-free steady-state (4) by setting the right side to zero.
From the model (4), we have,
α
S= ,
τ (1 − K d ) + α − β
(6)
( τ ( K d − 1) + β ) α
X= .
(τ (Kd − 1) − α + β)(Kd τ + β)

Substituting (6) into the third equation of (4), we obtain

C0 KIa τ ( τ ( K d − 1) + β ) α
C= − . (7)
KIa τ + β γ(τ (Kd − 1) − α + β)(KIa τ + β)

Next, we focus on analyzing the stability of the equilibrium of (4). The stability of the
two forms of the equilibrium, E0 and E1 , is studied.

β
Lemma 1. The steady-state solution E0 is locally asymptotically stable when τ < 1− K d , and is
β
unstable when τ > 1− K d .
Mathematics 2021, 9, 675 4 of 19

Proof. The Jacobian matrix of system (4) at E0 is given by

− τ −1 − α −1
 
0
 
β
J ( E0 ) =  0 − τ + 1 − Kd 0 . (8)
 
 
0 − γ −1 −[ τβ + KIa ]

The eigenvalues of (8) is expressed as

− τ −1
 
 
τ (1−Kd )− β
λi =  (9)
 
τ 
 
− KIa ττ + β
β
Note that λ1,3 are negative. After some algebra, we notice λ2 < 0 when τ < 1− K d .

Lemma 2. The solution E1 is locally asymptotically stable for a physically meaningful solution.

Proof. For the stability of the solution E1 , it is possible to write the Jacobian of the system as
 
− J11 − J12 0
 
J ( E1 ) = 
 J21 − J22 0 ,
 (10)
− J31 − J32 − J33

where,

1 X2 S2 X2 XS
J11 = + 2
, J12 = 2
, J21 = 2
, J22 =
τ α (S + X ) α (S + X ) (S + X ) ( S + X )2
X2 S2 KIa τ + β
J31 = ,J =
2 32 2
, J33 = .
γ (S + X ) γ (S + X ) τ

As a simple analysis, J (E1 ) characteristics are expressed as

λ3 + D1 λ2 + D2 λ + D3 , (11)
where
D1 = J33 + J22 + J11 ,
D2 = (J12 + J22 + J33 )J11 + J33 J22 ,
D3 = J33 J11 (J12 + J22 ).

We use the Rough–Hurtwiz criterion, that is, D1 > 0, D3 > 0, and D1 D2 − D3 > 0, to
illustrate the stability of the solution E1 . We have

D1 D2 − D3 = D11 J11 2 + D22 J11 + B33 > 0. (12)

D11 = (J12 + J22 + J33 ),
D22 = J22 J12 + J22 2 + 2 J33 J22 + J33 2 ,
D33 = (J33 + J22 )J33 J22 .

Thus, the condition of the Routh–Huewitz theorem for the stability of E1 is

satisfied.

3.2. Numerical Simulations

The model’s steady-state solution is evaluated in this section with respect to parame-
ters using the experimental values in Table 1. Figure 1 shows stable curves of the oxygen
Mathematics 2021, 9, 675 5 of 19

concentration, the substrate concentration, and the microorganism concentration as func-

tions of residence time for different values of the parameter β. Note that the parameter
β represents the configuration of the reactor. β = 1 corresponds to the continuous flow
reactor, while β = 0 corresponds to the membrane reactor.

Table 1. Parameters’ values.

Parameters Values Unit Reference

−1
µmax 0.9297 day [36]
Kd 0.0131 day−1 [36]
Yx/s 0.2116 gg−1 [36]
Ks 0.4818 g/L [36]
−1
K La 108 day [37]
Yx/O2 0.67 gg−1 [38]
S0 200 mg/L [39]
C0 10 mg/L [37]
β 0–1 - Range value of the parameter

Figure 1a shows that the oxygen concentration decreases as a function of residence

time to minimum value, called the minimum residence time. After this minimum value,
the oxygen concentration increases as the value of the residence time increases. Note that
the minimum residence time for β = 1 and β = 0.1 are 1.30025 and 0.201052, respectively.
Thus, the value of parameter β must be taken lower to minimise the oxygen concentration.
In Figure 1b, the substrate concentration decreases as the value of the residence time
increases for different values of β. In Figure 1c, the microorganism concentration increases
to the maximum value, then decreases as the value of the residence time.

(a) The oxygen concentration.

Figure 1. Cont.
Mathematics 2021, 9, 675 6 of 19

(b) The substrate concentration.

(c) The microorganism concentration.

Figure 1. A steady-state diagram of the oxygen concentration, the substrate concentration, and the
microorganism concentration varying the residence time τ for different values of parameter β. The
value of the parameters stated in Table 1. β = 0.1 (Black color), β = 0.5 (Blue color), β = 1 (Red color).

Figure 2 is an unfolding diagram for the minimum residence time in Figure 1a as a

function of β and the influent concentration (S0 ). We note that the minimum residence time
increases as the value of these parameters increases.
Mathematics 2021, 9, 675 7 of 19

(a)

(b)
Figure 2. Unfolding diagram for the minimum residence time in Figure 1a as a function of β and the
influent concentration (S0 ). The value of the parameters stated in Table 1. (a) Minimum residence
time versus parameter S0 . (b) Minimum residence time versus parameter S0 .
Mathematics 2021, 9, 675 8 of 19

4. Analysis of the Model under Caputo Fractional Operator

Here, we aim to extend the model (4) by incorporating the fractional order in the
sense of a Caputor fractional derivative. The model (4) in the sense of a Caputo fractional
operator is given by

1 SX
CD (1 − S ) − ψ
ψ
S= ,
τ ψ α (S + X )
− βψ X SX ψ
CD − Kd X,
ψ
X= + (13)
τψ (S + X )
− βψ C ψ ψ SX
CD C = + K Ia (C0 − C ) − ψ
ψ
.
τψ γ (S + X )

4.1. Results for Existence of Solutions under Caputo

In this chapter, we address the existence and uniqueness of the suggested model’s solu-
tion via fixed-point theorems. Let’s reformulate the proposed (13) model in the form below:
 ψ
C
 D0+ S = ζ 1 (t, S, X, C ),

C D ψ X = ζ ( t, S, X, C ), (14)
+ 2
C 0ψ

D0+ C = ζ 3 (t, S, X, C ),
with

ζ (t, S, X, C ) = τ1 (1 − S) − α(SSX ,
 1 +X)


βX SX
ζ 2 (t, S, X, C ) = − τ + (S+ X ) − Kd X, (15)

ζ (t, S, X, C ) = − βC + K (C − C ) −
 SX
3 τ Ia 0 γ(S+ X )
.

Now, (13) becomes

( ψ
C D v ( t ) = Λ ( t, v ( t )); t ∈ J = [0, b], 0 < α ≤ 1,
0 (16)
v (0) = v0 ≥ 0,

only if

T
v (t) = (S, X, C ) ,

v (0) = (S0 , X0 , C0 ) T , (17)

Λ(t, v (t)) = (ζ i (t, S, X, C )) T , i = 1, 2, 3,


(·)T implies the transpose operation. Now, (16) becomes

v (t) = v0 + J0υ+ Λ(t, v (t))

1
Z t (18)
= v0 + (t − τ )υ−1 Λ(τ, v (τ ))dτ.
Γ(υ) 0

Let E = C ([0, b]; R) be a Banach space for all the functions which are continuous from
R → [0, b] with the norm kv k = sup |v (t)|.
t∈ J

Theorem 1. Let the function Λ ∈ C ([ J, R]) and maps of the bounded subset of J × R5 be compact
subsets of R. Further, there is a constant LΛ > 0 whereby
(A1 ) |Λ(t, v1 (t)) − Λ(t, v2 (t))| ≤ LΛ |v1 (t) − v2 (t)|; ∀ t ∈ J and all v1 , v2 ∈ C ([J , R]).
Thus, (18), which is equivalent to (13), has a unique solution whenever ΩLΛ < 1, and

bψ
Ω= .
Γ ( ψ + 1)
Mathematics 2021, 9, 675 9 of 19

Proof. Consider that P : E → E is defined by

Z t
1
( Pv )(t) = v0 + (t − τ )ψ−1 Λ(τ, v (τ ))dτ. (19)
Γ(ψ) 0

Now, P is the unique solution of (13) and is well-defined and depicts the fixed-point
of P. Consider supt∈ J kΛ(t, 0)k = M1 and κ ≥ kv0 k + ΩM1 . Now, it suffices to justify
PHκ ⊂ Hκ , and Hκ = {v ∈ E : kv k ≤ κ }, as convex and closed. Thus, each v ∈ Hκ ,
we have
1 t Z
|( Pv )(t)| ≤ |v0 | + (t − τ )ψ−1 |Λ(τ, v (τ ))|dτ
Γ(ψ) 0
Z t
1
≤ v0 + (t − τ )ψ−1 [|Λ(τ, v (τ )) − Λ(τ, 0)| + |Λ(τ, 0)|]dτ
Γ(ψ) 0
(L κ + M1 ) t
Z
≤ v0 + Λ (t − τ )ψ−1 dτ (20)
Γ(ψ) 0
(L κ + M1 ) ψ
≤ v0 + Λ b
Γ ( ψ + 1)
≤ v0 + Ω(LΛ κ + M1 )
≤ κ.

We justify the results. Further, for v1 , v2 ∈ E, one attains

1 t Z
|( Pv1 )(t) − ( Pv2 )(t)| ≤ (t − τ )ψ−1 |Λ(τ, v1 (τ )) − Λ(τ, v2 (τ ))|dτ
Γ(ψ) 0
Z t
L (21)
≤ Λ (t − τ )ψ−1 |v1 (τ ) − v2 (τ )|dτ
Γ(ψ) 0
≤ ΩLΛ |v1 (t) − v2 (t)|.

This justifies that k( Pv1 ) − ( Pv2 )k ≤ ΩLΛ kv1 − v2 k. Thus, due to Banach contrac-
tion, the unique solution for (13) is reached.

Now, we go for the existence of (13) solutions by Krasnoselskii’s fixed-point justification.

Lemma 3. Let M 6= ∅ be a closed, bounded, and convex subset of a Banach Space E. Let two
operators that respect the given relation be P1 , P2 .
• P1 v1 + P2 v2 ∈ M, provided that v1 , v2 ∈ M;
• P1 is compact and continuous;
• P2 is a contraction mapping.
Then, there is u ∈ M as far as u = P1 u + P2 u.

Theorem 2. Surmising Λ : J × R5 → R is continuous and holds for the condition (A1 ).

Further, let
(A2 ) |Λ(t, v )| ≤ ψ(t), for all (t, v ) ∈ J × R5 and ψ ∈ C ([0, b], R+ ).
Thus, Equation (13) has at least one solution whenever

LK kv1 (t0 ) − v2 (t0 )k < 1.

Proof. Setting sup |ψ(t)| = kψk and η ≥ kv0 k + Ωkψk, we consider Bη = {v ∈ E :

t∈ J
kv k ≤ η }. Assume P1 , P2 operators on Bη are expressed as
Z t
1
( P1 v )(t) = (t − τ )ψ−1 Λ(τ, v (τ ))dτ t ∈ J,
Γ(ψ) 0
Mathematics 2021, 9, 675 10 of 19

and
( P2 v )(t) = v (t0 ), t ∈ J.
Now, each v1 , v2 ∈ Bη , gives
Z t
1
k( P1 v1 )(t) + ( P2 v2 )(t)k ≤ kv0 k + (t − τ )ψ−1 kΛ(τ, v1 (τ ))kdτ
Γ(ψ) 0
(22)
≤ k v0 k + Ω k ψ k
≤ η < ∞.

Thus, P1 v1 + P2 v2 ∈ Bη .
Further, the contraction of P2 will be proved.
For t ∈ J and v1 , v2 ∈ Bη , one reaches

k( P2 v1 )(t) − ( P2 v2 )(t)k ≤ kv1 (t0 ) − v2 (t0 )k. (23)

Having Λ as a continuous function, P1 will also be continuous. Moreover, for any

t ∈ J and v1 ∈ Bη ,
k P1 v k ≤ Ωkψk < +∞,
implies that P1 is uniformly bounded. At last, one can justify that P1 is compact. Defining
sup |Λ(t, v (t))| = Λ∗ , yields
(t,v )∈ J ×Bη

Z t
1 1
|( P1 v )(t2 ) − ( P1 v )(t2 )| = [(t2 − τ )ψ−1 − (t1 − τ )ψ−1 ]Λ(τ, v (τ ))dτ
Γ(ψ) 0
Z t2
+ (t2 − τ )ψ−1 Λ(τ, v (τ ))dτ (24)
t1

Λ∗ h ψ ψ
i
≤ 2( t2 − t1 ) ψ + ( t2 − t1 )
Γ(ψ)
→ 0, as t2 → t1 .

This shows that by the Arzelá Ascoli principle, (13) has at least one solution.

4.2. Iterative and Stability Analysis via Caputo Operator

Let (B, ||.||) be a Banach space and Y∗ be a self-map of B. Further, consider yn+1 =
h(Y∗ , yn ) and G(Y∗ ) as a fixed-point set of non-empty Y∗ . It can be interpreted that yn
has converged to q∗ ∈ G(Y∗ ). Defining ||z∗n+1 − h(Y∗ , z∗n )|| with {z∗n ⊆ B}. The iterative
scheme, yn+1 = h(Y∗ , yn ) is Y∗ -stable if limn→∞ cn = 0, that is, limn→∞ c∗n = p∗ . To get
convergence for zn , it has to possess an upper bound. Having satisfied the aforementioned
conditions for yn+1 = Y∗ with n representing Picard’s iteration as in [40], then the iteration
is Y∗ -stable. Thus:

Theorem 3 ([40]). Let (B, ||.||) be a Banach space and Y∗ be a self-map on B. Then, for all
x, y ∈ B, the inequality as comes next holds

kY∗x − Y∗y k≤ K k x − Y∗x k+kk x − yk, (25)

with K ≥ 0, k ∈ [0, 1). If we assume that Y∗ is Picard Y∗ -stable, we can write

Mathematics 2021, 9, 675 11 of 19

  
1 1 Sn Xn
S n +1 ( t ) = S n ( t ) + L −1 L ( 1 − S n ) − ,
sα τψ α ψ ( Xn + Sn )
  
−1 1 β ψ Xn Sn Xn ψ
X n +1 ( t ) = X n ( t ) + L L − ψ + − K d Xn , (26)
sα τ ( Xn + Sn )
  
−1 1 βψ Cn ψ Sn Xn
Cn+1 (t) = Cn (t) + L L − ψ + K Id (C0 − Cn ) − ψ .
sα τ γ ( Xn + Sn )

Theorem 4. Let F be a self-map. Then, it is defined as follows:

  
1 1 Sn Xn
F [Sn (t)] = Sn+1 (t) = Sn (t) + L−1 L ( 1 − S n ) − ,
sα τψ α ψ ( Xn + Sn )
  
−1 1 β ψ Xn Sn Xn ψ
F [ Xn (t)] = Xn+1 (t) = Xn (t) + L L − ψ + − K d Xn , (27)
sα τ ( Xn + Sn )
  
1 βψ Cn Sn Xn
F [Cn (t)] = Cn+1 (t) = Cn (t) + L−1 α L − ψ + K Id (C0 − Cn ) − ψ
ψ
,
s τ γ ( Xn + Sn )

that is F -stable in the space of L1 ( a, b), only if

 !
 1 κ 1 ( ρ ) κ 2 ( ρ )
1 + ψ (1 − κ1 (ρ)) − ψ < 1,






 τ α ( P + Q)

 !

 β ψ κ2 ( ρ ) κ1 ( ρ )κ2 ( ρ ) ψ
1− + − Kd κ2 (ρ) < 1, (28)

 τψ ( P + Q)



 !

 βψ κ3 ( ρ ) ψ κ1 ( ρ )κ2 ( ρ )
 1− + K Id (C0 − κ3 (ρ)) − ψ < 1.


τψ γ ( P + Q)


Proof. Since F is a fixed point, thus, for each (m, n) ∈ N × N, one attains
  
1 1 Sn Xn
F [Sn (t)] − F [Sm (t)] = Sn+1 (t) = Sn (t) + L−1 α L ψ (1 − Sn ) − ψ
s τ α ( Xn + Sn )
  
1 1 S X
m m
− L −1 α L ψ (1 − S m ) − ψ ,
s τ α ( Xm + Sm )
  
1 β ψ Xn Sn Xn
F [ Xn (t)] − F [ Xm (t)] = Xn+1 (t) = Xn (t) + L−1 α L − ψ +
ψ
− K d Xn
s τ ( Xn + Sn )
  ψX  (29)
1 β m S m X m
− L −1 α L − ψ +
ψ
− K d Xm ,
s τ ( Xm + Sm )
  
1 βψ Cn Sn Xn
F [Cn (t)] − F [Cm (t)] = Cn+1 (t) = Cn (t) + L−1 α L − ψ + K Id (C0 − Cn ) − ψ
ψ
s τ γ ( Xn + Sn )
  ψC 
1 β m S m X m
− L−1 α L − ψ + K Id (C0 − Cm ) − ψ
ψ
.
s τ γ ( Xn + Sn )

Imposing norms on both sides of the first equation in (29), we attain

kF (Sn (t)) − F (Sm (t))k= Sn (t) − Sm (t)

  
1 1 Sn Xn
+ L −1 α L ψ (1 − S n ) − ψ (30)
s τ α ( Xn + Sn )
  
1 1 Sm Xm
− L −1 α L ψ (1 − S m ) − ψ , ,
s τ α ( Xm + Sm )
Mathematics 2021, 9, 675 12 of 19

and imposing a triangular inequality gives

kF (Sn (t)) − F (Sm (t))k= kSn (t) − Sm (t)k

  
1 1 Sn Xn
+ L −1 α L ψ (1 − S n ) − ψ
s τ α ( Xn + Sn ) (31)
  
−1 1 1 Sm Xm
−L L ψ (1 − Sm ) − ψ , .
sα τ α ( Xm + Sm )

Thus, (31) becomes

kF (Sn (t)) − F (Sm (t))k≤ kSn (t) − Sm (t)k

  
−1 1 1 Sn Xn
+L L k ψ (1 − Sn )k−k ψ k
sα τ α ( Xn + Sn ) (32)
  
−1 1 1 Sm Xm
−L L k ψ (1 − Sm )k−k ψ k .
sα τ α ( Xm + Sm )

To achieve our goal, we assume that

k Xn (t) − Xm (t)k∼
= kSn (t) − Sm (t)k,
(33)
∼
kCn (t) − Cm (t)k= kSn (t) − Sm (t)k.

Inserting (33) into the relation (32), one can have

kF (Sn (t)) − F (Sm (t))k≤ kSn (t) − Sm (t)k

  
1 1 Sn Xn
+ L−1 α L k ψ (1 − Sn )k−k ψ k
s τ α ( Xn + Sn ) (34)
  
1 1 Sm Xm
− L−1 α L k ψ (1 − Sm )k−k ψ k .
s τ α ( Xm + Sm )

Since Sm (t), Xm (t), and Cm (t) get convergent and bounded, there exists P > 0, Q > 0,
and R > 0 for every t. Thereby, we get

kSm (t)k< P, k Xm (t)k< Q, kCm (t)k< R, (m, n) ∈ N × N. (35)

From (34) and (35), one reaches

!
1 κ1 ( ρ )κ2 ( ρ )
kF (Sn (t)) − F (Sm (t))k≤ 1 + ψ (1 − κ1 (ρ)) − ψ k(Sn (t) − Sm (t))k, (36)
τ α ( P + Q)

with κ1 , κ2 , and κ3 being acquired functions by the inverse Laplace transform in (34). In a
similar way, we achieve
!
βψ κ2 ( ρ ) κ1 ( ρ )κ2 ( ρ ) ψ
kF ( Qn (t)) − F ( Qm (t))k≤ 1− + − Kd κ2 (ρ) k Qn (t) − Qm (t)k, (37)
τψ ( P + Q)
!
βψ κ3 ( ρ ) ψ κ1 ( ρ )κ2 ( ρ )
kF ( En (t)) − F ( Em (t))k≤ 1− + K Id (C0 − κ3 (ρ)) − ψ k En (t) − Em (t)k, (38)
τψ γ ( P + Q)
Mathematics 2021, 9, 675 13 of 19

where the condition (28) is valid. Now, F has a fixed point. To prove that F holds for
Theorem 3.3, let (37) and (38) hold, and the following relations are satisfied:
 !
 1 κ1 ( ρ )κ2 ( ρ )
1 + ψ (1 − κ1 (ρ)) − ψ < 1,






 τ α ( P + Q)

 !

 β ψ κ2 ( ρ ) κ1 ( ρ )κ2 ( ρ ) ψ
k = (0, 0, 0), K = 1− + − Kd κ2 (ρ) < 1, (39)

 τψ ( P + Q)



 !
 β ψ κ (ρ) κ ( ρ ) κ ( ρ )
 3 ψ 1 2
 1− + K Id (C0 − κ3 (ρ)) − ψ < 1.


τψ γ ( P + Q)


Hence, we reach the intended result.

4.3. The Solution’s Positiveness

Let us decide the invariant region and indicate that for all t ≥ 0, all solutions of the
model’s fractional system are positive. The main objective is to implement the convenience
of the analyzed model solutions by observing whether they join the Υ invariant field.
Employing the benefits of the fractional derivative of Caputo, we assume that

Υ = (S, X, C ) ∈ R3+ , R3+ = {n ∈ R3+ : n ≥ 0} (40)

is a model solution with positive initial conditions.
Further, n = (S(t), X (t), C (t)) T . We must also show that on each hyperplane that is
bound by the non-negative hyperoctant, the vector field points to R+3 . Thus, one can write

1
CD
ψ
S(t) = (1 − S) ≥ 0,
τ
SX (41)
ψ
C D Q(t) = ≥ 0,
S+X
C D E ( t ) = K Ia (C0 − C ) ≥ 0.
ψ

Thus, we give the convenient region as follows:

Υ = {(S, X, C ) : S ≥ 0, X ≥ 0, C ≥ 0, S + X + C ≤ 1}. (42)

Thus, in the Υ area, the fractional model is biologically acceptable and mathematically
well-defined when t > 0. In addition, this area is positively invariant, that is, the underlying
system’s solutions are positive for all t.

5. Numerical Dynamics
This is where, with the aid of the best-suited parameters as obtained in the table, we
present different numerical results for the proposed model by using a fractional differential
operator known as the Caputo operator. To serve the purpose of model simulations, a
predictor-corrector numerical technique for the fractional type of dynamical systems as
implemented and analyzed in [41–43] has been used on an effective software called Matlab
2019a. In terms of the Caputo differential operator of order ψ, the following form of the
Cauchy ordinary differential equation is taken into consideration:

C ψ p
0 Dt ϕ ( t ) = v (t, ϕ(t)), ϕ ( p ) (0) = ϕ0 , 0 < ψ ≤ 1, 0 < t ≤ τ, (43)

with p = 0, 1, . . . , n − 1, n = dψe. The Volterra equation version of the Equation (43)

becomes:
n −1 p Z t
( p) t 1
ϕ ( t ) = ∑ ϕ0 + (t − s)ψ−1 v (s, ϕ(s))ds. (44)
p =0
p! Γ ( ψ ) 0
Mathematics 2021, 9, 675 14 of 19

Applying the described predictor-corrector related with the Adam-Bashforth-Moulton

solver [42], we assume h = τ/N, t j = jh, and j = 0, 1, . . . , N ∈ Z+ , by taking ϕ j ≈ ϕ t j ,

with ϕ j standing for the approximate solution and ϕ(t j ) for the true one. We can now write

k −1 j k
!
( j ) t k +1 hψ Sj Xj
 1 
S k +1 = ∑ S0
j!
+ ∑ q j,k+1 τ (1 − Sj ) − α(Sj + Xj )
Γ ( ψ + 2)
j =0 j =0
k p p !
h ψ 1 Sk +1 X k +1
Γ(ψ + 2) j∑
p
+ (qk+1,k+1 ) (1 − S k +1 ) − p p ,
=0
τ α ( S k +1 + X k +1 )
k −1 j k
!
( j) t hψ βX j Sj Xj
 
X k +1 = ∑ X0 k + 1
j!
+
Γ ( ψ + 2) ∑ q j,k+1 −
τ
+
(S j + X j )
− Kd X j
j =0 j =0
k p p p !
hψ βXk+1 S k +1 X k +1
∑ (qk+1,k+1 )
p
+ − + p p − K d X k +1 ,
Γ ( ψ + 2) j =0
τ ( Sk +1 + X k +1 )
k −1 j k
!
( j) t hψ βCj Sj Xj
 
Ck+1 = ∑ C0 k+1
j!
+
Γ ( ψ + 2) ∑ q j,k+1 −
τ
+ K Ia (C0 − Cj ) −
γ(S j + X j )
j =0 j =0
k p p p !
hψ βCk+1 S k +1 X k +1
∑ (qk+1,k+1 )
p
+ − + K Ia (C0 − Ck+1 ) − p p ,
Γ ( ψ + 2) j =0
τ γ ( S k +1 + X k +1 )

where

 ψ +1
 k − (k − ψ)(k + 1)ψ , if j = 0,
q j,k+1 = ( k − j + 2 ) ψ +1 + ( k − j ) ψ +1 − 2 ( k − j + 1 ) ψ +1 , if 1 ≤ j ≤ k, (45)

1, if j = k + 1.

Here, ( p) implies the predictor portion and it can be determined as:

k −1 j k
!
( j) t hψ Sj Xj
  1
S k +1 = ∑ S0 k + 1
j!
+
Γ ( ψ + 2) ∑ w j,k+1
τ
(1 − S j ) −
α(S j + X j )
j =0 j =0

k −1 j k
!
( j) t hψ βX j Sj Xj
 
X k +1 = ∑ X0 k + 1
j!
+
Γ ( ψ + 2) ∑ w j,k+1 −
τ
+
(S j + X j )
− Kd X j
j =0 j =0

k −1 j k
!
( j ) t k +1 hψ βCj Sj Xj
 
Ck+1 = ∑ C0
j!
+
Γ ( ψ + 2) ∑ w j,k+1 −
τ
+ K Ia (C0 − Cj ) −
γ(S j + X j )
j =0 j =0

where
w j,k+1 = (k + 1 − j)ψ − (k − j)ψ .
The numerical simulations of the Caputo model (13) using the parameters mentioned
in Table 2 are presented in this section. In order to solve the nonlinear model shown in
Equation (13), we used the predictor-corrector Adams technique and obtained graphical
results based on the given parameters values. The concentration of substrate S(t), of
oxygen C (t), and the microorganism concentration X (t) are investigated using different
values of the fractional order ψ in order to see the dynamical behavior of the Caputo
model (13) for all three compartments. It can be observed in Figure 3 that the behaviour of
the curves changes as the value of fractional order increases. Moreover, the values of one of
the most important and effective parameters in the model, which is the main experimental
control parameter τ, have been varied in Figure 4 to practically see its influence in the
model. Comparing Figures 3 and 4, the changes in the curves can be noticed as a result of
the effect of the main experimental control, τ.
Mathematics 2021, 9, 675 15 of 19

Table 2. Parameters’ values.

Parameters Values
Kd 0.0140905
α 0.10194888
K Ia 116.1665
γ 0.322806
1
C0 20
β 0–1

0.3
=1
=0.8
0.25 =0.6

0.2
S(t)

0.15

0.1

0.05

0
0 0.5 1 1.5 2 2.5 3 3.5 4
t

(a)
0.22
=1
=0.8
0.21 =0.6

0.2
X(t)

0.19

0.18

0.17

0.16
0 0.5 1 1.5 2 2.5 3 3.5 4
t

(b)

0.1
=1
=0.8
0.09 =0.6

0.08
C(t)

0.07

0.06

0.05

0.04
0 0.5 1 1.5 2 2.5 3 3.5 4
t

(c)
Figure 3. Numerical simulation of the oxygen saturation concentration, the substrate concentration,
and the microorganism concentration for different values of fractional order. The value of the
parameters is stated in Table 1. (a) Dynamic of the concentration of substrate. (b) Dynamic of the
substrate concentration. (c) Dynamic of the microorganism concentration.
Mathematics 2021, 9, 675 16 of 19

0.3
=8
=5
0.25 =2

0.2

S(t)
0.15

0.1

0.05

0
0 0.5 1 1.5 2 2.5 3 3.5 4
t

(a)
0.22
=8
=5
0.2 =2

0.18
X(t)

0.16

0.14

0.12

0.1
0 0.5 1 1.5 2 2.5 3 3.5 4
t

(b)
0.1
=8
=5
0.09 =2

0.08

0.07
C(t)

0.06

0.05

0.04

0.03
0 0.5 1 1.5 2 2.5 3 3.5 4
t

(c)
Figure 4. Numerical simulation of the oxygen saturation concentration, the substrate concentration,
and the microorganism concentration for different values of residence time. The value of the pa-
rameters is stated in Table 1. (a) Dynamic of the concentration of substrate with different values
of τ. (b) Dynamic of the substrate concentration with different values of τ. (c) Dynamic of the
microorganism concentration with different values of τ.

6. Conclusions
In the current analysis, the wastewater model has been studied from a classical
viewpoint and modelled by one of the most robust non-local fractional operators, called the
Caputo operator. Firstly, some important properties of the model in a classical sense were
reported. The dynamical features of the model were depicted through the steady-state
region in its classical sense. Since it has been proved in the literature that the fractional
operators are especially suitable for studying the dynamical behaviour of real-world
problems, it is of vehement importance to analyze more critically the dynamics of the model
Mathematics 2021, 9, 675 17 of 19

under consideration with the robust fractional operator called Caputo, and also make some
comparisons with the classical variant. The fractionalized order is ψ. As a consequence,
many important characteristics of the proposed fractional version of the model have been
reported, such as the formation of the model, the existence and uniqueness of the solution
through the fixed-point theorem, iterative solutions, and the solutions’ positivity. It should
be observed that the fractional version of the model under investigation more correctly
understands the behavior of the model than the variant of the integer order. This has been
done by numerous numerical simulations by using an efficient numerical scheme.

Author Contributions: Conceptualization, A.Y.; Data curation, R.T.A. and A.Y.; Formal analysis,
R.T.A. and A.Y.; Funding acquisition, R.T.A.; Investigation, R.-P.A.; Methodology, R.T.A. and A.Y.;
Visualization, R.-P. A.; Writing–original draft, R.-P.A.; Writing–review–editing, R.-P.A. All authors
have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: All data used in this work can be found within the manuscript.
Conflicts of Interest: The authors declare no conflict of interest.

Nomenclature
Symbol Description
F The flow rate through the reactor.
S The substrate concentration within the reactor.
S∗ The dimensionless substrate concentration within the reactor [S∗ = S/S0 ].
C The dissolved oxygen concentration .
C∗ - The dimensionless dissolved oxygen concentration SC0 .
KI a Volumetric oxygen mass transfer coefficient
C0 The dissolved oxygen solubility limit
C0∗ The dimensionless dissolved oxygen solubility limit
Kd The death coefficient
Kd∗ The dimensionless death coefficient
Ks The saturation constant
S0 The concentration of substrate flowing into the reactor.
V The volume of the reactor.
Yxls The yield factor for substrate
α∗ The dimensionless yield factor for substrate
Yxlo2 The yield coefficient for oxygen
γ∗ The dimensionless yield coefficient for oxygen.
µm The specific growth rate parameter
β The reactor model for the reactor.
X The biomass concentration within the reactor
X∗ The dimensionless biomass concentration within the reactor
t The time.
t∗ The dimensionless time.
µ(S, X ) The specific growth rate model.
τ The residence time.
τ∗ The dimensionless residence time.

References
1. Brauer, F.; Castillo-Chavez, C.; Feng, Z. Mathematical Models in Epidemiology; Springer: New York, NY, USA, 2019.
2. Qureshi, S.; Yusuf, A. Modeling chickenpox disease with fractional derivatives: From caputo to atangana-baleanu. Chaos Solitons
Fractals 2019, 122, 111–118. [CrossRef]
Mathematics 2021, 9, 675 18 of 19

3. Qureshi, S.; Yusuf, A. Mathematical modeling for the impacts of deforestation on wildlife species using Caputo differential
operator. Chaos Solitons Fractals 2019, 126, 32–40. [CrossRef]
4. Jajarmi, A.; Yusuf, A.; Baleanu, D.; Inc, M. A new fractional HRSV model and its optimal control: A non-singular operator
approach. Phys. A 2020, 547, 123860. [CrossRef]
5. Abdulhameed, M.; Muhammad, M.M.; Gital, A.Y.; Yakubu, D.G.; Khan, I. Effect of fractional derivatives on transient MHD flow
and radiative heat transfer in a micro-parallel channel at high zeta potentials. Phys. A 2019, 519, 42–71. [CrossRef]
6. VDubey, P.; Kumar, R.; Kumar, D. Analytical study of fractional Bratu-type equation arising in electro-spun organic nanofibers
elaboration. Phys. A 2019, 521, 762–772.
7. Chang, A.; Sun, H.; Zhang, Y.; Zheng, C.; Min, F. Spatial fractional Darcy’s law to quantify fluid flow in natural reservoirs. Phys.
A 2019, 519, 119–126. [CrossRef]
8. Goulart, A.G.; Lazo, M.J.; Suarez, J.M.S. A new parameterization for the concentration flux using the fractional calculus to model
the dispersion of contaminants in the Planetary Boundary Layer. Phys. A 2019, 518, 38–49. [CrossRef]
9. Ochoaa, F.G.; Gomez, E.; Santos, V.E.; Merchuk, J.C. Oxygen uptake rate in microbial processes: An overview. Biochem. Eng. J.
2010, 49, 289–307. [CrossRef]
10. Garcia-Ochoa, F.; Gomez, E.; Santos, V. Oxygen transfer and uptake rates during xanthan gum production. Enzyme Microb.
Technol. 2000, 27, 680–690. [CrossRef]
11. Garcia-Ochoa, F.; Gomez, E. Mass transfer coefficient in stirrer tank reactors for xanthan solutions. Biochem. Eng. J. 1998, 1, 1–10.
[CrossRef]
12. Ho, C.; Oldshue, J.Y. (Eds.) Biotechnology Processes: Scale-Up and Mixing; AIChE: New York, NY, USA, 1987.
13. Andrew, S.P.S. Gas–liquid mass transfer in microbiological reactors. Trans. IChemE 1982, 60, 3–13.
14. Aiba, S.; Humphrey, A.E.; Mills, N.F. (Eds.) Biochemical Engineering, 2nd ed.; Academic Press: New York, NY, USA, 1973.
15. Zou, X.; Hang, H.; Chu, J.; Zhuang, Y.; Zhang, S. Oxygen uptake rate optimization with nitrogen regulation for erythromycin
production and scale-up from 50 L to 372 m3 scale. Bioresour. Technol. 2009, 100, 1406–1412. [CrossRef] [PubMed]
16. Palomares, L.A.; Ramirez, O. The effect of dissolved oxygen tension and the utility of oxygen uptake rate in insect cell culture.
Cytotechnology 1996, 22, 225–237. [CrossRef] [PubMed]
17. Bhattacharya, S.C.; Khai, P.Q. Kinetics of Anaerobic Cowdung Digestion. Energy 1987, 12, 497–500. [CrossRef]
18. Serhani, M.; Gouze, J.L.; Raissi, N. Dynamical study and robustness of a nonlinear wastewater treatment problem. J. Nonlinear
Anal. RWA 2011, 12, 487–500. [CrossRef]
19. Haandel, A.V.; Lubbe, J.V. Handbook Biological Waste Water Treatment; IWA Publishing: Leidschendam, The Netherlands, 2007.
20. Alqahtani, T.R.; Nelson, I.M.; Worthy, L.A. Analysis of a chemostat model with variable yield coefficient: Contois kinetics.
ANZIAM J. 2012, 53, C155–C171. [CrossRef]
21. Alqahtani, T.R.; Nelson, I.M.; Worthy, L.A. A fundamental analysis of continuous flow bioreactor models with recycle around
each reactor governed by Contois kinetics. III. Two and three reactor cascades. Chem. Eng. J. 2012, 183, 422–432. [CrossRef]
22. Koumboulis, F.N.; Kouvakas, N.D.; King, R.E.; Stathaki, A. Two-stage robust control of substrate concentration for an activated
sludge process. ISA Trans. 2008, 47, 267–278. [CrossRef]
23. Jourani, A.; Serhani, M.; Boutoulout, A. Dynamic and controllability of a nonlinear wastewater treatment problem. J. Appl. Math.
Inform. 2012, 30, 883–902.
24. Khan, K.; Zarin, R.; Khan, A.; Yusuf, A.; Al-Shomrani, M.; Ullah, A. Stability analysis of five-grade Leishmania epidemic model
with harmonic mean-type incidence rate. Adv. Differ. Equ. 2021, 86. [CrossRef]
25. Serhani, M.; Cartigny, P.; Raissi, N. Robust feedback design of wastewater treatment problem. Math. Model. Nat. Phenom. 2009,
4, 1139–1143. [CrossRef]
26. Haegeman, B.; Lobry, C.; Harmand, J. Modeling bacteria flocculation as density-dependent growth. AIChE J. 2007, 53, 535–539.
[CrossRef]
27. Singh, J.; Kumar, D.; Baleanu, D. On the analysis of fractional diabetes model with exponential law. Adv. Differ. Equ. 2018,
2018, 231. [CrossRef]
28. Yusuf, A.; Acay, B.; Mustafa, U.T.; Inc, M.; Baleanu, D. Mathematical modeling of pine wilt disease with Caputo fractional
operatorChaos. Solitons Fractals 2021, 143, 110569. [CrossRef]
29. Musa, S.S.; Qureshi, S.; Zhao, S.; Yusuf, A.; Mustapha, U.T.; He, D. Mathematical modeling of COVID-19 epidemic with effect of
awareness programs. Infect. Dis. Model. 2021, 6, 448-460.
30. Baba, I.A.; Yusuf, A.; Al-Shomrani, M. A mathematical model for studying rape and its possible mode of control. Results Phys.
2021, 22, 103917. [CrossRef]
31. Acay, B.; Inc, M. Fractional modeling of temperature dynamics of a building with singular kernels. Chaos Solitons Fractals 2021,
142, 110482. [CrossRef]
32. Ahmed, I.; Goufo, E.F.D.; Yusuf, A.; Kumam, P.; Chaipanya, P.; Nonlaopon, K. An epidemic prediction from analysis of a
combined HIV-COVID-19 co-infection model via ABC fractional Operator. Alex. Eng. J. 2021, 60, 2979–2995 [CrossRef]
33. Acay, B.; Bas, E.; Abdeljawad, T. Non-local fractional calculus from different viewpoint generated by truncated M-derivative. J.
Comput. Appl. Math. 2020, 366, 112410. [CrossRef]
34. Kirtphaiboon, S.; Humphries, U.; Khan, A.; Yusuf, A. Model of rice blast disease under tropical climate conditions. Chaos Solitons
Fractals 2021, 143, 110530. [CrossRef]
Mathematics 2021, 9, 675 19 of 19

35. Noeiaghdam, S.; Sidorov, D.N. Caputo-Fabrizio Fractional Derivative to Solve the Fractional Model of Energy Supply-Demand
System. Math. Model. Eng. Probl. 2020, 7, 359–367. [CrossRef]
36. Hu, W.C.; Thayanithy, K.; Forster, C.F. A kinetic study of the anaerobic digestion of ice-cream wastewater. Process. Biochem. 2002,
37, 965–971. [CrossRef]
37. Fikar, M.; Chachuat, B.; Latifi, M.A. Optimal operation of alternating activated sludge processes. Control Eng. Pract. 2005, 13,
853–861. [CrossRef]
38. Henze, M.; Grady, C.P.L., Jr.; Gujer, W.; Marais, G.V.R.; Matsuo, T. A general model for single-sludge wastewater treatment
systems. Water Res. 1987, 21, 505–515. doi:10.1016/0043-1354(87)90058-3. [CrossRef]
39. Yoon, S.-H.; Lee, S. Critical operational parameters for zero sludge production in biological wastewater treatment processes
combined with sludge disintegration. Water Res. 2005, 39, 3738–3754. [CrossRef]
40. Wang, Z.; Yang, D.; Ma, T.; Sun, N. Stability analysis for nonlinear fractional-order systems based on comparison principle.
Nonlinear Dyn. 2014, 75, 387–402. [CrossRef]
41. Diethelma, K.; Ford, N.J.; Freed, A.D. A Predictor-Corrector Approach for the Numerical Solution of Fractional Differential
Equations. Nonlinear Dyn. 2002, 29, 3–22. [CrossRef]
42. Diethelm, K. An algorithm for the numerical solution of differential equations of fractional order. Electron. Trans. Numer. Anal.
1997, 5, 1–6. [CrossRef]
43. Diethelm, K.; Ford, N.J.; Freed, A.D. Detailed error analysis for a fractional Adams method. Numer. Algorithms 2004, 36, 31–52.