# Theoretical Study and Optimization of the Biochemical Reaction Process by Means of Feedback Control Strategy

Authors: Kaibiao Sun; Andrzej Kasperski; Yuan Tian
Journal: Journal of Chemistry
 (2013)

---

## Abstract

The aim of this work is to present a theoretical analysis and optimization of a biochemical reaction process by means of feedback control strategy. To begin with, a mathematical model of the biochemical reaction process with feedback control is formulated. Then, based on the formulated model, the analysis of system's dynamics is presented. The optimization of the bioprocess is carried out, in order to achieve maximal biomass productivity. It is shown that during the optimization, the bioprocess with impulse effects loses the possibility of synchronization and strives for a simple continuous bioprocess.
The analytical results presented in the work are validated by numerical simulations for the Tessier kinetics model.


---

## Body

## 1. Introduction

Chemical and bioprocess engineering play an important role in the production of many chemical products. In particular, bioreactor engineering as a branch of chemical engineering and biotechnology is an active area of research on bioprocesses, including among others development, control, and commercialization of new technology [1]. The reaching of optimal results and the obtainment of maximal profits require modern control strategies based on mathematical models or artificial intelligence methods [2]. There are many advantages with the use of mathematical models, and one of which is the possibility of testing process stability [3]. According to different reactions and differential control technologies, many dynamic biochemical models in a chemostat have been established [4–12]. However, there are many factors which affect the growth and reproduction of microorganisms in the bioprocess. The key variable is the biomass concentration, because, among other things, it provides information concerning biomass productivity. For that reason, almost all mathematical growth models contain the biomass as an important variable [13]. Moreover, for aerobic microbes, the dissolved oxygen concentration (DOC) is an important growth factor. During the growth of microorganisms, the dissolved oxygen concentration depends on several factors, among them the biomass concentration, the concentration and type of substrate used, and the bioprocess conditions. Because a low level of dissolved oxygen concentration decreases the biomass yield and the specific growth rate, it is necessary to monitor and control the DOC level so that it stays within the appropriate range [14].Since many biological phenomena such as bursting rhythm models in, for example, medicine, biology, pharmacokinetics, and frequency modulated systems exhibit impulsive effects [15], impulsive differential equations, which appear as a natural description of observed evolution phenomena of several real world problems, have been introduced in different kind of biological systems, for example, in population dynamics [16–19] and in biochemical process [20, 21]. It should be pointed out that in these works, the authors were concerned about the fixed time impulse effects, which has the rationality in describing the biological phenomena. While in some cases, using the state-dependent impulse effects to describe the biological phenomena is more appropriate. As far as the state-dependent impulse effect is concerned, Tang and Chen [22] introduced a Lotka-Voterra model, which is constructed according to the practices of IPM. By using analytical method, it is shown that there exists an orbitally asymptotically stable periodic solution with a maximum value no larger than the given economic threshold. Further, the complete expression of period of the periodic solution is given. More researches on the applications of the state-dependent impulse effects can be found in [23–36]. In these researches, the authors analyzed the proposed system's dynamic behavior (e.g., the existence and stability of period-1 solution and the existence of period-2 solution) by applying the Poincaré principle and Poincaré-Bendixson theorem of the impulsive differential equation. Recently, in microbioprocess engineering, the dynamic properties of the kinetic models with impulse effects characterized by the universal microorganism growth rate and two different kinds of biomass yield are also analyzed theoretically by Sun et al. [37–39].The objective of this work is to illuminate theoretical and practical aspects of the nonlinear analysis of a universal mathematical model of the biochemical reaction process. The paper is organized as follows. In Section2, a universal mathematical model of the biochemical reaction process with any characteristics of growth kinetics is formulated under impulse effects. In Section 3, a qualitative analysis of the proposed model is presented. In Section 4, simulations with Tessier's kinetics are presented in order to verify the theoretical results and discuss the biological essence. Moreover, in Section 4, in order to optimize the biochemical reaction process, the conception of objective function is presented, and next the bioprocess optimization is put forward. Finally, in Section 5, we offer our conclusions.


## 2. Mathematical Model and Preliminaries

