# Dynamics Analysis and Biomass Productivity Optimisation of a Microbial Cultivation Process through Substrate Regulation

Authors: Kaibiao Sun; Shan Liu; Andrzej Kasperski; Yuan Tian
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## Abstract

A microbial cultivation process model with variable biomass yield, control of substrate concentration, and biomass recycle is formulated, where the biochemical kinetics follows an extension of the Monod and Contois models. Control of substrate concentration allows for indirect monitoring of biomass and dissolved oxygen concentrations and consequently obtaining high yield and productivity of biomass. Dynamics analysis of the proposed model is carried out and the existence of order-1 periodic solution is deduced with a formulation of the period, which provides a theoretical possibility to convert the state-dependent control to a periodic one while keeping the dynamics unchanged. Moreover, the stability of the order-1 periodic solution is verified by a geometric method. The stability ensures a certain robustness of the adopted control; that is, even with an inaccurately detected substrate concentration or a deviation, the system will be always stable at the order-1 periodic solution under the control. The simulations are carried out to complement the theoretical results and optimisation of the biomass productivity is presented.


---

## Body

## 1. Introduction

Microorganisms play an important role in nature and their activities have numerous industrial applications [1]. For that reason, bioreactor engineering, as a branch of chemical engineering and biotechnology, is an active area of research on microbial cultivation process, including development, control, and commercialization of new technology [2]. Reaching optimal results and attaining maximal profits require modern control strategies based on mathematical models or artificial intelligence methods. Many dynamic models concerning microbial cultivation processes, employing several types of reactions and control technologies, have been established [3–6]. In microbial cultivation process, there are a lot of factors affecting the growth and reproduction of the microorganisms. For example, for some aerobic microorganisms, the dissolved oxygen content in the bioreactor medium is a key factor to microbial growth. In order to maintain the dissolved oxygen content in an appropriate range, it is necessary to monitor and control the dissolved oxygen concentration (DOC) in the bioreactor medium since a low level of DOC decreases biomass yield and specific growth rate.For a given microorganism, the oxygen demand is affected by several factors, for example, the biomass and substrate concentrations and cell metabolism state. The control of biomass and substrate concentrations is one of the most important ways to maintain an appropriate dissolved oxygen concentration. To obtain this goal, several approaches to bioprocess design have been presented in recent years. Especially attractive is the use of the self-cycling bioprocess. In such a bioprocess, impulsive controls are introduced once the monitored state reaches an established control level. In fact, many biological phenomena such as bursting rhythm models in biology and pharmacokinetics or frequency modulated systems exhibit impulsive effects, and periodic and state-dependent impulsive effects are two commonly occurring ones. In recent years, state-dependent impulsive control strategies have been used in the study on the interaction between wild and transgenic mosquito populations [7], predator-prey systems [8–11], management of pests and fisheries [12–14], pulse vaccination for human infectious diseases [15], and bioprocess models [16–20]. In bioprocess, self-cycling is a new approach to deal with important environmental problems [16], and for that reason it has become a subject of consideration presented in this work. In general, the self-cycling bioprocess is a computer-controlled system, which can be considered as a specific type of continuous bioprocess. The impulsive effect causes the removal of a certain fraction of the medium and the input of an equal volume of fresh medium [16, 21]. The aim of introduction of the impulsive effect is to provide among others appropriate oxygen conditions by limitation of the biomass concentration, which stems from the fact that insufficient dissolved oxygen concentration decreases the biomass yield, specific growth rate, and biomass productivity.In bioprocess control, a key question is how to monitor reactant concentrations in a reliable and cost-effective manner [22]. However, in many industrial applications, not all of concentrations of the components (which are involved and critical for quality control) are available for online measurement. In this work, the substrate concentration is assumed to be measured online [16, 23, 24], and the biomass concentration is controlled indirectly by substrate concentration regulation. The microorganisms grow by assimilation of the substrate; for that reason the biomass concentration increases and substrate concentration decreases. Thus, in order to keep the biomass concentration lower than a critical level, the substrate concentration should be kept higher than a given level. This is also reasonable since a low level of the substrate concentration is disadvantageous to the transformation of microorganisms.Based on above considerations, this work proposes a mathematical model of a microbial cultivation process with variable biomass yield and substrate regulation. The paper is organized as follows. In Section2, a microbial cultivation process model of the considered bioprocess is formulated, the biochemical kinetics following an extension of the Monod and Contois models. In Section 3, the dynamic analysis of the proposed model is carried out, and the existence of order-1 periodic solution is deduced. The complete expression of the period of the order-1 periodic solution is also given. In addition, it is shown that there does not exist any order-k(k
≥
2) periodic solution. Next, a stability criterion for a general semicontinuous dynamic system is presented by a geometric method, and the stability of the order-1 periodic solution is verified. In Section 4, numerical simulations are carried out to complement the theoretical results and optimisation of the biomass productivity is presented. The final conclusions are presented in Section 5.


## 2. Model Formulation and Preliminaries

### 2.1. Model Formulation

A general model to describe the microbial growth on a limiting substrate is of the form(1)
d
s
d
t
=
D
s
f
-
s
-
1
Y
x
/
s
μ
s
,
x
x
d
x
d
t
=
μ
s
,
x
x
-
D
x
s
t
0
=
s
0
,
x
t
0
=
x
0
,where x
(
t
) and s
(
t
) denote, respectively, the concentrations of biomass and substrate in the bioreactor medium at time t (g
/
L), s
f denotes the substrate concentration in the input flow (
g
/
L
), D denotes the dilution rate (h
-
1) (D
>
0: chemostat; D
=
0: batch process), Y
x
/
s is the biomass yield (
g
/
g
) defined as the ratio of the biomass produced to the amount of substrate assimilated, and the function μ
(
s
,
x
) describes the biochemical kinetics, which is characterized by the cell concentration (x) and the limiting substrate concentration (
s
).The key to the proper functioning of a bioreactor lies mainly in understanding the growth rate of the biomass, which plays an important role in the intrinsic oscillation mechanisms [25]. Many kinetics models that have been proposed in the literature, including Tessier [26], Monod [27], Moser [28], Contois [29], and Andrews [30], address the kinetics of the growth of the cell mass in the bioreactor and most of these models report expressions for specific growth rate of cell mass. In this study, an extension of the Monod and Contois models is introduced; that is,(2)
μ
s
,
x
=
μ
M
s
K
ν
x
ν
+
s
,where ν
≥
0 is a constant and K
ν
>
0 is constant in the extended Contois model (g
/
L) and x is dimensionless biomass concentration.The biomass yield expression plays an important role for the generation of oscillatory behaviour in continuous bioreactor models [25]. In many models of bioreactors, it is assumed that the biomass yield coefficient is a constant during the course of the reaction. However, several researchers also suggest that the experimentally found oscillatory behaviour of microbial population can be explained when the biomass yield coefficient is dependent on the substrate concentration [31, 32] and cannot when the biomass yield coefficient is constant. Thus a variable biomass yield is presented as [25, 33–35]:(3)
Y
x
/
s
=
a
+
b
s
,
a
>
0
,
b
≥
0
.To effectively utilize the limiting substrate and avoid unnecessary waste, a modification on the chemostat is made by changing the inflow and outflow from continuous to impulsive. The sketch map of the apparatus is illustrated in Figure1. The apparatus includes an optical sensing device which continuously monitors the substrate concentration in the bioreactor medium and two switches which work in a synchronous way. When the substrate concentration is higher than a substrate low control level, the switches are closed. Once the substrate concentration s
(
t
) in the bioreactor medium decreases to a control level s
C
L (where 0
<
s
m
i
n
≤
s
C
L
≤
s
m
a
x
<
s
f), then part of the medium containing biomass and substrate is discharged from the bioreactor to a membrane filter, and the next portion of medium of a given substrate concentration is inputted impulsively. Therefore, system (1) can be modified as follows by introducing the impulsive state feedback control:(4)
d
s
d
t
=
-
1
a
+
b
s
μ
M
s
K
ν
x
ν
+
s
x
,
d
x
d
t
=
μ
M
s
K
ν
x
ν
+
s
x
,
s
>
s
C
L
,
Δ
s
=
1
-
k
w
f
s
f
-
s
,
Δ
x
=
-
w
f
1
-
1
-
θ
k
x
,
s
=
s
C
L
,
x
t
0
+
=
x
0
,
s
t
0
+
=
s
0
,where s
f is the concentration of the feed substrate which is inputted impulsively, 0
<
w
f
m
i
n
≤
w
f
≤
w
f
m
a
x
<
1 is the proportion of medium which is removed from the bioreactor in each substrate oscillation cycle, 0
<
θ
m
i
n
≤
θ
≤
1 represents the biomass filter efficiency constant, and 0
≤
k
≤
k
m
a
x
<
1 is the substrate recycle part.Figure 1
Schematic view of the impulsive bioprocess.


### 2.2. Preliminaries

Let us consider a general model with impulsive effects:(5)
d
s
d
t
=
F
s
,
x
,
d
x
d
t
=
G
s
,
x
,
i
f
ϕ
s
,
x
≠
0
,
Δ
s
=
α
s
,
x
,
Δ
x
=
β
s
,
x
,
i
f
ϕ
s
,
x
=
0
,where F
(
s
,
x
) and G
(
s
,
x
) are indefinitely differentiable with respect to (
s
,
x
)
∈
Ω
⊂
R
2, and ϕ, α, and β are linearly dependent on s and x, that is, ϕ
s, ϕ
x, α
s, α
x, β
s, and β
x are constant.Denote(6)
M
I
M
P
=
s
,
x
∣
ϕ
s
,
x
=
0
,
N
P
H
A
=
s
,
x
∣
s
=
s
′
+
α
s
′
,
x
′
,
x
=
x
′
+
β
s
′
,
x
′
,
s
′
,
x
′
∈
M
I
M
P
. Then M
I
M
P is referred to as the control set, which contains all the states at which the control strategy is taken on and α and β describe the effects of the control strategy. When the state variables (
s
,
x
) arrive at the control set M
I
M
P, the impulsive control measures are taken; then the state variables (
s
,
x
) jump from M
I
M
P to the set N
P
H
A.Letx
(
t
)
=
(
s
(
t
)
,
x
(
t
)
) be any solution of system (5). Let γ
t
(
x
,
x
0
)
=
x
(
t
), given the initial value x
(
0
)
=
x
0. Then, the (forward) orbit is the set of all values that this trajectory obtains γ
(
x
,
x
0
)
=
{
γ
t
x
,
x
0
∣
t
≥
t
0
}, also denoted as γ
(
x
) in short. Denote x
k
=
x
(
t
k
+
)
∈
γ
(
x
), where t
k
∈
∏
≜
{
t
k
∣
k
=
1,2
,
…
} with x
(
t
k
)
∈
M
I
M
P.LetS
⊆
R
2
=
(
-
∞
,
+
∞
)
2 be an arbitrary set, let P
∈
R
2 be an arbitrary point, and δ
>
0. Then the distance between P and S and the noncentral δ
-neighborhood of P are denoted by(7)
d
P
,
S
=
inf
P
0
∈
S
⁡
P
-
P
0
,
U
δ
o
P
=
Q
∣
0
<
d
P
,
Q
≤
δ
.Definition 1 (periodic orbit [36, 37]).
An orbitγ
(
x
) of system (5) is said to be periodic if there exists a positive integer m
⩾
1 such that x
m
=
x
0. Denote k
=
m
i
n
{
m
∣
x
m
=
x
0
}. Then the orbit γ
(
x
) is said to be an order-k periodic orbit.Definition 2 (ϵ-close [38]).
Two orbitsγ
(
x
) and γ
(
x
^
) are ϵ-close if there is a reparameterization of time (a smooth, monotonic function) t
^
(
t
) such that(8)
γ
t
x
-
γ
t
^
t
x
^
<
ϵ
∀
t
≥
t
0
.Definition 3 (orbitally stable [36–38]).
An orbitγ
(
x
^
) is said to be orbitally stable if, for any ϵ
>
0, there is a neighborhood V of x
^ so that, for all x in V, γ
(
x
) and γ
(
x
^
) are ϵ-close, as illustrated in Figure 2.Figure 2
Orbital stability.Definition 4 (asymptotic orbital stability [36–38]).
An orbitγ
(
x
^
) is said to be asymptotically orbitally stable if it is orbitally stable and additionally V may be chosen such that, for all x
∈
V, there exists a constant τ
(
x
) so that γ
t
x
-
γ
t
-
τ
x
x
^
→
0 as t
→
∞, as illustrated in Figure 3.Figure 3
Asymptotic orbital stability.Definition 5 (successor function [37]).
Suppose that the control setM
I
M
P and the phase set N
P
H
A are both lines, as illustrated in Figure 4. The coordinate on the phase set N
P
H
A is defined as follows: the origin is defined as the intersection point between N
P
H
A and the s-axis, denoted as O
′, and the coordinate of any point P on N
P
H
A is defined as the distance between P and O
′. Let P
3 denote the intersection point between the trajectory starting from P
1 and the control set M
I
M
P, and P
4 denote the phase point of P
3 after control measures. Let P
2 denote the intersection point if the trajectory starting from P
1 intersects the phase set N
P
H
A. Then P
4 is called as the successor point of P
1, and the mapping from P
1 to P
4 constructs the successor function. If P
2 exists, the successor function is called type-II; that is, f
S
O
R
I
I
(
P
1
)
=
d
(
P
4
,
O
′
)
-
d
(
P
1
,
O
′
); otherwise, it is called type-I; that is, f
S
O
R
I
(
P
1
)
=
d
(
P
4
,
O
′
)
-
d
(
P
1
,
O
′
). Moreover, the successor function is continuous on N
P
H
A.Figure 4
Illustration of successor function.Remark 6.
If there exists a pointP
(
s
P
,
x
P
)
∈
N
P
H
A such that f
S
O
R
I
(
P
)
=
0 or f
S
O
R
I
I
(
P
)
=
0, then the orbit starting from P forms a periodic orbit.


## 2.1. Model Formulation

A general model to describe the microbial growth on a limiting substrate is of the form(1)
d
s
d
t
=
D
s
f
-
s
-
1
Y
x
/
s
μ
s
,
x
x
d
x
d
t
=
μ
s
,
x
x
-
D
x
s
t
0
=
s
0
,
x
t
0
=
x
0
,where x
(
t
) and s
(
t
) denote, respectively, the concentrations of biomass and substrate in the bioreactor medium at time t (g
/
L), s
f denotes the substrate concentration in the input flow (
g
/
L
), D denotes the dilution rate (h
-
1) (D
>
0: chemostat; D
=
0: batch process), Y
x
/
s is the biomass yield (
g
/
g
) defined as the ratio of the biomass produced to the amount of substrate assimilated, and the function μ
(
s
,
x
) describes the biochemical kinetics, which is characterized by the cell concentration (x) and the limiting substrate concentration (
s
).The key to the proper functioning of a bioreactor lies mainly in understanding the growth rate of the biomass, which plays an important role in the intrinsic oscillation mechanisms [25]. Many kinetics models that have been proposed in the literature, including Tessier [26], Monod [27], Moser [28], Contois [29], and Andrews [30], address the kinetics of the growth of the cell mass in the bioreactor and most of these models report expressions for specific growth rate of cell mass. In this study, an extension of the Monod and Contois models is introduced; that is,(2)
μ
s
,
x
=
μ
M
s
K
ν
x
ν
+
s
,where ν
≥
0 is a constant and K
ν
>
0 is constant in the extended Contois model (g
/
L) and x is dimensionless biomass concentration.The biomass yield expression plays an important role for the generation of oscillatory behaviour in continuous bioreactor models [25]. In many models of bioreactors, it is assumed that the biomass yield coefficient is a constant during the course of the reaction. However, several researchers also suggest that the experimentally found oscillatory behaviour of microbial population can be explained when the biomass yield coefficient is dependent on the substrate concentration [31, 32] and cannot when the biomass yield coefficient is constant. Thus a variable biomass yield is presented as [25, 33–35]:(3)
Y
x
/
s
=
a
+
b
s
,
a
>
0
,
b
≥
0
.To effectively utilize the limiting substrate and avoid unnecessary waste, a modification on the chemostat is made by changing the inflow and outflow from continuous to impulsive. The sketch map of the apparatus is illustrated in Figure1. The apparatus includes an optical sensing device which continuously monitors the substrate concentration in the bioreactor medium and two switches which work in a synchronous way. When the substrate concentration is higher than a substrate low control level, the switches are closed. Once the substrate concentration s
(
t
) in the bioreactor medium decreases to a control level s
C
L (where 0
<
s
m
i
n
≤
s
C
L
≤
s
m
a
x
<
s
f), then part of the medium containing biomass and substrate is discharged from the bioreactor to a membrane filter, and the next portion of medium of a given substrate concentration is inputted impulsively. Therefore, system (1) can be modified as follows by introducing the impulsive state feedback control:(4)
d
s
d
t
=
-
1
a
+
b
s
μ
M
s
K
ν
x
ν
+
s
x
,
d
x
d
t
=
μ
M
s
K
ν
x
ν
+
s
x
,
s
>
s
C
L
,
Δ
s
=
1
-
k
w
f
s
f
-
s
,
Δ
x
=
-
w
f
1
-
1
-
θ
k
x
,
s
=
s
C
L
,
x
t
0
+
=
x
0
,
s
t
0
+
=
s
0
,where s
f is the concentration of the feed substrate which is inputted impulsively, 0
<
w
f
m
i
n
≤
w
f
≤
w
f
m
a
x
<
1 is the proportion of medium which is removed from the bioreactor in each substrate oscillation cycle, 0
<
θ
m
i
n
≤
θ
≤
1 represents the biomass filter efficiency constant, and 0
≤
k
≤
k
m
a
x
<
1 is the substrate recycle part.Figure 1
Schematic view of the impulsive bioprocess.


## 2.2. Preliminaries