The general model of continuously culturing microorganism in a chemostat is given by the following form of differential equations [40]:
(1)dSdt=D(Sin-S)-1Yx/Sμ(S)x,dxdt=(μ(S)-D)x,S(0)=S0,x(0)=x0,
where x=x(t) denotes the biomass concentration (g/L), and S=S(t) the substrate concentration (g/L) in the bioreactor medium at time t, x0, and S0 denotes the initial biomass concentration and substrate concentration, respectively; D is the dilution rate (h-1); Sin is the influent substrate concentration (g/L); Yx/S is the biomass yield (g/g) defined as the ratio of the biomass x produced to the amount of substrate S consumed; Yx/S<1 for biological constraints; the function μ(S) describes the biochemical kinetics, which are characterized by the cell concentration (x), depending on one limiting substrate with concentration S.Crooke et al. [6] pointed out that the biochemical kinetics expression plays an important role in the intrinsic oscillation mechanisms. So many different assumptions for the kinetics models are given in the literature, for example, Monod-type kinetics [3], that is, μ(S)=μmaxS/(KS+S), where μmax is the maximum specific growth rate and KS is the substrate saturation constant; Tessier-type kinetics [41], that is, μ(S)=μmax(1-e-S/KS); Moser-type kinetics [41], that is, μ(S)=μmax(1+KSS-ϑ)-1, where ϑ is a positive constant. It can be easily seen that these kinetic models satisfy the following properties:  (i)
(regularity)μ: ℝ+→ℝ+ is continuously differential and μ(0)=0,(ii)
(monotonicity)μ is monotonically increasing, that is, dμ/dS≥0 for all S∈ℝ+; μ′(S¯)=0 for some S¯∈ℝ+ means that μ′(S)=0 for S≥S¯,(iii)
(convexity)μ(θS1+(1-θ)S2)≥θμ(S1)+(1-θ)μ(S2) for all θ∈[0,1] and S1,S2∈ℝ+.Thus, in the study that the kinetic models only satisfy the regularity, monotonicity, and convexity are considered.According to the Herbert's and Pirt's models [41], if the substrate concentration is big enough (i.e., S≫KS), the biomass yield is constant. Because in continuous bioprocesses the optimal biomass productivity is obtained for bigger substrate concentration, in this work the constant biomass yield is assumed, that is, Yx/S=Yx/Smax=λ<1. In addition, according to the design ideas of the bioreactor, the biomass concentration in the chemostat should be controlled to a set level xset, where 0<xset≤xCRI, and xCRI is the critical level of biomass concentration in the bioreactor medium. When the biomass concentration reaches the set level, part of the medium containing biomass and substrate is discharged from the bioreactor, followed by the input of the “clear” medium (i.e., the medium without substrate). Therefore, system (1) can be modified by introducing the impulsive state feedback control as follows:
(2)dSdt=D(Sin-S)-1λμ(S)xdxdt=(μ(S)-D)x}x<xset,ΔS=-WfSΔx=-Wfx}x=xset,S(0+)=S0,x(0+)=x0,
where Wf is the part of “clear” medium inputted into the bioreactor in each biomass oscillation cycle, which satisfies 0≤Wf≤Wfmax<1, where Wfmax is the maximum part of “clear” medium inputted into the bioreactor in each biomass oscillation cycle.Before presenting the main results, we recall the following definitions and lemmas first [31, 42, 43].Definition 1.
(
ξ
(
t
)
,
η
(
t
)
) is said to be period-1 solution if in a minimum cycle time there is one impulse effect. Similarly, (ξ(t),η(t)) is said to be period-2 solution if in a minimum cycle time there are two impulse effects.Definition 2.
Γ is said to be orbitally stable, if for any ɛ>0, there exists δ>0, with the proviso that every solution (ξ(t),η(t)) of system (2) whose distance from Γ is less than δ at t=t0, will remain within a distance less than ɛ from Γ for all t≥t0. Such a Γ is said to be orbitally asymptotically stable if, in addition, the distance of (ξ(t),η(t)) from Γ tends to zero as t→∞. Moreover, if there exist positive constants α,β, and a real constant C such that ρ((ξ(t),η(t),Γ)<αe-βt for t>C, then Γ is said to be orbitally asymptotically stable and enjoys the property of asymptotic phase.Definition 3.
LambertW function is defined to be a multivalued inverse of the function z→zez satisfying
(3)LambertW(z)exp(LambertW(z))=z.In the followings, we simply denoteW≜LambertW. Then, it follows from (3) that
(4)W'(z)=W(z)z(1+W(z)),
when z≠0 and z≠-1/e. First of all, the function zexp(z) has the positive derivative (z+1)exp(z) if z>-1. Define the inverse function of zexp(z) restricted on the interval [-1,∞) to be W(0,z)≜W(z). Similarly, we define the inverse function of zexp(z) restricted on the interval (-∞,-1] to be W(-1,z)≜W(z). For this study, both W(0,z) and W(-1,z) will be employed only for z∈[-exp(-1),0) due to our practical problem (see Figure 1). For more details of the concepts and properties of the Lambert W function, see Corless et al. [43].Figure 1
The property of LambertW function.Lemma 4.
Assume that there exists a bounded closed regionG=ABCDA, as shown in Figure 2. The boundaries AD^ and BC^ are the no cut-arcs of model (5), and on AD and BC, the direction of the orientation field determined by model (5) is pointing to the interior of G. In addition, there is no singularity in the interior of G, and the boundary. One of the boundary CD is the impulse set of the model (5); the corresponding phase set satisfies I(CD)⊂AB; AB is also the no cut-arcs of model (5), and on AB, the direction of the orientation field determined by model (5) is pointing to the interior of G. Then, there must exist at least one period-1 solution of model (5) in region G.Figure 2
Illustration of Bendixson field.Lemma 5 (analogue of Poincaré’ criterion).
TheT-periodic solution S=ξ(t), x=η(t) of system
(5)dSdt=P(S,x),dxdt=R(S,x),ifϕ(S,x)≠0,ΔS=α(S,x),Δx=β(S,x),ifϕ(S,x)=0,
is orbitally asymptotically stable and enjoying the property of asymptotic phase if the multiplier μ2 satisfies the condition |μ2|<1 and unstable if |μ2|>1, where
(6)μ2=∏k=1q‍Δkexp(∫0T‍[∂P∂S(ξ(t),η(t))+∂R∂x(ξ(t),η(t))]dt),Δk=P+((∂β/∂X)(∂ϕ/∂S)-(∂β/∂S)(∂ϕ/∂x)+∂ϕ/∂S)P(∂ϕ/∂S)+R(∂ϕ/∂x)+R+((∂α/∂S)(∂ϕ/∂x)-(∂α/∂x)(∂ϕ/∂S)+∂ϕ/∂x)P(∂ϕ/∂S)+R(∂ϕ/∂x),P+=P(ξ(τk+),η(τk+)), R+=R(ξ(τk+),η(τk+)), and P,R,∂α/∂S,∂α/∂x,∂β/∂S,∂β/∂x,∂ϕ/∂S,∂ϕ/∂x are calculated at the point (ξ(τk),η(τk)).


## 3. Properties Analysis of the Bioprocess

We focus our discussion on the caseSin>μ-1(D), xset<λ(Sin-μ-1(D)), and (S0,x0)∈ΩT={(S,x)∣0<S≤Sin,0<x<xset}.


### 3.1. Existence of Period-1 Solution

The linex=xset intersects the isoclinal line dx/dt=0 at the point A, intersects the curve dS/dt=0 at the point M(SM,xset), and intersects the line S=Sin at the point B, where SM satisfies that
(7)λDSM+μ(SM)xset=λDSin.
The line x=(1-Wf)xset intersects the line dx/dt=0 at the point E and intersects the line S=Sin at the point H. Denote C=C(SC,(1-Wf)xset) and D=D(SD,(1-Wf)xset), where SC=(1-Wf)μ-1(D) and SD=(1-Wf)Sin. The impulsive set Mimp⊆AB¯={(S,x)∣μ-1(D)≤S≤Sin,x=xset}, and the phase set Npha=I(Mimp)⊆CD={(S,x)∣SC≤S≤SD,x=(1-Wf)xset}. Let the point K′ be the primary image of E, and denote the trajectories of the system starting from the points C and E that intersect the segment AB¯ at C' and E' respectively. Denote the trajectory of the system between C and C' by CC'^, the trajectory of the system between E and E' by EE'^.Theorem 6.
If one of the conditions holds (i)SD≤μ-1(D)<Sin and xset<λ(Sin-μ-1(D)); (ii) If μ-1(D)<SD, SC'≤SK' and xset<λ(Sin-μ-1(D)); (iii) If μ-1(D)<SD, SE'≥SK' and xset<λ(Sin-μ-1(D)), then system (2) has a period-1 solution with (S0,x0)∈ΩT.Proof.
Firstly, it is obvious that all trajectories of system (2) starting from the region ΩT must interest with the segment AB¯ and then jump to the segment CD¯ if μ-1(D)<Sin and xset<λ(Sin-μ-1(D)). Next, we construct the closed region G⊂ΩT such that all solutions of system (2) starting from ΩT enter into G and retain there.
Case (i) (SD≤μ-1(D)). It is obvious that the straight line AC is described by x=xsetS/μ-1(D), and BD is described by x=xsetS/Sin, and both lines pass through the point O(0,0). The derivative of the straight line L(S,x) passing through the points A(μ-1(D),xset) and D((1-Wf)Sin,(1-Wf)xset) along the trajectories of system (2) is
(8)dLdt|(2)=dxdt|(2)-Wfxsetμ-1(D)-(1-Wf)SindSdt|(2)<0,
in terms of dx/dt|(2)<0, dS/dt|(2)>0 for the segment AD¯. From the qualitative characteristic of system (2), we know that dL/dt|(2)>0 for AD¯ and dx/dt|(2)<0 for CD¯. Besides, dS/dt|(2)<0 for S=Sin. On the other hand, it follows from system (2) that x(t)=0 and S(t)=Sin-(Sin-S0)e-Dt for t∈(0,+∞). Especially, when S0=SC, FG→ is the semitrivial solution of system (2). Therefore, we have found a closed region G, the boundary of which consists of AD¯, DC¯, CF¯,  FG¯,  GB¯, and BA¯; see Figure 3(a).
Case (ii) (SD>μ-1(D) and SC′≤SK′). Similar to the discussion of Case (i), the boundary of G can be constructed by EC¯, CC'^, C'A¯, and AE¯; see Figure 3(b).
Case (iii) (SD>μ-1(D) and SE′≥SK′). In this case, the boundary of G can be constructed by ED¯, DH¯, HB¯, BE'¯, and EE'^; see Figure 3(c).
In addition, there is no singularity in the regionG. Therefore, it follows from the qualitative characteristics of the system that all the trajectories satisfying the conditions of theorem enter the closed region G and retain there. Then it follows from Lemma 4 that system (2) has a period-1 solution.Illustration of system (2) when μ-1(D)<Sin and xset<λ(Sin-μ-1(D)): (a) SD≤μ-1(D); (b) SD>μ-1(D) and SC'≤SK'; (c) SD>μ-1(D) and SE'≥SK'.
(a)
(b)
(c)


### 3.2. Position of Period-1 Solution

Suppose the period of the period-1 solution isT. Let S(t)=ξ(t),  x(t)=η(t) be such a T-periodic solution. Then ζ(t)=λξ(t)+η(t) is also T-periodic. Denote ξ0=ξ(0+)=ξ(T+),   ξ1=ξ(T),  η0=η(0+)=η(T+), and η1=η(T). Then from the T-periodicity of the periodic solution and the third and fourth equations of system (2), we have ξ0=(1-Wf)ξ1,  η0=(1-Wf)η1, and η1=xset. Without loss of generality, for t∈[0,T], ζ(t) satisfies the relation
(9)dζdt=D(λSin-ζ).
Then, we have
(10)ζ(t)=λSin+(ζ(0+)-λSin)exp(-Dt).
In particular, for t=T, we have ζ(T)=λSin(1-exp(-DT))+ζ(0+)exp(-DT).By (10), if ζ(0+)≥λSin, then we have λSin≤ζ(t)<λSin+(ζ(0+)-λSin)=ζ(0+) for all t∈(0,T], else we have λSin≥ζ(t)>ζ(0) for all t∈(0,T].LetSR>μ-1(D) be the root of the following equation:
(11)SR-∫μ-1(D)SR‍Dμ(ϑ)dϑ=μ-1(D)-xsetλ(Wf+ln(1-Wf)).
Then, we have the following results.Theorem 7.
Suppose thatSin>μ-1(D) and xset≤λ(Sin-μ-1(D)). Then the initial values of period-1 solution (ξ(t),η(t)) of system (2) satisfies (1-Wf)min{SR,SM}≤ξ(0+)<(1-Wf)(Sin-xset/λ) and η(0+)=(1-Wf)xset, where SM and SR are determined by (7) and (11), respectively.Proof.
Sinceζ(0+)=λξ(0+)+η(0+)=λ(1-Wf)ξ(T)+(1-Wf)η(T)=(1-Wf)(λξ(T)+η(T))=(1-Wf)ζ(T). Then, ζ(T)=λSin(1-exp(-DT))+ζ(0+)exp(-DT)=λSin(1-exp(-DT))+(1-Wf)ζ(T)exp(-DT), which implies that
(12)ζ(T)=1-exp(-DT)1-(1-Wf)exp(-DT)λSin<λSin,
that is, ξ1=ξ(T)<Sin-xset/λ.
Next, we consider the following comparison system of system (2) without impulsive feedback control:
(13)dSdt=1λμ(S)(xset-x),dxdt=[μ(S)-D]x.
It is easy to verify that system (13) has one positive equilibrium (μ-1(D),xset), which is a center point. Then all of its solutions are closed trajectories which satisfy
(14)V1(S,x)≜x-xsetlnx+λS-λ∫SMS‍Dμ(ϑ)dϑ-xset+xsetlnxset-λSM=C,
where C is an arbitrary constant. The derivative of V1(S,x) along the trajectories of system (2) is
(15)dV1dt|(2)=∂V1∂xdxdt+∂V1∂SdSdt=(1-xsetx)[μ(S)-D]x+λ(1-Dμ(S))[D(Sin-S)-1λμ(S)x]=λ[μ(S)-D]μ(S)[D(Sin-S)-1λμ(S)xset]>0,
for μ-1(D)<S<SM. This implies that the trajectories of system (2) intersecting with the trajectories of system (13) pass through the trajectories of system (13) from the left to the right. In addition, the trajectory of system (13) passing through M(SM,xset) (i.e., C=0 in (14)) intersects the line S=μ-1(D) at two points L1(μ-1(D),xL1) and L2(μ-1(D),xL2), where xLi(i=1,2) satisfy that (16)lnxLi-xLixset=Λ1≜[λμ-1(D)-λ∫SMμ-1(D)‍(D/μ(ϑ))dϑ-xset+xsetlnxset-λSM]xset,that is,
(17)-xLixsete-xLi/xset=-eΛ1xset.
Since Λ1=ln(xLi)-xLi/xset<lnxset-1, then we have -eΛ1/xset>-1/e. Thus, we have
(18)xL1=-xsetW(-1,-eΛ1xset),xL2=-xsetW(0,-eΛ1xset),
and xL2<xset<xL1, where W denotes Lambert W function defined in Definition 3.
IfxL2≥(1-Wf)xset, then the closed region G constructed in Theorem 6 can be reduced to G′, the boundary of which consists of ML2^, L2E¯, EC¯, CF¯, FG¯, GB¯, and BM¯. Thus, the period-1 solution satisfies that ξ1=ξ(T)≥xM>μ-1(D). Else, we have xL2<(1-Wf)xset. Then, the trajectory of system (13) passing through E intersects the segment AM at the point R(SR,xset). In this case, the closed region G constructed in Theorem 6 can be reduced to G', the boundary of which consists of RL2^, L2E¯, EC¯, CF¯, FG¯, GB¯, and BR¯. Thus, the period-1 solution satisfies that ξ1=ξ(T)≥SR>μ-1(D).Corollary 8.
Suppose thatμ-1(D)<Sin≤μ-1(D)/(1-Wf) and xset≤λ(Sin-μ-1(D)). Then, the initial values of the period-1 solution (ξ(t),η(t)) of system (2) satisfies that η(0+)=(1-Wf)xset and (1-Wf)min{SR,SM}≤ξ(0+)<μ-1(D).Corollary 9.
Suppose thatSin>μ-1(D)/(1-Wf) and max{λ[Sin-μ-1(D)/(1-Wf)],0}≤xset≤λ(SF-μ-1(D)). Then the initial values of the period-1 solution (ξ(t),η(t)) of system (2) satisfies that (1-Wf)min{SR,SM}≤ξ(0+)<μ-1(D) and η(0+)=(1-Wf)xset.Corollary 10.
Suppose thatSin>μ-1(D), xset≤λ(SF-μ-1(D)) and Wf≥1-μ-1(D)/(Sin-xset/λ). Then, the initial values of the period-1 solution (ξ(t),η(t)) of system (2) satisfies that (1-Wf)min{SR,SM}≤ξ(0+)<μ-1(D) and η(0+)=(1-Wf)xset.


### 3.3. Stability of Period-1 Solution

LetS*=μ-1((2-Wf)D/(1+(1-Wf)2)). Denote
(19)x*=λD(Sin-S*)μ(S*),x¯=-x*W(0,-eΛ2x*),
where Λ2≜[λμ-1(D)-λ∫S*μ-1(D)‍(D/μ(ϑ))dϑ-xset+x*lnxset-λS*]/x*, and W denotes Lambert W function defined in Definition 3. Then, we have the following result.Theorem 11.
Suppose thatSin>μ-1(D), xset≤min{x*,x¯/(1-Wf)}. Then, the period-1 solution (ξ(t),η(t)) of system (2) is orbitally asymptotically stable and enjoys the property of asymptotic phase.Proof.
LetS*=μ-1((2-Wf)D/(1+(1-Wf)2)). Denote
(20)x*=λD(Sin-S*)μ(S*),x¯=-x*W(0,-eΛ2x*),
where Λ2≜[λμ-1(D)-λ∫S*μ-1(D)‍(D/μ(ϑ))dϑ-xset+x*lnxset-λS*]/x*, and W denotes Lambert W function defined in Definition 3. According to Lemma 5, we calculate the multiplier of system (2) in variations corresponding to the T-periodic solution (ξ(t),η(t)). Denote P0(ξ0,η0), P1(ξ1,η1), where ξ0=(1-Wf)ξ1, η0=(1-Wf)xset, μ-1(D)<ξ1<Sin-xset/λ, and η1=xset. In system (2), since P(S,x)=D(Sin-S)-μ(S)x/λ, R(S,x)=(μ(S)-D)x, α(S,x)=-WfS, β(S,x)=-Wfx, and ϕ(S,x)=x-xset, then
(21)∂P∂S=-D-1λμ′(S)x,∂R∂x=μ(S)-D,∂α∂S=-Wf,∂α∂x=0,∂β∂S=0,∂β∂x=-Wf,∂ϕ∂S=0,∂ϕ∂x=1.
Therefore,
(22)Λ1=P+((∂β/∂x)(∂ϕ/∂S)-(∂β/∂S)(∂ϕ/∂x)+∂ϕ/∂S)P(∂ϕ/∂S)+R(∂ϕ/∂x)+R+((∂α/∂S)(∂ϕ/∂x)-(∂α/∂x)(∂ϕ/∂S)+∂ϕ/∂x)P(∂ϕ/∂S)+R(∂ϕ/∂x)=R+(1-Wf)R=(1-Wf)2μ(ξ0)-Dμ(ξ1)-D,μ2=Λ1exp(∫0T‍[∂P∂S(ξ(t),η(t))+∂R∂x(ξ(t),η(t))]dt)=Λ1exp(∫0T‍[-D-1λμ'(S)x]dt+∫0T‍[μ(S)-D]dt)=Λ1exp(∫η0η1‍1xdx-∫0T‍[D+1λμ′(S)x]dt)=(1-Wf)μ(ξ0)-Dμ(ξ1)-Dexp(-∫0T‍[D+1λμ′(S)x]dt).
Note that ξ1=ξ(T)≥min{SM,SR}>μ-1(D), that is, μ(ξ1)>D. If ξ0=(1-Wf)ξ1≥μ-1(D), that is, the periodic solution lies in the region, where dx/dt>0, then we have 0≤μ2<1, and the periodic solution (ξ(t),η(t)) is orbitally asymptotically stable and enjoying the property of asymptotic phase. But for ξ0<μ-1(D), that is, ξ1<μ-1(D)/(1-Wf), we have μ(ξ0)-D<0. Since the expression of ξ1 is hard to be given, we cannot determine whether the periodic solution is stable or not. However, we can give some sufficient conditions to guarantee -1<μ2<0 when ξ1<μ-1(D)/(1-Wf), that is, the periodic solution is orbitally asymptotically stable and enjoying the property of asymptotic phase.
Since0<Wf<1, then we have 1<(2-Wf)/(1+(1-Wf)2)<1/(1-Wf), that is, D<(2-Wf)D/(1+(1-Wf)2)=μ(S*)<D/(1-Wf), which implies that μ-1(D)<S*<μ-1(D)/(1-Wf). If ξ1≥S*, then we have μ(ξ1)+(1-Wf)μ(ξ0)≥(1+(1-Wf)2)μ(ξ1)≥(2-Wf)D, which is equivalent to (1-Wf)(D-μ(ξ0))≤μ(ξ1)-D or (1-Wf)(D-μ(ξ0))/(μ(ξ1)-D)≤1. Thus, we get |μ2|<1, that is, the periodic solution is orbitally asymptotically stable and enjoying the property of asymptotic phase. In the following, we will give the conditions under which ξ1≥S*.
From the expression of the isoclinal linedS/dt=0, it follows that for S=S*, we have x=λD(Sin-S*)/μ(S*):=x*.
Consider the following comparison system of system (2) without impulsive feedback control
(23)dSdt=1λμ(S)(x*-x),dxdt=[μ(S)-D]x.
It is easy to verify that system (23) has one positive equilibrium (μ-1(D),x*), which is a center point. Then, all of its solutions are closed trajectories which satisfy
(24)V2(S,x)≜x-x*lnx+λS-λ∫S*S‍Dμ(ϑ)dϑ-xset+x*lnxset-λS*=C,
where C is an arbitrary constant. The derivative of V2(S,x) along the trajectories of system (2) is
(25)dV2dt|(2)=∂V2∂xdxdt+∂V2∂SdSdt=(1-x*x)[μ(S)-D]x+λ(1-Dμ(S))[D(Sin-S)-1λμ(S)x]=λ[μ(S)-D]μ(S)[D(Sin-S)-1λμ(S)x*]>0,
for μ-1(D)<S<S*. This implies that the trajectories of system (2) intersecting with the trajectories of system (23) pass through the trajectories of system (23) from the left to the right. In addition, the trajectory of system (23) passing through N(S*,xset) (i.e., C=0 in (24)) intersects the line S=μ-1(D) at two points K1(μ-1(D),xK1) and K2(μ-1(D),xK2), where xKi(i=1,2) satisfy that(26)lnxKi-xKix*=Λ2≜[λμ-1(D)-λ∫S*μ-1(D)‍(D/μ(ϑ))dϑ-xset+x*lnxset-λS*]x*,that is,
(27)-xKix*e-xKi/x*=-eΛ2x*.
Since Λ2=ln(xKi)-xKi/x*<lnx*-1, then we have -eΛ2/x*>-1/e. Thus, by Definition 3, we have
(28)xK1=-x*W(-1,-eΛ2x*),xK2=-x*W(0,-eΛ2x*),
and xK2<x*<xK1.
Next, we will prove thatξ1≥S* when the conditions of theorem are satisfied. Let N be the intersection point of the line x=x* and dS/dt=0, K1, and let K2 be the intersection points of the line S=μ-1(D) and the trajectory of system (23) passing through the point M. If xset≤min{x*,x¯/(1-Wf)}, then we have (1-Wf)xset≤x¯. Since dV2/dt|(2)>0, dS/dt|K1E¯>0, dx/dt|EC¯<0, dS/dt|CF¯>0, dS/dt|GB¯>0, and FG→ is semi-trivial solution of system (2). Therefore, we get a closed region G, the boundary of which consists of NK1^, K1E¯, EC¯, CF¯, FG¯, GB¯, and BN¯. In addition, there is no singularity in the region G. Therefore, it follows from the qualitative characteristics of the system that all the trajectories satisfying the conditions of theorem enter the closed region G and retain there. Then, it follows from Theorem 6 that system (2) has a period-1 solution in G. Furthermore, we have S*≤ξ1≤Sin. If ξ1≥μ-1(D), then ξ0≥μ-1(D), in this case, we have 0≤μ2<1. Else S*≤ξ1<μ-1(D)/(1-Wf), in this case, we have -1<μ2<0. Therefore, if the conditions of the theorem are satisfied, then the period-1 solution (ξ(t),η(t)) of system (2) is orbitally asymptotically stable and enjoying the property of asymptotic phase.


### 3.4. Discussion on the Period-k (k≥2) Solution

For the predator-prey model concerning IPM strategies, Tang and Cheke [23] proved that there is no periodic solution with order larger than or equal to three, except for one special case, by using the properties of the Lambert W function and Poincaré map. Moreover, They showed that the existence of an order two periodic solution implies the existence of an order one periodic solution. Next, the existence of the period-k (k≥3) solution is discussed.Theorem 12.
Suppose thatSin>μ-1(D)/(1-Wf) and λ(Sin-μ-1(D)/(1-Wf))≤xset<λ(Sin-μ-1(D)). Then, system (2) does not have period-k (k≥3) solution.Proof.
Let(S,x) be an arbitrary solution of system (2). Denote the first intersection point of the trajectory to the impulsive set M by (S1,xset) (here without loss of generality, we assume that S1>ξ1), and the corresponding consecutive points are (S2,xset), (S3,xset), …, respectively. Consequently, under the effect of impulsive function I, the corresponding points after pulse are (S1+,(1-Wf)xset),(S2+,(1-Wf)xset),(S3+,(1-Wf)xset),…. By the qualitative analysis of system (2), we know that x˙>0 for S>μ-1(D) and x˙<0 for S<μ-1(D). If there exists the jth impulse effect such that Sk+≥μ-1(D) for k≥j, then the sequence will be monotone, that is, Sj+≥Sj+1+≥⋯≥μ-1(D) or μ-1(D)≤Sj+≤Sj+1+≤⋯. In this case, system (2) has no period-k(k≥2) solution. Let l1 be the line passing the point G(Sin,0) with the slope -λ:
(29)l1:x+λS-λSin=0.
Since dl1/dt|(2)=dx/dt+λdS/dt=D(λSin-λS-x)=0, then the trajectory starting from the point on the segment CD either does not intersect l1 or retain on l1 before reaching the segment AB. Therefore, when xset≥λ(Sin-μ-1(D)/(1-Wf)), the period-1 solution determined in Theorem 6 satisfies that (1-Wf)min{SR,SM}≤ξ0<(1-Wf)(Sin-xset/λ)≤μ-1(D), which also implies that all Sk+≤μ-1(D).
The tr
ajectory of the solution(S,x) jumps to the point (S1+,(1-Wf)xset) from the point (S1,xset). By the dynamics of system (2), we have the following two sequences according to (S3,xset) lies on the right or left of (S1,xset): (a)  ⋯≤S4≤S2<ξ1<S1≤S3≤⋯; (b)  S2≤S4≤⋯≤ξ1≤⋯≤S3≤S1. Since the subsequences {S2k-1} and {S2k}   (k=1,2,…) are monotone and bounded, then there exists a limit for them, respectively. Let S_=limk→∞S2k and S¯=limk→∞S2k-1. When sequence (a) occurs, it can be shown that S_<ξ1<S¯. Therefore, system (2) has a period-2 solution. When sequence (b) occurs, we have S_≤ξ1≤S¯. If S_<ξ1<S¯, then system (2) has a period-2 solution. Else we have S_=ξ1=S¯, in this case, the trajectory tends to the period-1 solution (ξ,η). But it can be obtained that there is no order k solution (k≥3) in system (2), thus the system is not chaos by Li-Yorke chaos Theorem [44].Remark 13.
From the proof of Theorem12, it can be seen that there is a possibility to exist a period-2 solution. Since it is difficult to obtain the exact expression of the solution to (2), then it will be challenging to give the sufficient conditions under which system (2) has a period-2 solution or not.


## 3.1. Existence of Period-1 Solution

The linex=xset intersects the isoclinal line dx/dt=0 at the point A, intersects the curve dS/dt=0 at the point M(SM,xset), and intersects the line S=Sin at the point B, where SM satisfies that
(7)λDSM+μ(SM)xset=λDSin.
The line x=(1-Wf)xset intersects the line dx/dt=0 at the point E and intersects the line S=Sin at the point H. Denote C=C(SC,(1-Wf)xset) and D=D(SD,(1-Wf)xset), where SC=(1-Wf)μ-1(D) and SD=(1-Wf)Sin. The impulsive set Mimp⊆AB¯={(S,x)∣μ-1(D)≤S≤Sin,x=xset}, and the phase set Npha=I(Mimp)⊆CD={(S,x)∣SC≤S≤SD,x=(1-Wf)xset}. Let the point K′ be the primary image of E, and denote the trajectories of the system starting from the points C and E that intersect the segment AB¯ at C' and E' respectively. Denote the trajectory of the system between C and C' by CC'^, the trajectory of the system between E and E' by EE'^.Theorem 6.
If one of the conditions holds (i)SD≤μ-1(D)<Sin and xset<λ(Sin-μ-1(D)); (ii) If μ-1(D)<SD, SC'≤SK' and xset<λ(Sin-μ-1(D)); (iii) If μ-1(D)<SD, SE'≥SK' and xset<λ(Sin-μ-1(D)), then system (2) has a period-1 solution with (S0,x0)∈ΩT.Proof.
Firstly, it is obvious that all trajectories of system (2) starting from the region ΩT must interest with the segment AB¯ and then jump to the segment CD¯ if μ-1(D)<Sin and xset<λ(Sin-μ-1(D)). Next, we construct the closed region G⊂ΩT such that all solutions of system (2) starting from ΩT enter into G and retain there.
Case (i) (SD≤μ-1(D)). It is obvious that the straight line AC is described by x=xsetS/μ-1(D), and BD is described by x=xsetS/Sin, and both lines pass through the point O(0,0). The derivative of the straight line L(S,x) passing through the points A(μ-1(D),xset) and D((1-Wf)Sin,(1-Wf)xset) along the trajectories of system (2) is
(8)dLdt|(2)=dxdt|(2)-Wfxsetμ-1(D)-(1-Wf)SindSdt|(2)<0,
in terms of dx/dt|(2)<0, dS/dt|(2)>0 for the segment AD¯. From the qualitative characteristic of system (2), we know that dL/dt|(2)>0 for AD¯ and dx/dt|(2)<0 for CD¯. Besides, dS/dt|(2)<0 for S=Sin. On the other hand, it follows from system (2) that x(t)=0 and S(t)=Sin-(Sin-S0)e-Dt for t∈(0,+∞). Especially, when S0=SC, FG→ is the semitrivial solution of system (2). Therefore, we have found a closed region G, the boundary of which consists of AD¯, DC¯, CF¯,  FG¯,  GB¯, and BA¯; see Figure 3(a).
Case (ii) (SD>μ-1(D) and SC′≤SK′). Similar to the discussion of Case (i), the boundary of G can be constructed by EC¯, CC'^, C'A¯, and AE¯; see Figure 3(b).
Case (iii) (SD>μ-1(D) and SE′≥SK′). In this case, the boundary of G can be constructed by ED¯, DH¯, HB¯, BE'¯, and EE'^; see Figure 3(c).
In addition, there is no singularity in the regionG. Therefore, it follows from the qualitative characteristics of the system that all the trajectories satisfying the conditions of theorem enter the closed region G and retain there. Then it follows from Lemma 4 that system (2) has a period-1 solution.Illustration of system (2) when μ-1(D)<Sin and xset<λ(Sin-μ-1(D)): (a) SD≤μ-1(D); (b) SD>μ-1(D) and SC'≤SK'; (c) SD>μ-1(D) and SE'≥SK'.
(a)
(b)
(c)


## 3.2. Position of Period-1 Solution

Suppose the period of the period-1 solution isT. Let S(t)=ξ(t),  x(t)=η(t) be such a T-periodic solution. Then ζ(t)=λξ(t)+η(t) is also T-periodic. Denote ξ0=ξ(0+)=ξ(T+),   ξ1=ξ(T),  η0=η(0+)=η(T+), and η1=η(T). Then from the T-periodicity of the periodic solution and the third and fourth equations of system (2), we have ξ0=(1-Wf)ξ1,  η0=(1-Wf)η1, and η1=xset. Without loss of generality, for t∈[0,T], ζ(t) satisfies the relation
(9)dζdt=D(λSin-ζ).
Then, we have
(10)ζ(t)=λSin+(ζ(0+)-λSin)exp(-Dt).
In particular, for t=T, we have ζ(T)=λSin(1-exp(-DT))+ζ(0+)exp(-DT).By (10), if ζ(0+)≥λSin, then we have λSin≤ζ(t)<λSin+(ζ(0+)-λSin)=ζ(0+) for all t∈(0,T], else we have λSin≥ζ(t)>ζ(0) for all t∈(0,T].LetSR>μ-1(D) be the root of the following equation:
(11)SR-∫μ-1(D)SR‍Dμ(ϑ)dϑ=μ-1(D)-xsetλ(Wf+ln(1-Wf)).
Then, we have the following results.Theorem 7.
Suppose thatSin>μ-1(D) and xset≤λ(Sin-μ-1(D)). Then the initial values of period-1 solution (ξ(t),η(t)) of system (2) satisfies (1-Wf)min{SR,SM}≤ξ(0+)<(1-Wf)(Sin-xset/λ) and η(0+)=(1-Wf)xset, where SM and SR are determined by (7) and (11), respectively.Proof.
Sinceζ(0+)=λξ(0+)+η(0+)=λ(1-Wf)ξ(T)+(1-Wf)η(T)=(1-Wf)(λξ(T)+η(T))=(1-Wf)ζ(T). Then, ζ(T)=λSin(1-exp(-DT))+ζ(0+)exp(-DT)=λSin(1-exp(-DT))+(1-Wf)ζ(T)exp(-DT), which implies that
(12)ζ(T)=1-exp(-DT)1-(1-Wf)exp(-DT)λSin<λSin,
that is, ξ1=ξ(T)<Sin-xset/λ.
Next, we consider the following comparison system of system (2) without impulsive feedback control:
(13)dSdt=1λμ(S)(xset-x),dxdt=[μ(S)-D]x.
It is easy to verify that system (13) has one positive equilibrium (μ-1(D),xset), which is a center point. Then all of its solutions are closed trajectories which satisfy
(14)V1(S,x)≜x-xsetlnx+λS-λ∫SMS‍Dμ(ϑ)dϑ-xset+xsetlnxset-λSM=C,
where C is an arbitrary constant. The derivative of V1(S,x) along the trajectories of system (2) is
(15)dV1dt|(2)=∂V1∂xdxdt+∂V1∂SdSdt=(1-xsetx)[μ(S)-D]x+λ(1-Dμ(S))[D(Sin-S)-1λμ(S)x]=λ[μ(S)-D]μ(S)[D(Sin-S)-1λμ(S)xset]>0,
for μ-1(D)<S<SM. This implies that the trajectories of system (2) intersecting with the trajectories of system (13) pass through the trajectories of system (13) from the left to the right. In addition, the trajectory of system (13) passing through M(SM,xset) (i.e., C=0 in (14)) intersects the line S=μ-1(D) at two points L1(μ-1(D),xL1) and L2(μ-1(D),xL2), where xLi(i=1,2) satisfy that (16)lnxLi-xLixset=Λ1≜[λμ-1(D)-λ∫SMμ-1(D)‍(D/μ(ϑ))dϑ-xset+xsetlnxset-λSM]xset,that is,
(17)-xLixsete-xLi/xset=-eΛ1xset.
Since Λ1=ln(xLi)-xLi/xset<lnxset-1, then we have -eΛ1/xset>-1/e. Thus, we have
(18)xL1=-xsetW(-1,-eΛ1xset),xL2=-xsetW(0,-eΛ1xset),
and xL2<xset<xL1, where W denotes Lambert W function defined in Definition 3.
IfxL2≥(1-Wf)xset, then the closed region G constructed in Theorem 6 can be reduced to G′, the boundary of which consists of ML2^, L2E¯, EC¯, CF¯, FG¯, GB¯, and BM¯. Thus, the period-1 solution satisfies that ξ1=ξ(T)≥xM>μ-1(D). Else, we have xL2<(1-Wf)xset. Then, the trajectory of system (13) passing through E intersects the segment AM at the point R(SR,xset). In this case, the closed region G constructed in Theorem 6 can be reduced to G', the boundary of which consists of RL2^, L2E¯, EC¯, CF¯, FG¯, GB¯, and BR¯. Thus, the period-1 solution satisfies that ξ1=ξ(T)≥SR>μ-1(D).Corollary 8.
Suppose thatμ-1(D)<Sin≤μ-1(D)/(1-Wf) and xset≤λ(Sin-μ-1(D)). Then, the initial values of the period-1 solution (ξ(t),η(t)) of system (2) satisfies that η(0+)=(1-Wf)xset and (1-Wf)min{SR,SM}≤ξ(0+)<μ-1(D).Corollary 9.
Suppose thatSin>μ-1(D)/(1-Wf) and max{λ[Sin-μ-1(D)/(1-Wf)],0}≤xset≤λ(SF-μ-1(D)). Then the initial values of the period-1 solution (ξ(t),η(t)) of system (2) satisfies that (1-Wf)min{SR,SM}≤ξ(0+)<μ-1(D) and η(0+)=(1-Wf)xset.Corollary 10.
Suppose thatSin>μ-1(D), xset≤λ(SF-μ-1(D)) and Wf≥1-μ-1(D)/(Sin-xset/λ). Then, the initial values of the period-1 solution (ξ(t),η(t)) of system (2) satisfies that (1-Wf)min{SR,SM}≤ξ(0+)<μ-1(D) and η(0+)=(1-Wf)xset.


## 3.3. Stability of Period-1 Solution

LetS*=μ-1((2-Wf)D/(1+(1-Wf)2)). Denote
(19)x*=λD(Sin-S*)μ(S*),x¯=-x*W(0,-eΛ2x*),
where Λ2≜[λμ-1(D)-λ∫S*μ-1(D)‍(D/μ(ϑ))dϑ-xset+x*lnxset-λS*]/x*, and W denotes Lambert W function defined in Definition 3. Then, we have the following result.Theorem 11.
Suppose thatSin>μ-1(D), xset≤min{x*,x¯/(1-Wf)}. Then, the period-1 solution (ξ(t),η(t)) of system (2) is orbitally asymptotically stable and enjoys the property of asymptotic phase.Proof.
LetS*=μ-1((2-Wf)D/(1+(1-Wf)2)). Denote
(20)x*=λD(Sin-S*)μ(S*),x¯=-x*W(0,-eΛ2x*),
where Λ2≜[λμ-1(D)-λ∫S*μ-1(D)‍(D/μ(ϑ))dϑ-xset+x*lnxset-λS*]/x*, and W denotes Lambert W function defined in Definition 3. According to Lemma 5, we calculate the multiplier of system (2) in variations corresponding to the T-periodic solution (ξ(t),η(t)). Denote P0(ξ0,η0), P1(ξ1,η1), where ξ0=(1-Wf)ξ1, η0=(1-Wf)xset, μ-1(D)<ξ1<Sin-xset/λ, and η1=xset. In system (2), since P(S,x)=D(Sin-S)-μ(S)x/λ, R(S,x)=(μ(S)-D)x, α(S,x)=-WfS, β(S,x)=-Wfx, and ϕ(S,x)=x-xset, then
(21)∂P∂S=-D-1λμ′(S)x,∂R∂x=μ(S)-D,∂α∂S=-Wf,∂α∂x=0,∂β∂S=0,∂β∂x=-Wf,∂ϕ∂S=0,∂ϕ∂x=1.
Therefore,
(22)Λ1=P+((∂β/∂x)(∂ϕ/∂S)-(∂β/∂S)(∂ϕ/∂x)+∂ϕ/∂S)P(∂ϕ/∂S)+R(∂ϕ/∂x)+R+((∂α/∂S)(∂ϕ/∂x)-(∂α/∂x)(∂ϕ/∂S)+∂ϕ/∂x)P(∂ϕ/∂S)+R(∂ϕ/∂x)=R+(1-Wf)R=(1-Wf)2μ(ξ0)-Dμ(ξ1)-D,μ2=Λ1exp(∫0T‍[∂P∂S(ξ(t),η(t))+∂R∂x(ξ(t),η(t))]dt)=Λ1exp(∫0T‍[-D-1λμ'(S)x]dt+∫0T‍[μ(S)-D]dt)=Λ1exp(∫η0η1‍1xdx-∫0T‍[D+1λμ′(S)x]dt)=(1-Wf)μ(ξ0)-Dμ(ξ1)-Dexp(-∫0T‍[D+1λμ′(S)x]dt).
Note that ξ1=ξ(T)≥min{SM,SR}>μ-1(D), that is, μ(ξ1)>D. If ξ0=(1-Wf)ξ1≥μ-1(D), that is, the periodic solution lies in the region, where dx/dt>0, then we have 0≤μ2<1, and the periodic solution (ξ(t),η(t)) is orbitally asymptotically stable and enjoying the property of asymptotic phase. But for ξ0<μ-1(D), that is, ξ1<μ-1(D)/(1-Wf), we have μ(ξ0)-D<0. Since the expression of ξ1 is hard to be given, we cannot determine whether the periodic solution is stable or not. However, we can give some sufficient conditions to guarantee -1<μ2<0 when ξ1<μ-1(D)/(1-Wf), that is, the periodic solution is orbitally asymptotically stable and enjoying the property of asymptotic phase.
Since0<Wf<1, then we have 1<(2-Wf)/(1+(1-Wf)2)<1/(1-Wf), that is, D<(2-Wf)D/(1+(1-Wf)2)=μ(S*)<D/(1-Wf), which implies that μ-1(D)<S*<μ-1(D)/(1-Wf). If ξ1≥S*, then we have μ(ξ1)+(1-Wf)μ(ξ0)≥(1+(1-Wf)2)μ(ξ1)≥(2-Wf)D, which is equivalent to (1-Wf)(D-μ(ξ0))≤μ(ξ1)-D or (1-Wf)(D-μ(ξ0))/(μ(ξ1)-D)≤1. Thus, we get |μ2|<1, that is, the periodic solution is orbitally asymptotically stable and enjoying the property of asymptotic phase. In the following, we will give the conditions under which ξ1≥S*.
From the expression of the isoclinal linedS/dt=0, it follows that for S=S*, we have x=λD(Sin-S*)/μ(S*):=x*.
Consider the following comparison system of system (2) without impulsive feedback control
(23)dSdt=1λμ(S)(x*-x),dxdt=[μ(S)-D]x.
It is easy to verify that system (23) has one positive equilibrium (μ-1(D),x*), which is a center point. Then, all of its solutions are closed trajectories which satisfy
(24)V2(S,x)≜x-x*lnx+λS-λ∫S*S‍Dμ(ϑ)dϑ-xset+x*lnxset-λS*=C,
where C is an arbitrary constant. The derivative of V2(S,x) along the trajectories of system (2) is
(25)dV2dt|(2)=∂V2∂xdxdt+∂V2∂SdSdt=(1-x*x)[μ(S)-D]x+λ(1-Dμ(S))[D(Sin-S)-1λμ(S)x]=λ[μ(S)-D]μ(S)[D(Sin-S)-1λμ(S)x*]>0,
for μ-1(D)<S<S*. This implies that the trajectories of system (2) intersecting with the trajectories of system (23) pass through the trajectories of system (23) from the left to the right. In addition, the trajectory of system (23) passing through N(S*,xset) (i.e., C=0 in (24)) intersects the line S=μ-1(D) at two points K1(μ-1(D),xK1) and K2(μ-1(D),xK2), where xKi(i=1,2) satisfy that(26)lnxKi-xKix*=Λ2≜[λμ-1(D)-λ∫S*μ-1(D)‍(D/μ(ϑ))dϑ-xset+x*lnxset-λS*]x*,that is,
(27)-xKix*e-xKi/x*=-eΛ2x*.
Since Λ2=ln(xKi)-xKi/x*<lnx*-1, then we have -eΛ2/x*>-1/e. Thus, by Definition 3, we have
(28)xK1=-x*W(-1,-eΛ2x*),xK2=-x*W(0,-eΛ2x*),
and xK2<x*<xK1.
Next, we will prove thatξ1≥S* when the conditions of theorem are satisfied. Let N be the intersection point of the line x=x* and dS/dt=0, K1, and let K2 be the intersection points of the line S=μ-1(D) and the trajectory of system (23) passing through the point M. If xset≤min{x*,x¯/(1-Wf)}, then we have (1-Wf)xset≤x¯. Since dV2/dt|(2)>0, dS/dt|K1E¯>0, dx/dt|EC¯<0, dS/dt|CF¯>0, dS/dt|GB¯>0, and FG→ is semi-trivial solution of system (2). Therefore, we get a closed region G, the boundary of which consists of NK1^, K1E¯, EC¯, CF¯, FG¯, GB¯, and BN¯. In addition, there is no singularity in the region G. Therefore, it follows from the qualitative characteristics of the system that all the trajectories satisfying the conditions of theorem enter the closed region G and retain there. Then, it follows from Theorem 6 that system (2) has a period-1 solution in G. Furthermore, we have S*≤ξ1≤Sin. If ξ1≥μ-1(D), then ξ0≥μ-1(D), in this case, we have 0≤μ2<1. Else S*≤ξ1<μ-1(D)/(1-Wf), in this case, we have -1<μ2<0. Therefore, if the conditions of the theorem are satisfied, then the period-1 solution (ξ(t),η(t)) of system (2) is orbitally asymptotically stable and enjoying the property of asymptotic phase.


## 3.4. Discussion on the Period-k (k≥2) Solution

For the predator-prey model concerning IPM strategies, Tang and Cheke [23] proved that there is no periodic solution with order larger than or equal to three, except for one special case, by using the properties of the Lambert W function and Poincaré map. Moreover, They showed that the existence of an order two periodic solution implies the existence of an order one periodic solution. Next, the existence of the period-k (k≥3) solution is discussed.Theorem 12.
Suppose thatSin>μ-1(D)/(1-Wf) and λ(Sin-μ-1(D)/(1-Wf))≤xset<λ(Sin-μ-1(D)). Then, system (2) does not have period-k (k≥3) solution.Proof.
Let(S,x) be an arbitrary solution of system (2). Denote the first intersection point of the trajectory to the impulsive set M by (S1,xset) (here without loss of generality, we assume that S1>ξ1), and the corresponding consecutive points are (S2,xset), (S3,xset), …, respectively. Consequently, under the effect of impulsive function I, the corresponding points after pulse are (S1+,(1-Wf)xset),(S2+,(1-Wf)xset),(S3+,(1-Wf)xset),…. By the qualitative analysis of system (2), we know that x˙>0 for S>μ-1(D) and x˙<0 for S<μ-1(D). If there exists the jth impulse effect such that Sk+≥μ-1(D) for k≥j, then the sequence will be monotone, that is, Sj+≥Sj+1+≥⋯≥μ-1(D) or μ-1(D)≤Sj+≤Sj+1+≤⋯. In this case, system (2) has no period-k(k≥2) solution. Let l1 be the line passing the point G(Sin,0) with the slope -λ:
(29)l1:x+λS-λSin=0.
Since dl1/dt|(2)=dx/dt+λdS/dt=D(λSin-λS-x)=0, then the trajectory starting from the point on the segment CD either does not intersect l1 or retain on l1 before reaching the segment AB. Therefore, when xset≥λ(Sin-μ-1(D)/(1-Wf)), the period-1 solution determined in Theorem 6 satisfies that (1-Wf)min{SR,SM}≤ξ0<(1-Wf)(Sin-xset/λ)≤μ-1(D), which also implies that all Sk+≤μ-1(D).
The tr
ajectory of the solution(S,x) jumps to the point (S1+,(1-Wf)xset) from the point (S1,xset). By the dynamics of system (2), we have the following two sequences according to (S3,xset) lies on the right or left of (S1,xset): (a)  ⋯≤S4≤S2<ξ1<S1≤S3≤⋯; (b)  S2≤S4≤⋯≤ξ1≤⋯≤S3≤S1. Since the subsequences {S2k-1} and {S2k}   (k=1,2,…) are monotone and bounded, then there exists a limit for them, respectively. Let S_=limk→∞S2k and S¯=limk→∞S2k-1. When sequence (a) occurs, it can be shown that S_<ξ1<S¯. Therefore, system (2) has a period-2 solution. When sequence (b) occurs, we have S_≤ξ1≤S¯. If S_<ξ1<S¯, then system (2) has a period-2 solution. Else we have S_=ξ1=S¯, in this case, the trajectory tends to the period-1 solution (ξ,η). But it can be obtained that there is no order k solution (k≥3) in system (2), thus the system is not chaos by Li-Yorke chaos Theorem [44].Remark 13.
From the proof of Theorem12, it can be seen that there is a possibility to exist a period-2 solution. Since it is difficult to obtain the exact expression of the solution to (2), then it will be challenging to give the sufficient conditions under which system (2) has a period-2 solution or not.


## 4. Applied Instance

Here, we will take the Tessier kinetics model [41], that is, μ(S)=μmax(1-exp(-S/KS)), as the example of specific growth rate to verify the theoretical results. In addition, μmax=0.5[1/h], KS=0.2[g/L], and λ=Yx/Smax=0.5[g/g] are used for the demonstration of system behavior.


### 4.1. Numerical Simulations

The simulations are carried out by changing one main parameter and fixing all other parameters. The following numerical simulations are given by the programs of Maple and Matlab softwares. The time interval is set two days, that is,48[h].We assume in the following thatS0=S(0)=0.2[g/L] and x0=x(0)=2[g/L]. The vector graph of system (1) for D=0.4[1/h] and Sin=20[g/L] are shown in Figure 4, from which it can be seen that the equilibrium (0.32,9.84) is a stable node.Figure 4
Illustration of vector graph of system (1) when μmax=0.5[1/h], KS=0.2[g/L], λ=Yx/Smax=0.5[g/g], D=0.4[1/h], and Sin=20[g/L].The time series and phase portrait forxset=10[g/L] is shown in Figure 5, it can be seen that no impulse occurs when xset>λ(Sin+KSln(1-D/μmax))=9.84[g/L].The time series and phase portrait whenD=0.4[1/h], Sin=20[g/L] and xset=10[g/L].
(a)
(b)
(c)Figure6 gives the time series and phase portrait when xset=6[g/L] and Wf=0.3. It can be seen that the trajectory tends to be periodic.The time series and phase portrait whenD=0.4[1/h], Wf=0.3, Sin=20[g/L], and xset=6[g/L].
(a)
(b)
(c)Figure7 shows the different positions of the period-1 solution for D=0.4[1/h] and Wf=0.1,0.3,0.5.Dependence of the biomass and substrate concentration onWf. The phase portraits when D=0.4[1/h], Sin=20[g/L], xset=6[g/L], and Wf=0.1,0.3,0.5.
(a)
(b)
(c)Figure8 shows the different positions of the period-1 solution for Wf=0.3, Sin=20[g/L], and D=0.3[1/h],0.4[1/h],0.45[1/h].Dependence of the biomass and substrate concentration onD. The phase portraits when Wf=0.3, Sin=20[g/L], xset=6[g/L], and D=0.3[1/h],0.4[1/h],0.45[1/h].
(a)
(b)
(c)The numerical simulations are consistent with the theoretical results obtained and presented in Section3.


### 4.2. Optimization of the Bioprocess

Next, we discuss aspects of the bioprocess optimization.


#### 4.2.1. The Proposed Objective Function and the Constraints

(
1
)  The Objective Function(30)XOUT(D,Wf,Sin,xset)=1T∫0T‍Dxdt+WfxsetT,
where XOUT is the biomass productivity in the steady state.(
2
)  The Constraints(a)
0≤Wf≤Wfmax<1
Wf—the part of “clear” medium inputted into the bioreactor in each biomass oscillation cycle,Wfmax—the maximum part of “clear” medium inputted into the bioreactor in each biomass oscillation cycle.(b)
Smin≤Sin≤SCRI:
Sin—the concentration of the inputted substrate,Smin—the minimum concentration of the inputted substrate,SCRI—the critical level of the dosaged substrate concentration.(c)
0<Dminφ≤D≤Dmax=μmax(1-e-Sin/KS):
D—the dilution rate of the chemostat,Dminφ—the minimum dilution rate. If the used microorganisms growth kinetics model satisfies regularity and monotonicity, then for Sin≫Smin there is
(31)φ=1,i.e.,forD≥Dmin1,Yx/S=Yx/Smax=λ,φ=2,thatis,forD≥Dmin2,Yx/S=Yx/Smax=λ,μ(S)=μmax,whereDmin2>Dmin1,μmax—the maximum specific growth rate,KS—the substrate saturation constant.(d)
0<xset≤λ(Sin+KSln(1-D/μmax))≤xCRI:
xset—the set level of the biomass concentration in the bioreactor medium,xCRI—the critical level of biomass concentration in the bioreactor medium.The aim of the optimization is to find the maximum of the objective function (30) under constraints (a)–(d), that is,
(32)XOUT→max.


#### 4.2.2. The Analysis of the Optimization

(A) Case of a Chemostat with Impulse Effect.For D>Dmin2 and Wf>0, a continuous bioprocess (a chemostat with impulse effect) is obtained, that is,
(33)dSdt=D(Sin-S)-μmaxλxdxdt=(μmax-D)x}x<xset,ΔS=-WfSΔx=-Wfx}x=xset.
In the steady state, that is, for the period-1 solution (ξ(t),η(t)), by the second equation of system (33), we have
(34)dt=1μmax-D1ηdη.
Then, travelling along η(ξ) from the point P0(ξ0,η0), with t=tP0 to the point P1(ξ1,η1), with t=tP1 in the counterclockwise direction yields the period T, that is,
(35)T=tP1-tP0=∫η0η1‍1μmax-D1ηdη=-ln(1-Wf)μmax-D.
It is obviouse that T→0 when Wf→0. In addition, by  (34), for t∈(0,T], we have η(t)=(1-Wf)xsete(μmax-D)t, thus,
(36)∫0T‍Dη(t)dt=Dμmax-DWfxset.
Therefore, the objective function can be formulated as
(37)XOUT(2)(D,Wf,Sin,xset)=Dμmax-DWfxsetT+WfxsetT.
For given Wf=Wfset>0 and xset, we have
(38)XOUToptimal(2)=limD→μmaxXOUT(2)(D,Wfset,Sin,xset)=limD→μmaxxsetWfsetT(1-D/μmax)→-μmaxxsetWfsetln(1-Wfset),
where xset<λSCRI and Wfset∈(0,Wfmax].(B) Case of Synchronization Loss.In this case the optimization of the function (37) for any value of Wf is presented. We should find the maximum of the objective function (37) under the constraints 0≤Wf≤Wfmax, Dmin≤D≤Dmax, Smin≤Sin≤SCRI, and 0<xset<λSin.Letφ(Wf)=-Wf/ln(1-Wf). Take the first derivative of φ(Wf) with respect to Wf, we have
(39)φ'(Wf)=-ln(1-Wf)+Wf/(1-Wf)ln2(1-Wf).
Let ψ(Wf)=ln(1-Wf)-Wf/(1-Wf). Then, we have ψ(0)=0. Take the first derivative of ψ(Wf) with respect to Wf, we have
(40)ψ(Wf)=Wf(1-Wf)2>0,
then ψ(Wf)≥ψ(0)=0, which implies that φ'(Wf)≤0, that is, φ(Wf) is decreasing on Wf.Therefore, for givenxset, we have
(41)XOUToptimal(3)=limWf→0XOUToptimal(2)=limD→μmax,Wf→0XOUT(2)→μmaxxset,
where xset<λSCRI.


## 4.1. Numerical Simulations

The simulations are carried out by changing one main parameter and fixing all other parameters. The following numerical simulations are given by the programs of Maple and Matlab softwares. The time interval is set two days, that is,48[h].We assume in the following thatS0=S(0)=0.2[g/L] and x0=x(0)=2[g/L]. The vector graph of system (1) for D=0.4[1/h] and Sin=20[g/L] are shown in Figure 4, from which it can be seen that the equilibrium (0.32,9.84) is a stable node.Figure 4
Illustration of vector graph of system (1) when μmax=0.5[1/h], KS=0.2[g/L], λ=Yx/Smax=0.5[g/g], D=0.4[1/h], and Sin=20[g/L].The time series and phase portrait forxset=10[g/L] is shown in Figure 5, it can be seen that no impulse occurs when xset>λ(Sin+KSln(1-D/μmax))=9.84[g/L].The time series and phase portrait whenD=0.4[1/h], Sin=20[g/L] and xset=10[g/L].
(a)
(b)
(c)Figure6 gives the time series and phase portrait when xset=6[g/L] and Wf=0.3. It can be seen that the trajectory tends to be periodic.The time series and phase portrait whenD=0.4[1/h], Wf=0.3, Sin=20[g/L], and xset=6[g/L].
(a)
(b)
(c)Figure7 shows the different positions of the period-1 solution for D=0.4[1/h] and Wf=0.1,0.3,0.5.Dependence of the biomass and substrate concentration onWf. The phase portraits when D=0.4[1/h], Sin=20[g/L], xset=6[g/L], and Wf=0.1,0.3,0.5.
(a)
(b)
(c)Figure8 shows the different positions of the period-1 solution for Wf=0.3, Sin=20[g/L], and D=0.3[1/h],0.4[1/h],0.45[1/h].Dependence of the biomass and substrate concentration onD. The phase portraits when Wf=0.3, Sin=20[g/L], xset=6[g/L], and D=0.3[1/h],0.4[1/h],0.45[1/h].
(a)
(b)
(c)The numerical simulations are consistent with the theoretical results obtained and presented in Section3.


## 4.2. Optimization of the Bioprocess

Next, we discuss aspects of the bioprocess optimization.


### 4.2.1. The Proposed Objective Function and the Constraints

(
1
)  The Objective Function(30)XOUT(D,Wf,Sin,xset)=1T∫0T‍Dxdt+WfxsetT,
where XOUT is the biomass productivity in the steady state.(
2
)  The Constraints(a)
0≤Wf≤Wfmax<1
Wf—the part of “clear” medium inputted into the bioreactor in each biomass oscillation cycle,Wfmax—the maximum part of “clear” medium inputted into the bioreactor in each biomass oscillation cycle.(b)
Smin≤Sin≤SCRI:
Sin—the concentration of the inputted substrate,Smin—the minimum concentration of the inputted substrate,SCRI—the critical level of the dosaged substrate concentration.(c)
0<Dminφ≤D≤Dmax=μmax(1-e-Sin/KS):
D—the dilution rate of the chemostat,Dminφ—the minimum dilution rate. If the used microorganisms growth kinetics model satisfies regularity and monotonicity, then for Sin≫Smin there is
(31)φ=1,i.e.,forD≥Dmin1,Yx/S=Yx/Smax=λ,φ=2,thatis,forD≥Dmin2,Yx/S=Yx/Smax=λ,μ(S)=μmax,whereDmin2>Dmin1,μmax—the maximum specific growth rate,KS—the substrate saturation constant.(d)
0<xset≤λ(Sin+KSln(1-D/μmax))≤xCRI:
xset—the set level of the biomass concentration in the bioreactor medium,xCRI—the critical level of biomass concentration in the bioreactor medium.The aim of the optimization is to find the maximum of the objective function (30) under constraints (a)–(d), that is,
(32)XOUT→max.


### 4.2.2. The Analysis of the Optimization

(A) Case of a Chemostat with Impulse Effect.For D>Dmin2 and Wf>0, a continuous bioprocess (a chemostat with impulse effect) is obtained, that is,
(33)dSdt=D(Sin-S)-μmaxλxdxdt=(μmax-D)x}x<xset,ΔS=-WfSΔx=-Wfx}x=xset.
In the steady state, that is, for the period-1 solution (ξ(t),η(t)), by the second equation of system (33), we have
(34)dt=1μmax-D1ηdη.
Then, travelling along η(ξ) from the point P0(ξ0,η0), with t=tP0 to the point P1(ξ1,η1), with t=tP1 in the counterclockwise direction yields the period T, that is,
(35)T=tP1-tP0=∫η0η1‍1μmax-D1ηdη=-ln(1-Wf)μmax-D.
It is obviouse that T→0 when Wf→0. In addition, by  (34), for t∈(0,T], we have η(t)=(1-Wf)xsete(μmax-D)t, thus,
(36)∫0T‍Dη(t)dt=Dμmax-DWfxset.
Therefore, the objective function can be formulated as
(37)XOUT(2)(D,Wf,Sin,xset)=Dμmax-DWfxsetT+WfxsetT.
For given Wf=Wfset>0 and xset, we have
(38)XOUToptimal(2)=limD→μmaxXOUT(2)(D,Wfset,Sin,xset)=limD→μmaxxsetWfsetT(1-D/μmax)→-μmaxxsetWfsetln(1-Wfset),
where xset<λSCRI and Wfset∈(0,Wfmax].(B) Case of Synchronization Loss.In this case the optimization of the function (37) for any value of Wf is presented. We should find the maximum of the objective function (37) under the constraints 0≤Wf≤Wfmax, Dmin≤D≤Dmax, Smin≤Sin≤SCRI, and 0<xset<λSin.Letφ(Wf)=-Wf/ln(1-Wf). Take the first derivative of φ(Wf) with respect to Wf, we have
(39)φ'(Wf)=-ln(1-Wf)+Wf/(1-Wf)ln2(1-Wf).
Let ψ(Wf)=ln(1-Wf)-Wf/(1-Wf). Then, we have ψ(0)=0. Take the first derivative of ψ(Wf) with respect to Wf, we have
(40)ψ(Wf)=Wf(1-Wf)2>0,
then ψ(Wf)≥ψ(0)=0, which implies that φ'(Wf)≤0, that is, φ(Wf) is decreasing on Wf.Therefore, for givenxset, we have
(41)XOUToptimal(3)=limWf→0XOUToptimal(2)=limD→μmax,Wf→0XOUT(2)→μmaxxset,
where xset<λSCRI.


## 4.2.1. The Proposed Objective Function and the Constraints

(
1
)  The Objective Function(30)XOUT(D,Wf,Sin,xset)=1T∫0T‍Dxdt+WfxsetT,
where XOUT is the biomass productivity in the steady state.(
2
)  The Constraints(a)
0≤Wf≤Wfmax<1
Wf—the part of “clear” medium inputted into the bioreactor in each biomass oscillation cycle,Wfmax—the maximum part of “clear” medium inputted into the bioreactor in each biomass oscillation cycle.(b)
Smin≤Sin≤SCRI:
Sin—the concentration of the inputted substrate,Smin—the minimum concentration of the inputted substrate,SCRI—the critical level of the dosaged substrate concentration.(c)
0<Dminφ≤D≤Dmax=μmax(1-e-Sin/KS):
D—the dilution rate of the chemostat,Dminφ—the minimum dilution rate. If the used microorganisms growth kinetics model satisfies regularity and monotonicity, then for Sin≫Smin there is
(31)φ=1,i.e.,forD≥Dmin1,Yx/S=Yx/Smax=λ,φ=2,thatis,forD≥Dmin2,Yx/S=Yx/Smax=λ,μ(S)=μmax,whereDmin2>Dmin1,μmax—the maximum specific growth rate,KS—the substrate saturation constant.(d)
0<xset≤λ(Sin+KSln(1-D/μmax))≤xCRI:
xset—the set level of the biomass concentration in the bioreactor medium,xCRI—the critical level of biomass concentration in the bioreactor medium.The aim of the optimization is to find the maximum of the objective function (30) under constraints (a)–(d), that is,
(32)XOUT→max.


## 4.2.2. The Analysis of the Optimization

(A) Case of a Chemostat with Impulse Effect.For D>Dmin2 and Wf>0, a continuous bioprocess (a chemostat with impulse effect) is obtained, that is,
(33)dSdt=D(Sin-S)-μmaxλxdxdt=(μmax-D)x}x<xset,ΔS=-WfSΔx=-Wfx}x=xset.
In the steady state, that is, for the period-1 solution (ξ(t),η(t)), by the second equation of system (33), we have
(34)dt=1μmax-D1ηdη.
Then, travelling along η(ξ) from the point P0(ξ0,η0), with t=tP0 to the point P1(ξ1,η1), with t=tP1 in the counterclockwise direction yields the period T, that is,
(35)T=tP1-tP0=∫η0η1‍1μmax-D1ηdη=-ln(1-Wf)μmax-D.
It is obviouse that T→0 when Wf→0. In addition, by  (34), for t∈(0,T], we have η(t)=(1-Wf)xsete(μmax-D)t, thus,
(36)∫0T‍Dη(t)dt=Dμmax-DWfxset.
Therefore, the objective function can be formulated as
(37)XOUT(2)(D,Wf,Sin,xset)=Dμmax-DWfxsetT+WfxsetT.
For given Wf=Wfset>0 and xset, we have
(38)XOUToptimal(2)=limD→μmaxXOUT(2)(D,Wfset,Sin,xset)=limD→μmaxxsetWfsetT(1-D/μmax)→-μmaxxsetWfsetln(1-Wfset),
where xset<λSCRI and Wfset∈(0,Wfmax].(B) Case of Synchronization Loss.In this case the optimization of the function (37) for any value of Wf is presented. We should find the maximum of the objective function (37) under the constraints 0≤Wf≤Wfmax, Dmin≤D≤Dmax, Smin≤Sin≤SCRI, and 0<xset<λSin.Letφ(Wf)=-Wf/ln(1-Wf). Take the first derivative of φ(Wf) with respect to Wf, we have
(39)φ'(Wf)=-ln(1-Wf)+Wf/(1-Wf)ln2(1-Wf).
Let ψ(Wf)=ln(1-Wf)-Wf/(1-Wf). Then, we have ψ(0)=0. Take the first derivative of ψ(Wf) with respect to Wf, we have
(40)ψ(Wf)=Wf(1-Wf)2>0,
then ψ(Wf)≥ψ(0)=0, which implies that φ'(Wf)≤0, that is, φ(Wf) is decreasing on Wf.Therefore, for givenxset, we have
(41)XOUToptimal(3)=limWf→0XOUToptimal(2)=limD→μmax,Wf→0XOUT(2)→μmaxxset,
where xset<λSCRI.