Let us consider a general model with impulsive effects:(5)
d
s
d
t
=
F
s
,
x
,
d
x
d
t
=
G
s
,
x
,
i
f
ϕ
s
,
x
≠
0
,
Δ
s
=
α
s
,
x
,
Δ
x
=
β
s
,
x
,
i
f
ϕ
s
,
x
=
0
,where F
(
s
,
x
) and G
(
s
,
x
) are indefinitely differentiable with respect to (
s
,
x
)
∈
Ω
⊂
R
2, and ϕ, α, and β are linearly dependent on s and x, that is, ϕ
s, ϕ
x, α
s, α
x, β
s, and β
x are constant.Denote(6)
M
I
M
P
=
s
,
x
∣
ϕ
s
,
x
=
0
,
N
P
H
A
=
s
,
x
∣
s
=
s
′
+
α
s
′
,
x
′
,
x
=
x
′
+
β
s
′
,
x
′
,
s
′
,
x
′
∈
M
I
M
P
. Then M
I
M
P is referred to as the control set, which contains all the states at which the control strategy is taken on and α and β describe the effects of the control strategy. When the state variables (
s
,
x
) arrive at the control set M
I
M
P, the impulsive control measures are taken; then the state variables (
s
,
x
) jump from M
I
M
P to the set N
P
H
A.Letx
(
t
)
=
(
s
(
t
)
,
x
(
t
)
) be any solution of system (5). Let γ
t
(
x
,
x
0
)
=
x
(
t
), given the initial value x
(
0
)
=
x
0. Then, the (forward) orbit is the set of all values that this trajectory obtains γ
(
x
,
x
0
)
=
{
γ
t
x
,
x
0
∣
t
≥
t
0
}, also denoted as γ
(
x
) in short. Denote x
k
=
x
(
t
k
+
)
∈
γ
(
x
), where t
k
∈
∏
≜
{
t
k
∣
k
=
1,2
,
…
} with x
(
t
k
)
∈
M
I
M
P.LetS
⊆
R
2
=
(
-
∞
,
+
∞
)
2 be an arbitrary set, let P
∈
R
2 be an arbitrary point, and δ
>
0. Then the distance between P and S and the noncentral δ
-neighborhood of P are denoted by(7)
d
P
,
S
=
inf
P
0
∈
S
⁡
P
-
P
0
,
U
δ
o
P
=
Q
∣
0
<
d
P
,
Q
≤
δ
.Definition 1 (periodic orbit [36, 37]).
An orbitγ
(
x
) of system (5) is said to be periodic if there exists a positive integer m
⩾
1 such that x
m
=
x
0. Denote k
=
m
i
n
{
m
∣
x
m
=
x
0
}. Then the orbit γ
(
x
) is said to be an order-k periodic orbit.Definition 2 (ϵ-close [38]).
Two orbitsγ
(
x
) and γ
(
x
^
) are ϵ-close if there is a reparameterization of time (a smooth, monotonic function) t
^
(
t
) such that(8)
γ
t
x
-
γ
t
^
t
x
^
<
ϵ
∀
t
≥
t
0
.Definition 3 (orbitally stable [36–38]).
An orbitγ
(
x
^
) is said to be orbitally stable if, for any ϵ
>
0, there is a neighborhood V of x
^ so that, for all x in V, γ
(
x
) and γ
(
x
^
) are ϵ-close, as illustrated in Figure 2.Figure 2
Orbital stability.Definition 4 (asymptotic orbital stability [36–38]).
An orbitγ
(
x
^
) is said to be asymptotically orbitally stable if it is orbitally stable and additionally V may be chosen such that, for all x
∈
V, there exists a constant τ
(
x
) so that γ
t
x
-
γ
t
-
τ
x
x
^
→
0 as t
→
∞, as illustrated in Figure 3.Figure 3
Asymptotic orbital stability.Definition 5 (successor function [37]).
Suppose that the control setM
I
M
P and the phase set N
P
H
A are both lines, as illustrated in Figure 4. The coordinate on the phase set N
P
H
A is defined as follows: the origin is defined as the intersection point between N
P
H
A and the s-axis, denoted as O
′, and the coordinate of any point P on N
P
H
A is defined as the distance between P and O
′. Let P
3 denote the intersection point between the trajectory starting from P
1 and the control set M
I
M
P, and P
4 denote the phase point of P
3 after control measures. Let P
2 denote the intersection point if the trajectory starting from P
1 intersects the phase set N
P
H
A. Then P
4 is called as the successor point of P
1, and the mapping from P
1 to P
4 constructs the successor function. If P
2 exists, the successor function is called type-II; that is, f
S
O
R
I
I
(
P
1
)
=
d
(
P
4
,
O
′
)
-
d
(
P
1
,
O
′
); otherwise, it is called type-I; that is, f
S
O
R
I
(
P
1
)
=
d
(
P
4
,
O
′
)
-
d
(
P
1
,
O
′
). Moreover, the successor function is continuous on N
P
H
A.Figure 4
Illustration of successor function.Remark 6.
If there exists a pointP
(
s
P
,
x
P
)
∈
N
P
H
A such that f
S
O
R
I
(
P
)
=
0 or f
S
O
R
I
I
(
P
)
=
0, then the orbit starting from P forms a periodic orbit.


## 3. Main Results

It is visible thatM
I
M
P
=
{
(
s
,
x
)
∣
s
=
s
C
L
,
x
≥
0
} and(9)
N
P
H
A
=
s
,
x
∣
s
=
s
U
P
≜
1
-
k
w
f
s
f
+
1
-
1
-
k
w
f
s
C
L
,
x
≥
0
.


### 3.1. Qualitative Analysis of the Bioprocess Model

By (4) it is easily obtained that(10)
d
x
=
-
a
+
b
s
d
s
.Integrating both sides of (10) from t
0 to t yields (11)
x
t
-
x
t
0
=
a
s
t
0
+
b
s
2
t
0
2
-
a
s
t
-
b
s
2
t
2
.Now it is assumed that (
s
(
t
0
)
,
x
(
t
0
)
)
∈
N
P
H
A; that is, s
(
t
0
)
=
s
U
P. Let t
1 be the first time for the trajectory to reach the control set M
I
M
P; that is, s
(
t
1
)
=
s
C
L. Then(12)
x
t
1
-
x
t
0
=
a
s
t
0
+
b
s
2
t
0
2
-
a
s
t
1
-
b
s
2
t
1
2
=
a
s
U
P
-
s
C
L
+
b
s
U
P
2
-
s
C
L
2
2
≜
Δ
s
. Denote x
1
=
x
(
t
1
). Due to the effects of the impulsive control, one has(13)
x
t
1
+
=
x
t
1
-
1
-
1
-
θ
k
w
f
x
t
1
. Thus(14)
x
t
1
+
-
x
t
0
=
x
t
1
-
1
-
1
-
θ
k
w
f
x
t
1
-
x
t
0
=
Δ
s
-
1
-
1
-
θ
k
w
f
x
t
1
=
1
-
1
-
1
-
θ
k
w
f
Δ
s
-
1
-
1
-
θ
k
w
f
x
t
0
.Theorem 7.
System (4) admits a unique order-1 periodic solution.Proof.
For anyP
(
s
P
,
x
P
)
∈
N
P
H
A, by Definition 5, one has(15)
f
S
O
R
I
s
P
,
x
P
=
1
-
1
-
1
-
θ
k
w
f
Δ
s
-
1
-
1
-
θ
k
w
f
x
P
. Define(16)
x
¯
1
≜
a
s
U
P
-
s
C
L
+
b
s
U
P
2
-
s
C
L
2
/
2
1
-
1
-
θ
k
w
f
,
x
¯
0
≜
1
-
1
-
1
-
θ
k
w
f
x
¯
1
.Denote P
¯
(
s
U
P
,
x
¯
0
). Then there is f
S
O
R
I
(
P
¯
)
=
0, which implies that the trajectory of (4) starting from P
¯ first intersects the control set M
I
M
P at Q
¯
(
s
C
L
,
x
¯
1
) and then jumps to P
¯ again, which forms an order-1 periodic solution. In addition, f
S
O
R
I
(
P
)
<
0 for x
P
>
x
¯
0 and f
S
O
R
I
(
P
)
>
0 for x
P
<
x
¯
0. This means that f
S
O
R
I
(
P
)
=
0 if and only if x
P
=
x
¯
0; that is, the order-1 periodic solution is unique.Theorem 8.
System (4) does not admit any order-k(k
≥
2) periodic solution.Proof.
Define(17)
f
x
=
x
+
Δ
s
,
g
x
=
x
-
1
-
1
-
θ
k
w
f
x
,
x
≥
0
. Then f
(
x
¯
0
)
=
x
¯
1 and g
(
x
¯
1
)
=
x
¯
0. In addition, there is f
S
O
R
I
(
s
U
P
,
x
)
=
g
(
f
(
x
)
)
-
x. For any point P
(
s
P
,
x
P
)
∈
N
P
H
A, if x
P
>
x
¯
0, then f
(
x
P
)
>
f
(
x
¯
0
)
=
x
¯
1; thus x
P
>
g
(
f
(
x
P
)
)
>
g
(
x
¯
1
)
=
x
¯
0. Similarly, if x
P
<
x
¯
0, then f
(
x
P
)
<
f
(
x
¯
0
)
=
x
¯
1; thus x
P
<
g
(
f
(
x
P
)
)
<
g
(
x
¯
1
)
=
x
¯
0. This implies that the order-k (k
≥
2) periodic solution does not exist.Lets
=
ξ
(
t
)
,
x
=
η
(
t
) be the order-1 periodic solution determined in Theorem 7; that is, ξ
(
t
0
)
=
s
U
P
,
η
(
t
0
)
=
x
¯
0 and ξ
(
t
0
+
k
T
)
=
s
C
L
,
η
(
t
0
+
k
T
)
=
x
¯
1
,
k
≥
0.Theorem 9.
The period of the order-1 periodic solution is(18)
T
=
∫
s
C
L
s
U
P
a
+
b
s
K
ν
Γ
0
-
a
s
-
b
s
2
/
2
ν
+
s
μ
M
s
Γ
0
-
a
s
-
b
s
2
/
2
d
s
,where Γ
0
≜
x
¯
0
+
a
s
U
P
+
b
s
U
P
2
/
2.Proof.
For the order-1 periodic solutions
=
ξ
(
t
)
,
x
=
η
(
t
), by (11), one has(19)
η
t
=
η
t
0
+
a
ξ
t
0
+
b
ξ
2
t
0
2
-
a
ξ
t
-
b
ξ
2
t
2
,
t
0
≤
t
≤
t
0
+
T
.Meanwhile, by (4), one has(20)
d
t
=
-
a
+
b
s
K
ν
η
ν
+
ξ
μ
M
ξ
η
d
ξ
.Substituting (19) into (20) yields (21)
d
t
=
-
a
+
b
ξ
K
ν
Γ
0
-
a
ξ
-
b
ξ
2
/
2
ν
+
ξ
μ
M
ξ
Γ
0
-
a
ξ
-
b
ξ
2
/
2
d
ξ
.Integrating both sides of (21) from t
0 to t
0
+
T yields (22)
T
=
∫
s
C
L
s
U
P
a
+
b
s
K
ν
Γ
0
-
a
s
-
b
s
2
/
2
ν
+
s
μ
M
s
Γ
0
-
a
s
-
b
s
2
/
2
d
s
.Besides the qualitative analysis of model (4) presented above, the stability of the controlled system is also an important issue to be considered. As far as the stability is concerned, most of the works in literature rely on the Analogue of Poincaré Criterion [36]. Recently, a geometric method is used to verify the stability of the order-1 periodic solution for a general semicontinuous dynamic system [13]. For the convenience of readers, a detailed overview is presented.


### 3.2. Stability Criterion Representation

To test the stability of the order-1 periodic solutions
=
ξ
(
t
)
,
x
=
η
(
t
), the behaviour of the orbit of system (5) starting from U
δ
o
(
A
)
∩
N
P
H
A for given small δ
>
0 should be characterized. Let R be a point on the order-1 periodic orbit with rectangular coordinates (
ϕ
(
l
)
,
φ
(
l
)
), where 0
≤
l
≤
L is the arc length between A and R. To begin with, assume R
=
A; that is, l
=
0. For a given point D
∈
U
δ
o
(
A
)
∩
N
P
H
A, here assume that D locates above A, where the rectangular coordinates of A are (
φ
(
l
)
,
ψ
(
l
)
). The trajectory starting from D will intersect the normal of A and the intersection point is denoted by A
k. When t increases, this trajectory intersects the normal of B and the control set M
I
M
P at B
k and C, respectively.Define(23)
ρ
Γ
≜
max
⁡
ρ
D
∣
D
∈
U
δ
o
A
∩
N
P
H
A
. If ρ
Γ
<
1, then the order-1 periodic solution s
=
ξ
(
t
)
,
x
=
η
(
t
) is orbitally asymptotically stable.Theorem 10 (stability criterion [13]).
The order-1T-periodic solution s
=
ξ
(
t
)
,
x
=
η
(
t
) of system (5) is orbitally asymptotically stable if(24)
∫
0
T
∂
F
∂
s
ξ
t
,
η
t
+
∂
G
∂
x
ξ
t
,
η
t
d
t
<
Ξ
,where(25)
Ξ
≜
l
n
F
I
ϕ
s
+
G
I
ϕ
x
F
+
I
1
+
β
x
ϕ
s
-
β
s
ϕ
x
+
G
+
I
1
+
α
s
ϕ
x
-
α
x
ϕ
s
,where F
I
(
G
I
) represents the value of F
(
G
) at C
(
ξ
(
T
)
,
η
(
T
)
)
∈
M
I
M
P, and F
+
I
(
G
+
I
) represents the value of F
(
G
) at A
(
ξ
(
0
)
,
η
(
0
)
)
∈
N
P
H
A.Proof.
As illustrated in Figure5, the slope of the control set is K
M
I
M
P
=
-
ϕ
s
′
/
ϕ
x
′ (in case of ϕ
x
′
=
0, K
M
I
M
P is equal to ∞). Due to the specific forms of the impulsive control, for any R
(
s
R
,
x
R
)
∈
M
I
M
P, if α
s
′, α
x
′, β
s
′, and β
x
′ are constants, one has(26)
s
R
+
x
R
+
=
1
+
α
s
′
α
x
′
β
s
′
1
+
β
x
′
s
R
x
R
+
α
0,0
β
0,0
.
Denote(27)
Φ
=
1
+
α
s
′
α
x
′
β
s
′
1
+
β
x
′
. Then(28)
ϕ
s
′
,
ϕ
x
′
Φ
-
1
Φ
s
R
+
x
R
+
-
ϕ
s
′
,
ϕ
x
′
Φ
-
1
Φ
α
0,0
β
0,0
=
-
ϕ
0,0
, which implies that the phase set N
P
H
A is also a straight line with slope K
N
P
H
A, where(29)
K
N
P
H
A
=
-
ϕ
s
′
1
+
β
x
′
-
ϕ
x
′
β
s
′
ϕ
x
′
1
+
α
s
′
-
ϕ
s
′
α
x
′
=
β
s
′
+
1
+
β
x
′
K
M
I
M
P
1
+
α
s
′
+
α
x
′
K
M
I
M
P
.
On the other hand, by the effect of the impulsive controls, one has(30)
d
A
,
C
¯
=
1
+
α
s
′
+
K
M
I
M
P
α
x
2
+
β
s
′
+
K
M
I
M
P
1
+
β
x
′
2
s
C
-
s
B
2
. Since d
(
B
,
C
)
=
(
s
C
-
s
B
)
2
+
(
x
C
-
x
B
)
2
=
(
1
+
K
M
I
M
P
2
)
(
s
C
-
s
B
)
2, then(31)
d
A
,
C
¯
d
A
,
D
=
1
+
α
s
′
+
K
M
I
M
P
α
x
′
1
+
K
N
P
H
A
2
1
+
K
M
I
M
P
2
d
B
,
C
d
A
,
D
.
Letθ denote the angle between the normal at A and the phase line and let ϑ denote the angle between the normal at B and the impulse line. When D is sufficiently close to A, one has d
(
A
,
A
k
)
≈
d
A
,
D
cos
⁡
θ
,
d
(
B
,
B
k
)
≈
d
B
,
C
cos
⁡
(
ϑ
), which implies that(32)
d
B
,
C
d
A
,
D
=
d
B
,
B
k
d
A
,
A
k
cos
⁡
θ
cos
⁡
ϑ
.
Denoteθ
1 (ϑ
1) as the angle between N
P
H
A (M
I
M
P) and the s-axis and θ
2 (ϑ
2) as the angle between the tangent at A (B) and the s-axis. Then there are t
a
n
(
θ
1
)
=
-
K
N
P
H
A
,
t
a
n
(
θ
2
)
=
x
A
′, t
a
n
(
ϑ
1
)
=
-
K
M
I
M
P, and t
a
n
(
ϑ
2
)
=
x
B
′. Thus(33)
cos
⁡
θ
cos
⁡
ϑ
=
x
A
′
-
K
N
P
H
A
x
B
′
-
K
M
I
M
P
1
+
x
B
′
2
1
+
K
M
I
M
P
2
1
+
x
A
′
2
1
+
K
N
P
H
A
2
.
Summing up (31)–(33) yields (34)
d
A
,
C
¯
d
A
,
D
=
1
+
α
s
′
+
K
M
I
M
P
α
x
′
x
A
′
-
K
N
P
H
A
x
B
′
-
K
M
I
M
P
1
+
x
B
′
2
1
+
x
A
′
2
d
B
,
B
k
d
A
,
A
k
.
To characterize the ratiod
(
B
,
B
k
)
/
d
(
A
,
A
k
), it needs to introduce the curvilinear coordinates (
l
,
n
); the description is as follows: l represents length of the arc starting from A and its increasing direction is consistent with that of time t; n represents the normal length and its positive direction is to the right side when traveling along the periodic orbit.
The relationship between the rectangular coordinates(
s
,
x
) and the curvilinear coordinates (
l
,
n
) on point R is s
=
ϕ
(
l
)
-
n
φ
′
(
l
)
,
x
=
φ
(
l
)
+
n
ϕ
′
(
l
), where(35)
ϕ
′
l
=
d
d
l
ϕ
=
F
0
F
0
2
+
G
0
2
,
φ
′
l
=
d
d
l
φ
=
G
0
F
0
2
+
G
0
2
.
F
0 and G
0 denote the values of F and G at the order-1 periodic orbit x
=
(
ϕ
(
l
)
,
φ
(
l
)
); that is, F
0
=
F
(
ϕ
(
l
)
,
φ
(
l
)
) and G
0
=
G
(
ϕ
(
l
)
,
φ
(
l
)
). It can be inferred that(36)
d
d
l
n
=
G
ϕ
′
-
F
φ
′
F
ϕ
′
+
G
φ
′
-
n
F
ϕ
′
′
+
G
φ
′
′
F
ϕ
′
+
G
φ
′
≜
J
1
l
,
n
-
J
2
l
,
n
=
J
l
,
n
.
Sinceϕ and φ are periodic functions of l, then J
(
l
,
n
) in (36) is also a periodic coefficient nonlinear equation with l as the variable and n as the undetermined function. By (35) it can be obtained that n
=
0 is a solution of (36) corresponding to the order-1 periodic orbit. Since F and G are continuously differentiable, J
(
l
,
n
) also has continuous partial derivatives with respect to n, and (36) can be rewritten as(37)
d
d
l
n
=
J
n
′
l
,
n
n
=
0
·
n
+
o
n
.
In order to calculateJ
n
′
(
l
,
n
)
|
n
=
0, it is noted that (
ϕ
′
(
l
)
)
2
+
(
φ
′
(
l
)
)
2
=
1, which implies that ϕ
′
ϕ
′
′
+
φ
′
φ
′
′
=
0; that is, F
0
ϕ
′
′
+
G
0
φ
′
′
=
0. Thus, (38)
J
2
n
n
=
0
=
F
ϕ
′
′
+
G
φ
′
′
F
ϕ
′
+
G
φ
′
n
=
0
+
n
∂
∂
n
F
ϕ
′
′
+
G
φ
′
′
F
ϕ
′
+
G
φ
′
n
=
0
=
0
. Therefore,(39)
J
n
′
l
,
n
n
=
0
=
J
1
n
′
l
,
n
n
=
0
=
F
0
2
G
x
0
-
F
0
G
0
F
x
0
+
G
s
0
+
G
0
2
F
s
0
F
0
2
+
G
0
2
3
/
2
≜
H
l
,where H
(
l
) represents the curvature of the orthogonal trajectory of (5) at point R. The first-order approximation of (37) is d
n
/
d
l
=
H
(
l
)
n, with the solution(40)
n
l
=
n
0
e
∫
0
l
H
ς
d
ς
n
0
=
n
0
.
Note thatd
l
=
ϕ
˙
2
+
φ
˙
2
d
t
=
(
F
0
)
2
+
(
G
0
)
2
d
t. Then(41)
∫
0
L
H
ζ
d
ζ
=
∫
0
+
t
L
1
F
0
2
+
G
0
2
F
0
2
G
x
0
-
F
0
G
0
F
x
0
+
G
s
0
+
G
0
2
F
s
0
d
t
, which results in the fact that(42)
n
L
n
0
=
e
∫
0
L
H
ζ
d
ζ
=
F
+
I
2
+
G
+
I
2
F
I
2
+
G
I
2
e
∫
0
+
t
L
∂
F
/
∂
s
+
∂
G
/
∂
x
ϕ
t
,
φ
t
d
t
.Combined with (34), this yields the conclusion.Figure 5
Illustration of the disturbance nearby the periodic orbit of model (5).Theorem 11.
The order-1 periodic solutions
=
ξ
(
t
), x
=
η
(
t
) determined in Theorem 7 is orbitally asymptotically stable.Proof.
DenoteA
(
s
U
P
,
x
¯
0
) and B
(
s
C
L
,
x
¯
1
), where s
m
a
x is defined by (9) and x
¯
1, x
¯
0 are defined by (16). In system (4), there are(43)
F
s
,
x
=
-
1
a
+
b
s
μ
s
,
x
x
,
G
s
,
x
=
μ
s
,
x
x
,
α
s
,
x
=
1
-
k
w
f
s
f
-
s
,
β
s
,
x
=
-
w
f
1
-
1
-
θ
k
x
,
ϕ
s
,
x
=
s
-
s
C
L
, where(44)
μ
s
,
x
=
μ
M
s
K
ν
x
ν
+
s
. Then(45)
∂
F
∂
s
=
-
1
a
+
b
s
∂
μ
s
,
x
∂
s
x
+
b
a
+
b
s
2
μ
s
,
x
x
,
∂
G
∂
x
=
μ
s
,
x
+
∂
μ
s
,
x
∂
x
x
,
α
s
′
=
-
1
-
k
w
f
,
α
x
′
=
0
,
β
s
′
=
0
,
β
x
′
=
-
w
f
1
-
1
-
θ
k
,
ϕ
s
′
=
1
,
ϕ
x
′
=
0
. Thus(46)
∂
F
∂
s
+
∂
G
∂
x
=
-
1
a
+
b
s
∂
μ
s
,
x
∂
s
x
+
∂
μ
s
,
x
∂
x
x
+
b
a
+
b
s
2
μ
s
,
x
x
+
μ
s
,
x
. Hence, the left side of (24) is(47)
∫
0
T
∂
F
∂
s
+
∂
G
∂
x
ξ
,
η
d
t
=
∫
0
T
-
1
a
+
b
ξ
∂
μ
ξ
,
η
∂
s
η
+
∂
μ
ξ
,
η
∂
x
η
+
b
η
a
+
b
ξ
2
+
1
μ
ξ
,
η
d
t
=
∫
η
0
η
1
1
η
d
η
-
∫
ξ
0
ξ
1
b
a
+
b
ξ
d
ξ
+
∫
μ
ξ
0
,
η
0
μ
ξ
1
,
η
1
1
μ
d
μ
=
ln
⁡
η
1
η
0
a
+
b
ξ
0
a
+
b
ξ
1
μ
ξ
1
,
η
1
μ
ξ
0
,
η
0
=
ln
⁡
F
ξ
1
,
η
1
F
ξ
0
,
η
0
<
ln
⁡
F
ξ
1
,
η
1
1
-
w
f
1
-
1
-
θ
k
F
ξ
0
,
η
0
=
Ξ
≜
ln
⁡
F
I
ϕ
s
′
+
G
I
ϕ
x
′
F
+
I
1
+
β
x
′
ϕ
s
′
-
β
s
′
ϕ
x
′
+
G
+
I
1
+
α
s
′
ϕ
x
′
-
α
x
′
ϕ
s
′
. Therefore, the order-1 periodic solution s
=
ξ
(
t
)
,
x
=
η
(
t
) is orbitally asymptotically stable.