## 5. Conclusions

In the paper the dynamic behavior of a universal biochemical reaction process with feedback control was analyzed, and it was shown that the stability of the bioprocess (i.e., the existence of the positive period-1 solution), depended on both the biomass yield and the microorganism growth rate. Furthermore, the conditions for the existence and stability of the system's period-1 solution were obtained (i.e., Theorems6 and 11). It also pointed out that the system (2) was not chaotic according to an analysis of period-2 solution existence. According to the theoretical results, the production of the microorganisms tended to be periodic. The microorganisms in the considered chemostat always kept the suitable growth rate. The biomass concentration could have been controlled to a given level between the critical and the optimal dissolved oxygen concentration. Moderation of the main bioprocess parameters in the selected biomass oscillation cycles for the Tessier kinetics model was presented.After the analysis of the system stability, bioprocess optimization was covered in the work. In the first case, the optimization of a simple continuous bioprocess (i.e., bioprocess in whichD>0 and Wf=0) was presented. In the second case of optimization, it was demonstrated that theoretically it was possible to obtain almost optimal biomass productivity without lost of synchronization of simultaneously performed bioprocesses (see (38) for small values of Wfset). This was possible, when the substrate was dosed by both continuous (i.e., D>0) and pulse (i.e., Wf>0) method. In this case, the impulsive chemostat was obtained, where the biomass productivity mainly depended on D, and the synchronization of simultaneously performed bioprocesses depended on Wf. The third case of optimization showed that during the optimization, the impulsive chemostat lost possibility of synchronization and strived for a simple continuous bioprocess (see (41)). Moreover, it was shown that improving the biomass productivity and synchronization of simultaneously performed bioprocesses were conflict criteria.



---

*Source: 105765-2012-12-18.xml*