## 3.1. Qualitative Analysis of the Bioprocess Model

By (4) it is easily obtained that(10)
d
x
=
-
a
+
b
s
d
s
.Integrating both sides of (10) from t
0 to t yields (11)
x
t
-
x
t
0
=
a
s
t
0
+
b
s
2
t
0
2
-
a
s
t
-
b
s
2
t
2
.Now it is assumed that (
s
(
t
0
)
,
x
(
t
0
)
)
∈
N
P
H
A; that is, s
(
t
0
)
=
s
U
P. Let t
1 be the first time for the trajectory to reach the control set M
I
M
P; that is, s
(
t
1
)
=
s
C
L. Then(12)
x
t
1
-
x
t
0
=
a
s
t
0
+
b
s
2
t
0
2
-
a
s
t
1
-
b
s
2
t
1
2
=
a
s
U
P
-
s
C
L
+
b
s
U
P
2
-
s
C
L
2
2
≜
Δ
s
. Denote x
1
=
x
(
t
1
). Due to the effects of the impulsive control, one has(13)
x
t
1
+
=
x
t
1
-
1
-
1
-
θ
k
w
f
x
t
1
. Thus(14)
x
t
1
+
-
x
t
0
=
x
t
1
-
1
-
1
-
θ
k
w
f
x
t
1
-
x
t
0
=
Δ
s
-
1
-
1
-
θ
k
w
f
x
t
1
=
1
-
1
-
1
-
θ
k
w
f
Δ
s
-
1
-
1
-
θ
k
w
f
x
t
0
.Theorem 7.
System (4) admits a unique order-1 periodic solution.Proof.
For anyP
(
s
P
,
x
P
)
∈
N
P
H
A, by Definition 5, one has(15)
f
S
O
R
I
s
P
,
x
P
=
1
-
1
-
1
-
θ
k
w
f
Δ
s
-
1
-
1
-
θ
k
w
f
x
P
. Define(16)
x
¯
1
≜
a
s
U
P
-
s
C
L
+
b
s
U
P
2
-
s
C
L
2
/
2
1
-
1
-
θ
k
w
f
,
x
¯
0
≜
1
-
1
-
1
-
θ
k
w
f
x
¯
1
.Denote P
¯
(
s
U
P
,
x
¯
0
). Then there is f
S
O
R
I
(
P
¯
)
=
0, which implies that the trajectory of (4) starting from P
¯ first intersects the control set M
I
M
P at Q
¯
(
s
C
L
,
x
¯
1
) and then jumps to P
¯ again, which forms an order-1 periodic solution. In addition, f
S
O
R
I
(
P
)
<
0 for x
P
>
x
¯
0 and f
S
O
R
I
(
P
)
>
0 for x
P
<
x
¯
0. This means that f
S
O
R
I
(
P
)
=
0 if and only if x
P
=
x
¯
0; that is, the order-1 periodic solution is unique.Theorem 8.
System (4) does not admit any order-k(k
≥
2) periodic solution.Proof.
Define(17)
f
x
=
x
+
Δ
s
,
g
x
=
x
-
1
-
1
-
θ
k
w
f
x
,
x
≥
0
. Then f
(
x
¯
0
)
=
x
¯
1 and g
(
x
¯
1
)
=
x
¯
0. In addition, there is f
S
O
R
I
(
s
U
P
,
x
)
=
g
(
f
(
x
)
)
-
x. For any point P
(
s
P
,
x
P
)
∈
N
P
H
A, if x
P
>
x
¯
0, then f
(
x
P
)
>
f
(
x
¯
0
)
=
x
¯
1; thus x
P
>
g
(
f
(
x
P
)
)
>
g
(
x
¯
1
)
=
x
¯
0. Similarly, if x
P
<
x
¯
0, then f
(
x
P
)
<
f
(
x
¯
0
)
=
x
¯
1; thus x
P
<
g
(
f
(
x
P
)
)
<
g
(
x
¯
1
)
=
x
¯
0. This implies that the order-k (k
≥
2) periodic solution does not exist.Lets
=
ξ
(
t
)
,
x
=
η
(
t
) be the order-1 periodic solution determined in Theorem 7; that is, ξ
(
t
0
)
=
s
U
P
,
η
(
t
0
)
=
x
¯
0 and ξ
(
t
0
+
k
T
)
=
s
C
L
,
η
(
t
0
+
k
T
)
=
x
¯
1
,
k
≥
0.Theorem 9.
The period of the order-1 periodic solution is(18)
T
=
∫
s
C
L
s
U
P
a
+
b
s
K
ν
Γ
0
-
a
s
-
b
s
2
/
2
ν
+
s
μ
M
s
Γ
0
-
a
s
-
b
s
2
/
2
d
s
,where Γ
0
≜
x
¯
0
+
a
s
U
P
+
b
s
U
P
2
/
2.Proof.
For the order-1 periodic solutions
=
ξ
(
t
)
,
x
=
η
(
t
), by (11), one has(19)
η
t
=
η
t
0
+
a
ξ
t
0
+
b
ξ
2
t
0
2
-
a
ξ
t
-
b
ξ
2
t
2
,
t
0
≤
t
≤
t
0
+
T
.Meanwhile, by (4), one has(20)
d
t
=
-
a
+
b
s
K
ν
η
ν
+
ξ
μ
M
ξ
η
d
ξ
.Substituting (19) into (20) yields (21)
d
t
=
-
a
+
b
ξ
K
ν
Γ
0
-
a
ξ
-
b
ξ
2
/
2
ν
+
ξ
μ
M
ξ
Γ
0
-
a
ξ
-
b
ξ
2
/
2
d
ξ
.Integrating both sides of (21) from t
0 to t
0
+
T yields (22)
T
=
∫
s
C
L
s
U
P
a
+
b
s
K
ν
Γ
0
-
a
s
-
b
s
2
/
2
ν
+
s
μ
M
s
Γ
0
-
a
s
-
b
s
2
/
2
d
s
.Besides the qualitative analysis of model (4) presented above, the stability of the controlled system is also an important issue to be considered. As far as the stability is concerned, most of the works in literature rely on the Analogue of Poincaré Criterion [36]. Recently, a geometric method is used to verify the stability of the order-1 periodic solution for a general semicontinuous dynamic system [13]. For the convenience of readers, a detailed overview is presented.


## 3.2. Stability Criterion Representation

To test the stability of the order-1 periodic solutions
=
ξ
(
t
)
,
x
=
η
(
t
), the behaviour of the orbit of system (5) starting from U
δ
o
(
A
)
∩
N
P
H
A for given small δ
>
0 should be characterized. Let R be a point on the order-1 periodic orbit with rectangular coordinates (
ϕ
(
l
)
,
φ
(
l
)
), where 0
≤
l
≤
L is the arc length between A and R. To begin with, assume R
=
A; that is, l
=
0. For a given point D
∈
U
δ
o
(
A
)
∩
N
P
H
A, here assume that D locates above A, where the rectangular coordinates of A are (
φ
(
l
)
,
ψ
(
l
)
). The trajectory starting from D will intersect the normal of A and the intersection point is denoted by A
k. When t increases, this trajectory intersects the normal of B and the control set M
I
M
P at B
k and C, respectively.Define(23)
ρ
Γ
≜
max
⁡
ρ
D
∣
D
∈
U
δ
o
A
∩
N
P
H
A
. If ρ
Γ
<
1, then the order-1 periodic solution s
=
ξ
(
t
)
,
x
=
η
(
t
) is orbitally asymptotically stable.Theorem 10 (stability criterion [13]).
The order-1T-periodic solution s
=
ξ
(
t
)
,
x
=
η
(
t
) of system (5) is orbitally asymptotically stable if(24)
∫
0
T
∂
F
∂
s
ξ
t
,
η
t
+
∂
G
∂
x
ξ
t
,
η
t
d
t
<
Ξ
,where(25)
Ξ
≜
l
n
F
I
ϕ
s
+
G
I
ϕ
x
F
+
I
1
+
β
x
ϕ
s
-
β
s
ϕ
x
+
G
+
I
1
+
α
s
ϕ
x
-
α
x
ϕ
s
,where F
I
(
G
I
) represents the value of F
(
G
) at C
(
ξ
(
T
)
,
η
(
T
)
)
∈
M
I
M
P, and F
+
I
(
G
+
I
) represents the value of F
(
G
) at A
(
ξ
(
0
)
,
η
(
0
)
)
∈
N
P
H
A.Proof.
As illustrated in Figure5, the slope of the control set is K
M
I
M
P
=
-
ϕ
s
′
/
ϕ
x
′ (in case of ϕ
x
′
=
0, K
M
I
M
P is equal to ∞). Due to the specific forms of the impulsive control, for any R
(
s
R
,
x
R
)
∈
M
I
M
P, if α
s
′, α
x
′, β
s
′, and β
x
′ are constants, one has(26)
s
R
+
x
R
+
=
1
+
α
s
′
α
x
′
β
s
′
1
+
β
x
′
s
R
x
R
+
α
0,0
β
0,0
.
Denote(27)
Φ
=
1
+
α
s
′
α
x
′
β
s
′
1
+
β
x
′
. Then(28)
ϕ
s
′
,
ϕ
x
′
Φ
-
1
Φ
s
R
+
x
R
+
-
ϕ
s
′
,
ϕ
x
′
Φ
-
1
Φ
α
0,0
β
0,0
=
-
ϕ
0,0
, which implies that the phase set N
P
H
A is also a straight line with slope K
N
P
H
A, where(29)
K
N
P
H
A
=
-
ϕ
s
′
1
+
β
x
′
-
ϕ
x
′
β
s
′
ϕ
x
′
1
+
α
s
′
-
ϕ
s
′
α
x
′
=
β
s
′
+
1
+
β
x
′
K
M
I
M
P
1
+
α
s
′
+
α
x
′
K
M
I
M
P
.
On the other hand, by the effect of the impulsive controls, one has(30)
d
A
,
C
¯
=
1
+
α
s
′
+
K
M
I
M
P
α
x
2
+
β
s
′
+
K
M
I
M
P
1
+
β
x
′
2
s
C
-
s
B
2
. Since d
(
B
,
C
)
=
(
s
C
-
s
B
)
2
+
(
x
C
-
x
B
)
2
=
(
1
+
K
M
I
M
P
2
)
(
s
C
-
s
B
)
2, then(31)
d
A
,
C
¯
d
A
,
D
=
1
+
α
s
′
+
K
M
I
M
P
α
x
′
1
+
K
N
P
H
A
2
1
+
K
M
I
M
P
2
d
B
,
C
d
A
,
D
.
Letθ denote the angle between the normal at A and the phase line and let ϑ denote the angle between the normal at B and the impulse line. When D is sufficiently close to A, one has d
(
A
,
A
k
)
≈
d
A
,
D
cos
⁡
θ
,
d
(
B
,
B
k
)
≈
d
B
,
C
cos
⁡
(
ϑ
), which implies that(32)
d
B
,
C
d
A
,
D
=
d
B
,
B
k
d
A
,
A
k
cos
⁡
θ
cos
⁡
ϑ
.
Denoteθ
1 (ϑ
1) as the angle between N
P
H
A (M
I
M
P) and the s-axis and θ
2 (ϑ
2) as the angle between the tangent at A (B) and the s-axis. Then there are t
a
n
(
θ
1
)
=
-
K
N
P
H
A
,
t
a
n
(
θ
2
)
=
x
A
′, t
a
n
(
ϑ
1
)
=
-
K
M
I
M
P, and t
a
n
(
ϑ
2
)
=
x
B
′. Thus(33)
cos
⁡
θ
cos
⁡
ϑ
=
x
A
′
-
K
N
P
H
A
x
B
′
-
K
M
I
M
P
1
+
x
B
′
2
1
+
K
M
I
M
P
2
1
+
x
A
′
2
1
+
K
N
P
H
A
2
.
Summing up (31)–(33) yields (34)
d
A
,
C
¯
d
A
,
D
=
1
+
α
s
′
+
K
M
I
M
P
α
x
′
x
A
′
-
K
N
P
H
A
x
B
′
-
K
M
I
M
P
1
+
x
B
′
2
1
+
x
A
′
2
d
B
,
B
k
d
A
,
A
k
.
To characterize the ratiod
(
B
,
B
k
)
/
d
(
A
,
A
k
), it needs to introduce the curvilinear coordinates (
l
,
n
); the description is as follows: l represents length of the arc starting from A and its increasing direction is consistent with that of time t; n represents the normal length and its positive direction is to the right side when traveling along the periodic orbit.
The relationship between the rectangular coordinates(
s
,
x
) and the curvilinear coordinates (
l
,
n
) on point R is s
=
ϕ
(
l
)
-
n
φ
′
(
l
)
,
x
=
φ
(
l
)
+
n
ϕ
′
(
l
), where(35)
ϕ
′
l
=
d
d
l
ϕ
=
F
0
F
0
2
+
G
0
2
,
φ
′
l
=
d
d
l
φ
=
G
0
F
0
2
+
G
0
2
.
F
0 and G
0 denote the values of F and G at the order-1 periodic orbit x
=
(
ϕ
(
l
)
,
φ
(
l
)
); that is, F
0
=
F
(
ϕ
(
l
)
,
φ
(
l
)
) and G
0
=
G
(
ϕ
(
l
)
,
φ
(
l
)
). It can be inferred that(36)
d
d
l
n
=
G
ϕ
′
-
F
φ
′
F
ϕ
′
+
G
φ
′
-
n
F
ϕ
′
′
+
G
φ
′
′
F
ϕ
′
+
G
φ
′
≜
J
1
l
,
n
-
J
2
l
,
n
=
J
l
,
n
.
Sinceϕ and φ are periodic functions of l, then J
(
l
,
n
) in (36) is also a periodic coefficient nonlinear equation with l as the variable and n as the undetermined function. By (35) it can be obtained that n
=
0 is a solution of (36) corresponding to the order-1 periodic orbit. Since F and G are continuously differentiable, J
(
l
,
n
) also has continuous partial derivatives with respect to n, and (36) can be rewritten as(37)
d
d
l
n
=
J
n
′
l
,
n
n
=
0
·
n
+
o
n
.
In order to calculateJ
n
′
(
l
,
n
)
|
n
=
0, it is noted that (
ϕ
′
(
l
)
)
2
+
(
φ
′
(
l
)
)
2
=
1, which implies that ϕ
′
ϕ
′
′
+
φ
′
φ
′
′
=
0; that is, F
0
ϕ
′
′
+
G
0
φ
′
′
=
0. Thus, (38)
J
2
n
n
=
0
=
F
ϕ
′
′
+
G
φ
′
′
F
ϕ
′
+
G
φ
′
n
=
0
+
n
∂
∂
n
F
ϕ
′
′
+
G
φ
′
′
F
ϕ
′
+
G
φ
′
n
=
0
=
0
. Therefore,(39)
J
n
′
l
,
n
n
=
0
=
J
1
n
′
l
,
n
n
=
0
=
F
0
2
G
x
0
-
F
0
G
0
F
x
0
+
G
s
0
+
G
0
2
F
s
0
F
0
2
+
G
0
2
3
/
2
≜
H
l
,where H
(
l
) represents the curvature of the orthogonal trajectory of (5) at point R. The first-order approximation of (37) is d
n
/
d
l
=
H
(
l
)
n, with the solution(40)
n
l
=
n
0
e
∫
0
l
H
ς
d
ς
n
0
=
n
0
.
Note thatd
l
=
ϕ
˙
2
+
φ
˙
2
d
t
=
(
F
0
)
2
+
(
G
0
)
2
d
t. Then(41)
∫
0
L
H
ζ
d
ζ
=
∫
0
+
t
L
1
F
0
2
+
G
0
2
F
0
2
G
x
0
-
F
0
G
0
F
x
0
+
G
s
0
+
G
0
2
F
s
0
d
t
, which results in the fact that(42)
n
L
n
0
=
e
∫
0
L
H
ζ
d
ζ
=
F
+
I
2
+
G
+
I
2
F
I
2
+
G
I
2
e
∫
0
+
t
L
∂
F
/
∂
s
+
∂
G
/
∂
x
ϕ
t
,
φ
t
d
t
.Combined with (34), this yields the conclusion.Figure 5
Illustration of the disturbance nearby the periodic orbit of model (5).Theorem 11.
The order-1 periodic solutions
=
ξ
(
t
), x
=
η
(
t
) determined in Theorem 7 is orbitally asymptotically stable.Proof.
DenoteA
(
s
U
P
,
x
¯
0
) and B
(
s
C
L
,
x
¯
1
), where s
m
a
x is defined by (9) and x
¯
1, x
¯
0 are defined by (16). In system (4), there are(43)
F
s
,
x
=
-
1
a
+
b
s
μ
s
,
x
x
,
G
s
,
x
=
μ
s
,
x
x
,
α
s
,
x
=
1
-
k
w
f
s
f
-
s
,
β
s
,
x
=
-
w
f
1
-
1
-
θ
k
x
,
ϕ
s
,
x
=
s
-
s
C
L
, where(44)
μ
s
,
x
=
μ
M
s
K
ν
x
ν
+
s
. Then(45)
∂
F
∂
s
=
-
1
a
+
b
s
∂
μ
s
,
x
∂
s
x
+
b
a
+
b
s
2
μ
s
,
x
x
,
∂
G
∂
x
=
μ
s
,
x
+
∂
μ
s
,
x
∂
x
x
,
α
s
′
=
-
1
-
k
w
f
,
α
x
′
=
0
,
β
s
′
=
0
,
β
x
′
=
-
w
f
1
-
1
-
θ
k
,
ϕ
s
′
=
1
,
ϕ
x
′
=
0
. Thus(46)
∂
F
∂
s
+
∂
G
∂
x
=
-
1
a
+
b
s
∂
μ
s
,
x
∂
s
x
+
∂
μ
s
,
x
∂
x
x
+
b
a
+
b
s
2
μ
s
,
x
x
+
μ
s
,
x
. Hence, the left side of (24) is(47)
∫
0
T
∂
F
∂
s
+
∂
G
∂
x
ξ
,
η
d
t
=
∫
0
T
-
1
a
+
b
ξ
∂
μ
ξ
,
η
∂
s
η
+
∂
μ
ξ
,
η
∂
x
η
+
b
η
a
+
b
ξ
2
+
1
μ
ξ
,
η
d
t
=
∫
η
0
η
1
1
η
d
η
-
∫
ξ
0
ξ
1
b
a
+
b
ξ
d
ξ
+
∫
μ
ξ
0
,
η
0
μ
ξ
1
,
η
1
1
μ
d
μ
=
ln
⁡
η
1
η
0
a
+
b
ξ
0
a
+
b
ξ
1
μ
ξ
1
,
η
1
μ
ξ
0
,
η
0
=
ln
⁡
F
ξ
1
,
η
1
F
ξ
0
,
η
0
<
ln
⁡
F
ξ
1
,
η
1
1
-
w
f
1
-
1
-
θ
k
F
ξ
0
,
η
0
=
Ξ
≜
ln
⁡
F
I
ϕ
s
′
+
G
I
ϕ
x
′
F
+
I
1
+
β
x
′
ϕ
s
′
-
β
s
′
ϕ
x
′
+
G
+
I
1
+
α
s
′
ϕ
x
′
-
α
x
′
ϕ
s
′
. Therefore, the order-1 periodic solution s
=
ξ
(
t
)
,
x
=
η
(
t
) is orbitally asymptotically stable.


## 4. Verifications and Biomass Productivity Optimisation

The microbial cultivation process for pulse dosage supply of substrates and removal of bioreactor medium has been analyzed theoretically. In order to complement the theoretical results, simulations of system (4)’s dynamics are presented firstly. The model parameters used in simulations are as follows: μ
M
=
0.5
[
1
/
h
], ν
=
0.5, K
ν
=
1
[
g
/
L
], s
f
=
16
[
g
/
L
], a
=
0.5, b
=
0.025, θ
=
0.8, k
=
0.4, and w
f
=
15
%
=
0.15. The low substrate control level is a decision parameter, which should be set according to an actual demand. Firstly, in order to verify the theoretical results, s
C
L is set to 10
%
s
f
=
1.6
[
g
/
L
]. Next, an optimal s
C
L
∗ is determined by a biomass productivity minimization problem.


### 4.1. Verification of the Main Results

The change of substrate concentration (s), biomass concentration (
x
), and phase diagram (
s
,
x
) starting from (
s
0
,
x
0
)
=
(
12,0.5
) are shown in Figure 6. It can be observed that the trajectories tend to be periodic, and this phenomenon is not dependent on the model parameters. Furthermore, from Figure 6(b), it can be concluded that the biomass concentration is controlled below a certain level through substrate regulation. Figure 7 illustrates the order-1 periodic solution with period T
=
0.6
[
h
]
=
36
[
m
i
n
]. The phase diagrams for different initial biomass concentrations on N
P
H
A are illustrated in Figure 8, from which the stability of the order-1 periodic solution can be observed.Figure 6
Verification of Theorem7. The change of substrate concentration (s), biomass concentration (
x
), and phase diagrams with (
s
0
,
x
0
)
=
(
6,2
).Figure 7
Verification of Theorem7. The change of substrate concentration (s), biomass concentration (
x
), and phase diagrams with (
s
0
,
x
0
)
=
(
2.9,4.5
).Figure 8
Verification of Theorem7. The phase diagrams for different initial biomass concentrations on N
P
H
A.


### 4.2. Biomass Productivity Optimisation

Besides analyzing the dynamics of the proposed bioprocess model, another important aspect is to find the optimal substrate control level such that the biomass productivity is maximized. In the proposed bioprocess, the biomass productivity(
B
P
x
[
g
/
L
h
]
) is formulated as follows:(48)
B
P
x
w
f
,
s
C
L
=
θ
w
f
x
O
U
T
T
w
f
,
s
C
L
, where x
O
U
T
=
x
¯
1 is determined by (16) and T
(
w
f
,
s
C
L
) is determined by (18). For the given model parameters, the dependence of B
P
x on w
f and the substrate control level s
C
L is shown in Figure 9. The maximum biomass productivity occurs for w
f
∗
=
0.1 and s
C
L
∗
=
6
[
g
/
L
]. Note that, for the maximum biomass productivity, the bioreactor medium removal part is also kept at the minimum value, which increases the frequency of the operation. For bigger w
f, the dependence of B
P
x on the substrate control level s
C
L is shown in Figure 10, and, for different substrate control level s
C
L, the dependence of B
P
x on w
f is illustrated in Figure 11.Figure 9
The dependence ofB
P
x on the control parameters w
f and s
C
L.Figure 10
The dependence ofB
P
x on the substrate control level s
C
L for different w
f.Figure 11
The dependence ofB
P
x on w
f for different substrate control level s
C
L.The model parameters also have a certain impact on the biomass productivity. For givens
C
L
=
6
[
g
/
L
], w
f
=
0.1, the biomass productivity achieves its maximum at θ
∗
=
1 and k
∗
=
0, as illustrated in Figure 12. This is easily understood since in case of full filter efficiency the recycling operation is unnecessary. The recycling takes effect only for low filter efficiency. The dependence of B
P
x on the model parameters k for different θ is given in Figure 13, and the dependence of B
P
x on the model parameter θ for different k is given in Figure 14. It can be concluded that the optimal substrate recycling proportion k is dependent on θ; that is, when θ is small, for example, θ
≤
0.2, k should be set equal to k
m
a
x
=
90
%; when θ
=
0.4, k should be set equal to 40
%; while, for bigger θ, for example, θ
≥
0.5, k should be set equal to zero, which means that the substrate recycling operation is unnecessary. But it should be pointed that the recycling operation is indeed beneficial in saving the substrate loss.Figure 12
The dependence ofB
P
x on the model parameters k and θ for the given s
C
L
=
6
[
g
/
L
] and w
f
=
0.1.Figure 13
The dependence ofB
P
x on the model parameters k for the given s
C
L
=
6
[
g
/
L
], w
f
=
0.1.Figure 14
The dependence ofB
P
x on the model parameters θ for the given s
C
L
=
6
[
g
/
L
], w
f
=
0.1.


## 4.1. Verification of the Main Results

The change of substrate concentration (s), biomass concentration (
x
), and phase diagram (
s
,
x
) starting from (
s
0
,
x
0
)
=
(
12,0.5
) are shown in Figure 6. It can be observed that the trajectories tend to be periodic, and this phenomenon is not dependent on the model parameters. Furthermore, from Figure 6(b), it can be concluded that the biomass concentration is controlled below a certain level through substrate regulation. Figure 7 illustrates the order-1 periodic solution with period T
=
0.6
[
h
]
=
36
[
m
i
n
]. The phase diagrams for different initial biomass concentrations on N
P
H
A are illustrated in Figure 8, from which the stability of the order-1 periodic solution can be observed.Figure 6
Verification of Theorem7. The change of substrate concentration (s), biomass concentration (
x
), and phase diagrams with (
s
0
,
x
0
)
=
(
6,2
).Figure 7
Verification of Theorem7. The change of substrate concentration (s), biomass concentration (
x
), and phase diagrams with (
s
0
,
x
0
)
=
(
2.9,4.5
).Figure 8
Verification of Theorem7. The phase diagrams for different initial biomass concentrations on N
P
H
A.


## 4.2. Biomass Productivity Optimisation

Besides analyzing the dynamics of the proposed bioprocess model, another important aspect is to find the optimal substrate control level such that the biomass productivity is maximized. In the proposed bioprocess, the biomass productivity(
B
P
x
[
g
/
L
h
]
) is formulated as follows:(48)
B
P
x
w
f
,
s
C
L
=
θ
w
f
x
O
U
T
T
w
f
,
s
C
L
, where x
O
U
T
=
x
¯
1 is determined by (16) and T
(
w
f
,
s
C
L
) is determined by (18). For the given model parameters, the dependence of B
P
x on w
f and the substrate control level s
C
L is shown in Figure 9. The maximum biomass productivity occurs for w
f
∗
=
0.1 and s
C
L
∗
=
6
[
g
/
L
]. Note that, for the maximum biomass productivity, the bioreactor medium removal part is also kept at the minimum value, which increases the frequency of the operation. For bigger w
f, the dependence of B
P
x on the substrate control level s
C
L is shown in Figure 10, and, for different substrate control level s
C
L, the dependence of B
P
x on w
f is illustrated in Figure 11.Figure 9
The dependence ofB
P
x on the control parameters w
f and s
C
L.Figure 10
The dependence ofB
P
x on the substrate control level s
C
L for different w
f.Figure 11
The dependence ofB
P
x on w
f for different substrate control level s
C
L.The model parameters also have a certain impact on the biomass productivity. For givens
C
L
=
6
[
g
/
L
], w
f
=
0.1, the biomass productivity achieves its maximum at θ
∗
=
1 and k
∗
=
0, as illustrated in Figure 12. This is easily understood since in case of full filter efficiency the recycling operation is unnecessary. The recycling takes effect only for low filter efficiency. The dependence of B
P
x on the model parameters k for different θ is given in Figure 13, and the dependence of B
P
x on the model parameter θ for different k is given in Figure 14. It can be concluded that the optimal substrate recycling proportion k is dependent on θ; that is, when θ is small, for example, θ
≤
0.2, k should be set equal to k
m
a
x
=
90
%; when θ
=
0.4, k should be set equal to 40
%; while, for bigger θ, for example, θ
≥
0.5, k should be set equal to zero, which means that the substrate recycling operation is unnecessary. But it should be pointed that the recycling operation is indeed beneficial in saving the substrate loss.Figure 12
The dependence ofB
P
x on the model parameters k and θ for the given s
C
L
=
6
[
g
/
L
] and w
f
=
0.1.Figure 13
The dependence ofB
P
x on the model parameters k for the given s
C
L
=
6
[
g
/
L
], w
f
=
0.1.Figure 14
The dependence ofB
P
x on the model parameters θ for the given s
C
L
=
6
[
g
/
L
], w
f
=
0.1.


## 5. Conclusions

In the paper the dynamic behaviour of a microbial cultivation process model with variable biomass yield and substrate regulation was analyzed. For the proposed bioprocess model, it was shown that (i) the existence of the order-1 periodic solution does not depend on the biomass yield and the kinetic model (Theorem7), but the bioprocess behaviour (i.e., the order-1 periodic solution) is dependent on the biomass yield and the kinetic model (see (11)–(18)); (ii) the system is not chaotic due to nonexistence of the order-k periodic solution (Theorem 8) and orbit asymptotic stability of the order-1 periodic solution (Theorem 11); (iii) the stability of the order-1 periodic solution is independent from the biomass yield and the kinetic model. The analytical results offer the possibility of establishing general and more systematic operation and control strategies based on the counteraction of the mechanisms underlying the adverse effects of bioreactor dynamics. In addition, microorganisms in the considered bioprocess always kept a suitable growth rate through substrate regulation. The biomass concentration could also be controlled to a certain level, which should ensure an appropriate dissolved oxygen concentration. Moreover, through biomass productivity optimisation, the optimal substrate control level and bioreactor medium removal proportion were obtained. The study showed how the biomass productivity depended on the bioprocess parameters (i.e., the substrate recycling proportion k and the biomass filter efficiency θ).



---

*Source: 3685941-2016-05-16.xml*