
Citation: Ting-Hao Hsu, Tyler Meadows, LinWang, Gail S. K. Wolkowicz. Growth on two limiting essential resources in a self-cycling fermentor[J]. Mathematical Biosciences and Engineering, 2019, 16(1): 78-100. doi: 10.3934/mbe.2019004
[1] | Alma Mašić, Hermann J. Eberl . On optimization of substrate removal in a bioreactor with wall attached and suspended bacteria. Mathematical Biosciences and Engineering, 2014, 11(5): 1139-1166. doi: 10.3934/mbe.2014.11.1139 |
[2] | Tongqian Zhang, Ning Gao, Tengfei Wang, Hongxia Liu, Zhichao Jiang . Global dynamics of a model for treating microorganisms in sewage by periodically adding microbial flocculants. Mathematical Biosciences and Engineering, 2020, 17(1): 179-201. doi: 10.3934/mbe.2020010 |
[3] | Huidong Cheng, Hui Xu, Jingli Fu . Dynamic analysis of a phytoplankton-fish model with the impulsive feedback control depending on the fish density and its changing rate. Mathematical Biosciences and Engineering, 2023, 20(5): 8103-8123. doi: 10.3934/mbe.2023352 |
[4] | Zhenzhen Shi, Huidong Cheng, Yu Liu, Yanhui Wang . Optimization of an integrated feedback control for a pest management predator-prey model. Mathematical Biosciences and Engineering, 2019, 16(6): 7963-7981. doi: 10.3934/mbe.2019401 |
[5] | Tarahom Mesri Gundoshmian, Sina Ardabili, Mako Csaba, Amir Mosavi . Modeling and optimization of the oyster mushroom growth using artificial neural network: Economic and environmental impacts. Mathematical Biosciences and Engineering, 2022, 19(10): 9749-9768. doi: 10.3934/mbe.2022453 |
[6] | Meng Zhang, Kaiyuan Liu, Lansun Chen, Zeyu Li . State feedback impulsive control of computer worm and virus with saturated incidence. Mathematical Biosciences and Engineering, 2018, 15(6): 1465-1478. doi: 10.3934/mbe.2018067 |
[7] | Yan Zhang, Shujing Gao, Shihua Chen . Modelling and analysis of a stochastic nonautonomous predator-prey model with impulsive effects and nonlinear functional response. Mathematical Biosciences and Engineering, 2021, 18(2): 1485-1512. doi: 10.3934/mbe.2021077 |
[8] | Jian Meng, Bin Zhang, Tengda Wei, Xinyi He, Xiaodi Li . Robust finite-time stability of nonlinear systems involving hybrid impulses with application to sliding-mode control. Mathematical Biosciences and Engineering, 2023, 20(2): 4198-4218. doi: 10.3934/mbe.2023196 |
[9] | Yuan Tian, Sanyi Tang . Dynamics of a density-dependent predator-prey biological system with nonlinear impulsive control. Mathematical Biosciences and Engineering, 2021, 18(6): 7318-7343. doi: 10.3934/mbe.2021362 |
[10] | Agustín Gabriel Yabo, Jean-Baptiste Caillau, Jean-Luc Gouzé . Optimal bacterial resource allocation: metabolite production in continuous bioreactors. Mathematical Biosciences and Engineering, 2020, 17(6): 7074-7100. doi: 10.3934/mbe.2020364 |
The self-cycling fermentation (SCF) process can be described in two stages: In the first stage, a well-stirred tank is filled with resources and inoculated with microorganisms that consume the resources. When a threshold concentration of one or more indicator quantities is reached, the second stage is initiated. The first stage is a batch culture [5]. During the second stage, the tank is partially drained, and subsequently refilled with fresh resources before repeating the first stage.
SCF is most often applied to wastewater treatment processes, where the goal is to reduce the concentration of one or more harmful compounds [10,16]. In this application the concentration of harmful compounds is the most reasonable threshold quantity, since acceptable concentrations would typically be given by some government agency. More recently, the SCF process has been used as a means to improve production of some biologically derived compounds [20,23]. In these instances, dissolved O2 content, or dissolved CO2 content have been used as threshold quantities, since they are good indicators of when the microorganism approaches the stationary phase in its growth cycle. In both scenarios, the end goal is to maximize the amount of substrate processed by the reactor, while maintaining stable operating conditions. SCF has also been used to culture synchronized microbial cultures [17], where stability of the operating conditions is much more important than output of the reactor. It is the first scenario that we model, i.e. we consider the case that the microbial population is used to reduce two harmful compounds to an acceptable level.
Assuming that the time taken to empty and refill the tank is negligible, we can model the SCF process using a system of impulsive differential equations. Smith and Wolkowicz [18] used this approach to model the growth of a single species with one limiting resource. Fan and Wolkowicz [7] extended this model to include the possibility that the resource is also limiting at large concentrations. Córdova-Lepe, Del Valle, and Robledo [6] also modeled single species growth in the SCF process, but used impulse dependent impulse times instead of the state dependent impulses used by the other models. For references on the theory of impulsive differential equations, see e.g. [1,2,3,9,15].
When there are two (or more) resources in limited supply, it is important to think about how the resources interact to promote growth. If all of the resources fulfill the same requirement for growth, the resources are called substitutable resources. For instance, both glucose and fructose are carbon sources for many bacteria, and can fulfill the same purpose in bacterial growth. If all of the resources are required in some way for growth, and the bacteria would die out if any of them were missing, they are called essential resources. For instance, both carbon and nitrogen are required for growth of many bacteria, but glucose cannot be used as a nitrogen source, so some other compound such as nitrate is required. Growth and competition with two essential resources has been studied in the chemostat [4], in the chemostat with delay [12], and in the unstirred chemostat [24]. In all of the aforementioned studies, the interaction of essential resources is through Liebig's law of the minimum [22]. To illustrate the law of the minimum, consider a barrel with several staves of unequal length. Growth is limited by the resource in shortest supply in the same way that the capacity of the barrel is limited by the length of the shortest stave.
In this paper we investigate the dynamics of the self-cycling fermentation process in a semi-batch culture with two essential resources that are assumed to be pollutants. The goal is to reduce both pollutant concentrations to acceptable levels. In section 2 we introduce the model. In section 3 we analyze the system of ordinary differential equations (ODEs) associated with the model introduced in section 2. In section 4 we analyze the system of impulsive differential equations introduced in section 2, and obtain our main results: Theorem 4.6, which gives necessary and sufficient conditions for the existence of a unique periodic orbit, and Theorem 4.14, which summarizes all of the possible long term dynamics of the model. In section 5 we demonstrate numerically that the emptying/refilling fraction can be used to maximize the output of the SCF process. In section 6 we summarize our results and discuss the implications of our analysis. All figures were produced using Matlab [13].
For a given function y(t) and time τ, using the standard notation for impulsive equations we denote Δy(τ)=y(τ+)−y(τ−), where
y(τ+)≡limt→τ+y(t) and y(τ−)≡limt→τ−y(t). |
Our model takes the form
ds1(t)dt=−1Y1min{f1(s1(t)),f2(s2(t))}x(t),ds2(t)dt=−1Y2min{f1(s1(t)),f2(s2(t))}x(t),dx(t)dt=(−D+min{f1(s1(t)),f2(s2(t))})x(t),}t≠tkΔs1(tk)=−rs1(t−k)+rsin1,Δs2(tk)=−rs2(t−k)+rsin2,Δx(tk)=−rx(t−k),}t=tk | (1) |
where tk are the times at which
either(s1(tk)=ˉs1,s2(tk)≤ˉs2)or(s1(tk)≤ˉs1,s2(tk)=ˉs2). | (2) |
Here, t denotes time. The variables si, i=1,2 denote the concentration of the limiting resources (assumed to be pollutants) in the fermentor as a function of t, with associated parameters Yi, the cell yield constants, sini, the concentrations of each limiting resource in the medium added to the tank at the beginning of each new cycle, and ˉsi the threshold concentrations of limiting resource that trigger the emptying and refilling process. Since we are considering the scenario where both s1 and s2 are pollutants, the emptying and refilling process is only triggered when both concentrations reach the acceptable levels ˉs1 and ˉs2 set by some environmental protection agency. The variable x denotes the biomass concentration of the population of microorganisms that consume the resource at time t, assumed to have death rate D. The emptying/refilling fraction is denoted by r. It is assumed that D>0, 0<r<1 and for i=1,2, Yi>0, and sini>ˉsi>0.
We call the times tk>0, impulse times, and when they exist they form an increasing sequence that we denote {tk}Nk=1. If (2) is satisfied at t=0 or si(0)<ˉsi, i=1,2, then we assume that there is an immediate impulse at time t=0. We consider the process to be successful if N=∞ and the time between impulses, tk−tk−1 remains bounded. We consider the process a failure if either there are a finite number of impulses, and hence N is finite, or if the time between impulses becomes unbounded. The two resources are assumed to be limiting essential resources (see e.g., Tilman [21] or Grover [8]) also called complementary resources (see Leon and Tumpson [11]), and as in those studies we use Liebig's law of the minimum [22] to model the uptake and growth of the microbial population.
We assume that each response function fj(s), j=1,2, in (1) satisfies:
(ⅰ) fj:R+→R+ is continuously differentiable;
(ⅱ) fj(0)=0 and f′j(s)>0 for s>0;
Define λi, i=1,2, to be the value of each resource that satisfies fi(λi)=D, and refer to each λi as a "break-even concentration". If fj is bounded below D, then we define the corresponding λj=∞.
Between two consecutive impulses the system is governed by a system of ordinary differential equations (ODE) that models a batch fermentor [5],
ds1(t)dt=−1Y1min{f1(s1(t)),f2(s2(t))}x(t),ds2(t)dt=−1Y2min{f1(s1(t)),f2(s2(t))}x(t),dx(t)dt=(−D+min{f1(s1(t)),f2(s2(t))})x(t). | (3) |
We will refer to system (3) as the associated ODE system.
First we show that system (3) is well-posed.
Proposition 3.1. Given any positive initial conditions (s1(0),s2(0),x(0)), the solution (s1(t),s2(t),x(t)) of (3) is defined for all t≥0 and remains positive.
Furthermore, limt→∞(s1(t),s2(t),x(t)) exists, is initial condition dependent, limt→∞x(t)=0, and limt→∞si(t)<λi for at least one i∈{1,2}.
Proof. Since the vector field in (3) is locally Lipschitz, the positivity of (s1(t),s2(t),x(t)) follows from the standard theory for the existence and uniqueness of solutions of ODEs (see e.g., [14]). Also observe that
ddt(x(t)+Y12s1(t)+Y22s2(t))=−Dx(t)<0, |
so the solution (s1(t),s2(t),x(t)) exists and is bounded for t in [0,∞).
From (3), s′i(t)<0, i=1,2. By the positivity of si(t), limt→∞si(t)≥0 exists. Denote limt→∞si(t)=s∗i, i=1,2.
We claim that min{f1(s∗1),f2(s∗2)}<D. Suppose not. Then, min{f1(s∗1),f2(s∗2)}≥D. Since both s1(t) and s2(t) are strictly decreasing functions, it follows that min{f1(s1(t)),f2(s2(t))}>D for all t≥0. From the equation of x′(t) in (3), x(t) would be a strictly increasing function. Then, s′1(t)<−DY1x(0) for all t≥0, contradicting the positivity of s1(t). Hence, s∗i<λi for at least one i∈{1,2}.
Since min{f1(s∗1),f2(s∗2)}<D, define γ=−D+min{f1(s∗1),f2(s∗2)}<0. By the equation for x′(t) in (3), x′(t)<γ2x(t)<0, for all sufficiently large t. Hence x(t)→0 as t→∞.
Note that from (3), Y1ds1dt=Y2ds2dt. Define
R12=Y2Y1andR21=Y1Y2=1R12. | (4) |
Then, every trajectory of (3) satisfies ds2ds1=R21.
Lemma 3.2. Let (s1(t),s2(t),x(t)) be a solution of (3) on an interval t∈[t0,t1] with positive initial conditions. Then,
s1(t)=s1(t0)+R12(s2(t)−s2(t0)), | (5) |
or equivalently
s2(t)=s2(t0)+R21(s1(t)−s1(t0)), | (6) |
x(t1)−x(t0)=Y1∫s1(t0)s1(t1)(1−Dmin{f1(v),f2(s2(t1)+R21(v−s1(t1)))})dv, | (7) |
or equivalently
x(t1)−x(t0)=Y2∫s2(t0)s2(t1)(1−Dmin{f1(s1(t1)+R12(v−s2)),f2(v)})dv. | (8) |
Proof. Solving the separable ODE, ds1ds2=R12, yields (5), and solving ds2ds1=R21, yields (6). Dividing the s′1(t) equation in (3) by the x′(t) equation, substituting for s2(t) using (6), and then integrating both sides, yields (7). We obtain (8) similarly.
First we visualize solutions of (1) in the s1-s2 plane as illustrated in Figure 1. Given any solution (s1(t),s2(t),x(t)) of (1) with positive initial conditions, let t0=0 and let tk, k∈N, denote the kth impulse time if it exists.
Let u±k=(s1(t±k),s2(t±k)). By (6), the trajectory of (s1(t),s2(t)), t∈(tk,tk+1), is a line segment with slope R21 and endpoints u+k and u−k+1. The conditions for impulses to occur are given in (2). Therefore, each point u−k lies in the following union of the two horizontal and vertical line segments:
Γ−≡{(s1,ˉs2):s1∈[0,ˉs1]}∪{(ˉs1,s2):s2∈[0,ˉs2]}. |
Define
uin=(sin1,sin2)andˉs+i=(1−r)ˉsi+rsini, i=1,2. |
By the definition of Δsi given in (1),
u+k=(1−r)u−k+ruin. | (9) |
This implies that each u+k lies in the following union of horizontal and vertical line segments:
Γ+≡{(s1,ˉs+2):s1∈[0,ˉs+1]}∪{(ˉs+1,s2):s2∈[0,ˉs+2]}. |
Therefore, if impulses occur indefinitely, then the total trajectory of (s1(t),s2(t)), t∈[t1,∞), is a countable union of line segments with slope R21 and endpoints in Γ−∪Γ+, (i.e., u+k∈Γ+ and u−k∈Γ−.)
For any positive solution (s1(t),s2(t),x(t)) of (1) with (s1(0),s2(0)) lying between the coordinate axes and Γ−, i.e., 0≤s1(0)≤ˉs1 and 0≤s2(0)≤ˉs2, an impulse occurs immediately at t=0, and so, after at most a finite number of impulses, (s1(0+),s2(0+)) lies above or to the right of Γ−. In the rest of this section, we therefore assume (s1(0),s2(0)) lies above or to the right of Γ−, i.e., si(0)>ˉsi, for at least one i∈{1,2}.
The following proposition asserts that system (1) does not exhibit the phenomenon of beating. That is, the system possesses no solution with impulse times that form an increasing sequence with a finite accumulation point.
Proposition 4.1. Assume that (s1(t),s2(t),x(t)) is a positive solution of (1) with an infinite number of impulse times, {tk}∞k=1. Then limk→∞tk=∞.
Proof. bBetween impulses, s1 and s2 are strictly decreasing for all x(t)>0, and therefore we can solve the first equation in (3) for the time between impulses:
tk+1−tk=Y1∫s1(tk)s1(tk+1)1min{f1(v),f2(s2(tk+1)+R21(v−s1(tk+1))}Xk(v)dv | (10) |
where Xk(v) is defined by Xk(s1(t))=x(t) for t∈(tk,tk+1).
To show that {tk}∞k=1 has no finite accumulation point, it suffices to show that there exists positive constants, M1, M2 and m, such that
s1(t+k)−s1(t−k+1)>mfor k=1,2,… |
and
min{f1(s1(t)),f2(s2(t))}<M1,x(t)<M2,for t≥0, |
since then the difference tk+1−tk is greater than Y1mM1M2.
Since s1(t−k)≤ˉs1, s2(t−k)≤ˉs2, and either s1(t+k)=ˉs+1 or s2(t+k)=ˉs+2, we have
eithers1(t+k)−s1(t−k+1)≥ˉs+1−ˉs1ors2(t+k)−s2(t−k+1)≥ˉs+2−ˉs2. |
From (5) it follows that
s1(t+k)−s1(t−k+1)≥min{ˉs+1−ˉs1,R12(ˉs+2−ˉs2)}≡m. |
Since after the first impulse occurs, (s1(t),s2(t)) is bounded by Γ+, the existence of M1 follows from the continuity of f1 and f2. By (7), there exists M0>0 such that
x(t)<x(t+k)+M0,for t∈(tk,tk+1). | (11) |
From the relation that x(t+k)=rx(t−k), we obtain
x(t+k+1)<r(x(t+k)+M0),for k=1,2,… |
By the comparison principle applied to x(tk) and the sequence {yk} defined by y0=x(0) and yk+1=r(yk+M0), k=1,2,…,
lim supk→∞x(t+k)≤limk→∞yk=rM01−r. | (12) |
The existence of M2 follows from (11) and (12).
Define Ω1 to be the set of points (s01,s02) such that, for some x0>0, the forward trajectory of the solution of (3) with initial value (s01,s02,x0) intersects Γ−. Then the boundary of Ω1 is the union of Γ− and the two lines of slope R21 passing through (ˉs1,0) and (0,ˉs2), respectively. Then,
Ω1={(s1,s2):s1>ˉs1ors2>ˉs2,ands2−ˉs2<R21s1,s2>R21(s1−ˉs1)}. |
Define Ω0 to be the open set complementary to Ω1 in the first quadrant above and to the right of Γ−. Then,
Ω0={(s1,s2):s1>ˉs1ors2>ˉs2,ands2−ˉs2>R21s1ors2<R21(s1−ˉs1)}. |
Remark 4.2. The sets Ω0 and Ω1 do not include the marginal cases where (s1(0),s2(0)) lies on the lines of slope R21 passing through (ˉs1,0) or (0,ˉs1). If (s1(0),s2(0)) lies on one of these lines and si(0)>0 for i=1,2, then no impulses occur, since the solution curve does not reach Γ− in finite time. If (s1(0),s2(0))=(ˉs1,0) or (s1(0),s2(0))=(0,ˉs2), then an impulse occurs immediately, and (s1(0+),s2(0+)) may be in either Ω0 or Ω1, depending on the location of uin.
In the following case, the fermentation process fails.
Lemma 4.3. If (s1(t),s2(t),x(t)) is a solution of (1) with (s1(0),s2(0))∈Ω0, then no impulses occur.
Proof. Since Ω0 is complementary to Ω1, (s1(t),s2(t))∉Γ− for any t≥0, and hence no impulses occur.
Next we define a Lyapunov-type function
V(s1,s2)=Y2(sin2−s2)−Y1(sin1−s1). | (13) |
Then,
Ω1={(s1,s2):s1>ˉs1ors2>ˉs2,andV(0,ˉs2)<V(s1,s2)<V(ˉs1,0)} |
and
Ω0={(s1,s2):s1>ˉs1ors2>ˉs2,andV(s1,s2)<V(0,ˉs2)orV(s1,s2)>V(ˉs1,0)}. |
Note that the level sets of V are straight lines with slope R21 in the s1-s2 plane. For any fixed value of s1, V(s1,s2) is a decreasing function of s2, and for any fixed value of s2, V(s1,s2) is an increasing function of s1.
Lemma 4.4. Assume (s1(t),s2(t),x(t)) is a solution of (1) with positive initial conditions. Let t0=0 and tk, k∈N, denote the kth impulse time, if it exists; otherwise set tk=∞. Define V(t)=V(s1(t),s2(t)). Then,
(i) ddtV(t)=0 for t∈(tk,tk+1),
(ii) V(t+k)=(1−r)V(t−k), if tk<∞.
Proof. (ⅰ) By the equations for dsi/dt in (3),
ddtV(t)=(Y2Y2−Y1Y1)min{f1(s1(t)),f2(s2(t))}x(t)=0. |
(ⅱ) Substituting (9) into (13),
V(t+k)=(1−r)(s1(t−k)−sin1,s2(t−k)−sin2)⋅(Y1,−Y2)=(1−r)V(t−k), |
where ⋅ denotes the inner product in R2.
By the definition of V(s1,s2) in (13) and Lemma 4.4, the line
{(s1,s2):V(s1,s2)=0}={(sin1+vY2,sin2+vY1):v∈R} | (14) |
is invariant under (3) for all x(0). By symmetry, we may assume the point (ˉs1,ˉs2) lies on or above this invariant line, i.e., V(ˉs1,ˉs2)≤0, or
sin2−ˉs2≤R21(sin1−ˉs1). | (15) |
Next we show that in the case that (sin1,sin2) lies in Ω0, once again the fermentation process is doomed to fail.
Lemma 4.5. If (sin1,sin2)∈Ω0, then every solution of (1) with positive initial conditions has at most finitely many impulses.
Proof. Note that, V(sin1,sin2)=0, so the condition (sin1,sin2)∉Ω1 implies that either V(0,ˉs2)>0 or V(ˉs1,0)<0. Since V(0,ˉs2)<V(ˉs1,ˉs2)≤0, by assumption (15), V(ˉs1,0)<0.
We proceed using proof by contradiction. Suppose that a solution (s1(t),s2(t),x(t)) of (1) with positive initial conditions has infinitely many impulses. Denote the impulse times by t1<t2<⋯. By Lemma 4.4(ⅱ), V(s1(t−k),s2(t−k))→0 as k→∞. Hence, V(s1(t−k),s2(t−k))>V(ˉs1,0) for some k≥1. By Lemma 4.3, no more impulses can occur.
In the case of (sin1,sin2)∈Ω1, under assumption (15), the line given by (14) intersects Γ− at the point (ˉs1,ˆs2) given by
ˆs2=sin2−R21(sin1−ˉs1). |
Define L to be the portion of the line given by (14) from the point (ˉs1,ˆs2) to its image via the impulsive map, namely
L={(s1,s2):ˉs1≤s1≤ˉs+1,V(s1,s2)=0}={(s1,s2):ˉs1≤s1≤ˉs+1,s2=ˆs2+R21(s1−ˉs1)}. |
Then L is invariant under (1) for all x(0).
Next we investigate under what conditions the reactor has a periodic solution and the process has the potential to succeed.
We regard the emptying/refilling fraction r∈(0,1) as a variable. Without loss of generality, from now on we assume that (ˉs1,ˉs2) lies on or above L. Otherwise, from (15), by symmetry we can relabel the resources. Therefore, the right endpoint of L is (ˉs+1(r),ˆs+2(r)) given by
ˉs+1(r)=(1−r)ˉs1+rsin1andˆs+2(r)=(1−r)ˆs2+rsin2. |
By (7), the net change in x(t) over one cycle with impulse at (ˉs1,ˆs2) is
μ(r)=Y1∫ˉs+1(r)ˉs1(1−Dmin{f1(v),f2(ˆs2+R21(v−ˉs1))})dv. | (16) |
We prove the following theorem concerning the existence and the uniqueness of periodic solutions.
Theorem 4.6. Assume (sin1,sin2)∈Ω1. If r∈(0,1) and μ(r)>0, then system (1) has a periodic orbit that is unique up to time translation and has one impulse per period. On a periodic orbit, x(t+k)=(1−r)rμ(r) and x(t−k)=1rμ(r), for all k∈N.
If μ(r)≤0, then system (1) has no periodic orbits.
See Figure 2 for an illustration of the case with μ(r)<0.
Proof. Suppose there is a positive periodic solution (s1(t),s2(t),x(t)). Since (3) has no periodic orbits, the solution has at least one impulse. By periodicity there are infinitely many impulses. Denote the impulse times by t1<t2<⋯ and the number of impulses within a period by N. Then x(t±N+k)=x(t±k) for all k∈N. By (1),
x(t−k+1)=x(t+k)+μ(r)andx(t+k)=(1−r)x(t−k). |
Thus,
x(t−k+1)=(1−r)x(t−k)+μ(r). | (17) |
If x(t−1)>μ(r) (resp. x(t+1)<μ(r)), then from (17) it can be shown by induction that x(t−k) is a strictly decreasing (resp. strictly increasing) sequence, contradicting x(t−N+k)=x(t−k). Hence, on any periodic orbit there is only one impulse, and it follows from (17), on a periodic orbit x(t−k)=1rμ(r) for all k∈N. Therefore, μ(r)>0, since the solution must be positive, and if a periodic orbit exists it is unique.
If μ(r)>0, the solution of (1) with initial condition (ˉs+1,ˆs+2,1−rrμ(r)) is periodic, since if x(t−k+1)=1rμ(r), then x(t+k+1)=1−rrμ(r) and if x(t+k)=1−rrμ(r), then x(t−k+1)=1rμ(r).
Proposition 4.7. If μ(1)>0, then there exists a unique r∗∈[0,1) such that μ(r)>0 for all r∈(r∗,1] and μ(r)≤0 for all r∈(0,r∗].
Proof. Let
r∗=max{r∈[0,1]:μ(v)≤0∀v∈[0,r]}. | (18) |
Note that r∗ is well-defined, since μ(r) is continuous and μ(0)=0. By definition, μ(r)≤0 for all r∈[0,r∗]. Since μ(r) is continuous, if μ(1)>0 then r∗<1. Furthermore, μ(r)>0 for all r∈(r∗,1], since f1(s1) and f2(s2) are monotone increasing, and so the integrand in (16) with v=ˉs+1(r∗) must be positive. Otherwise, r∗ can be increased, contradicting definition (18). By the monotonicity of f1(s1) and f2(s2), the integrand in (16) with v=ˉs+1(r) remains positive for r∗<r≤1. Hence, μ(r)>0 for r∈(r∗,1].
Remark 4.8. If λ1≤ˉs1 and λ2≤ˆs2, then μ(r)>0 for all r∈(0,1), i.e., r∗=0, because the integrand in (16) is then positive for all v∈(ˉs1,sin).
In this subsection we fix an r∈(r∗,1), where r∗ is the number given in Theorem 4.6. Hence, μ(r)>0 and a unique periodic orbit exists.
For each point (s1,s2)∈Ω1, we denote by π−(s1,s2), the point of intersection of the line through (s1,s2) of slope R21 with Γ−, i.e.,
π−(s1,s2)={(s1+R12(ˉs2−s2),ˉs2),ifV(s1,s2)≤V(ˉs1,ˉs2),(ˉs1,s2+R21(ˉs1−s1)),ifV(s1,s2)≥V(ˉs1,ˉs2). |
For each point (s1,s2)∈Γ−, we denote by π+(s1,s2) the pre-image of (s1,s2) in Γ+ under π−. Hence, π+(s1,s2)∈Γ+ satisfies π−(π+(s1,s2))=(s1,s2) for all (s1,s2)∈Γ−. Let g:Γ−→Γ+, be the image of (s1,s2)∈Γ− under the impulsive map, i.e.,
g(s1,s2)=(1−r)(s1,s2)+r(sin1,sin2). |
For each point (s01,s02)∈Ω1, let I(s01,s02) denote the net change in x(t), over the time interval from t=0 to the first impulse time. Therefore, for a solution (s1(t),s2(t),x(t)) of (1) satisfying (s1(0),s2(0))=(s01,s02) that has at least one impulse, by Lemma 3.2,
I(s01,s02)={Y2∫s02ˉs2(1−Dmin{f1(s01+R12(v−s02),f2(v)})dv,ifV(s1,s2)≤V(ˉs1,ˉs2),Y1∫s01ˉs1(1−Dmin{f1(v),f2(s02+R21(v−s01)})dv,ifV(s1,s2)≥V(ˉs1,ˉs2). |
Note that μ(r)=I(π+(ˉs1,ˆs2)).
If (s1(t),s2(t),x(t)) is a solution of (3) with (s1(t0),s2(t0))=(s01,s02)∈Γ+ and (s1(t1),s2(t1))=(s11,s12)∈Γ− for some t0<t1, then (s01,s02)=π+(s11,s12) and the net change in x(t) over the time interval [t0,t1] is I(π+(s11,s12)). Since V(ˉs1,ˉs2)<V(ˉs+1,ˉs+2)<V(sin1,sin2), from V(sin1,sin2)=V(ˉs1,ˆs2) there exists ˜s2∈(ˆs2,ˉs2) such that V(ˉs1,˜s2)=V(ˉs+1,ˉs+2). Define
Γ−A={(ˉs1,s2):0<s2≤˜s2:I(π+(s1,s2))>0}. | (19) |
Lemma 4.9. Assume (sin1,sin2)∈Ω1 and μ(r)>0. Then,
Γ−A={ˉs1}×(s2♯,˜s2] | (20) |
for some s2♯∈(0,ˆs2).
Proof. By assumption (15), for any s2∈[0,˜s2],
I(π+(ˉs1,s2))=Y1∫ˉs+1ˉs1(1−Dmin{f1(v),f2(s2+R21(v−ˉs1)})dv. |
Let
Λ={s2∈(0,˜s2):I(π+(ˉs1,s2))>0}. |
By the monotonicity of f1(s1) and f2(s2), since I(π+(ˉs1,ˆs2))=μ(r)>0, Λ is an interval containing ˆs2 with right endpoint ˜s2. Hence, Γ−A takes the form (20) for some s2♯∈(0,ˆs2).
Let Ω1A be the set of points (s01,s02) such that, for some x0>0, the forward trajectory of the solution of (3) with initial value (s01,s02,x0) passes through Γ−A. Then,
Ω1A={(s1,s2)∈Ω1:π−(s1,s2)∈Γ−A}. | (21) |
In the case μ(r)>0, by (20),
Ω1A={(s1,s2)∈Ω1:V−<V(s1,s2)<V+}, | (22) |
where V−=V(ˉs1,˜s2) and V+=V(ˉs1,s2♯).
Lemma 4.10. Assume (sin1,sin2)∈Ω1 and μ(r)>0. Let (s1(t),s2(t),x(t)) be a solution of (1) with x(0)>0 and
(s1(0),s2(0))∈Ω1A. |
The solution converges to the unique periodic solution given by Theorem 4.6, if and only if x(0)>−I(s1(0),s2(0)). If x(0)≤−I(s1(0),s2(0)), then no impulses occur.
Proof. If x(0)≤−I(s1(0),s2(0)), then by Lemma 3.2 the value of x(t) approaches 0 before any impulses occur.
If x(0)>−I(s1(0),s2(0)), then the first impulse occurs at some finite time t1>0. The condition (s1(0),s2(0))∈Ω1A implies that (s1(t+1),s2(t+1))∈Ω1A. Hence, I(s1(t+1),s2(t+1))>0. This implies that the net change of x(t) is positive over any time interval from t1 to a time before the next impulse. Hence, another impulse occurs at some finite time t2>t1. Inductively, it follows that impulses occur indefinitely. By Lemma 4.4, limt→∞V(s1(t),s2(t))=0. Hence,
(s1(t−k),s1(t−k))→(ˉs1,ˆs2)and(s1(t+k),s1(t+k))→(ˉs+1,ˆs+2). |
Therefore, by Lemma 3.2 and the relation I(π+(ˉs1,ˆs2))=μ(r),
limk→∞(x(t−k+1)−x(t+k))=μ(r). |
On the other hand, the impulsive map in (1) gives limk→∞x(t+k)−(1−r)x(t−k)=0. This gives
limk→∞(x(t−k+1)−(1−r)x(t−k))=μ(r), |
which implies limk→∞x(t−k)=1rμ(r) and limk→∞x(t+k)=1−rrμ(r). We conclude that the solution converges to the periodic orbit.
Corollary 4.11. If (sin1,sin2)∈Ω1 and μ(r)>0, then all solutions to (1) with (s1(0),s2(0))=(sin1,sin2) and x(0)>0 converge to the periodic orbit given by Theorem 4.6.
Proof. By the definitions of Γ−A and Ω1A given in (19) and (21), respectively, (sin1,sin2)∈Ω1A. Since I(sin1,sin2)>I(π+(ˉs1,ˆs2))=μ(r)>0, the desired result follows from Lemma 4.10.
For each (s01,s02)∈Ω1, let N0=N0(s01,s02) be the smallest positive integer such that
(g∘π−)N0(s01,s02)∈Ω1A. | (23) |
In particular, N0(s01,s02)=1 for all (s01,s02)∈Ω1A. For (s01,s02)∈Ω1∖Ω1A, by the identities V(g(s1,s2))=(1−r)V(s1,s2) and V(π−(s1,s2))=V(s1,s2),
V((g∘π−)n(s01,s02))=(1−r)nV(s01,s02),n=1,2,⋯. | (24) |
If V(s01,s02)<0, then by (22), the condition (s01,s02)∉Ω1A is equivalent to V(s01,s02)≤V−. Thus, by (23), condition (24), is equivalent to
(1−r)N0V(s01,s02)>V−. |
Similarly, if V(s01,s02)>0 and (s1(0),s2(0))∈Ω1A, then condition (24) is equivalent to
(1−r)N0V(s01,s02)<V+. |
Hence,
N0(s01,s02)={⌈ln(V(s01,s02)/V−)−ln(1−r)⌉,ifV(s01,s02)<V−,⌈ln(V(s01,s02)/V+)−ln(1−r)⌉,ifV(s01,s02)>V+, | (25) |
where ⌈y⌉ is the least integer greater than or equal to y.
For any solution (s1(t),s2(t),x(t)) of (1) with (s1(0),s2(0))=(s01,s02)∈Ω1, and x(0)>0,
x(t−1)=x(0)+I(s01,s02). |
Note that x(t+k)=(1−r)x(t−k) by the impulsive map in (1), and x(t−k+1)=x(t+k)+I(s1(t+k),s2(t+k)) by Lemma (3.2). Hence, for any n=2,3,…, the left limit of x(t) at the nth impulse, if it exists, equals
x(t−n)=(1−r)x(t−n−1)+I((g∘π−)n−1(s01,s02)), |
where t0=0, and tk, k≥1, is the kth impulse time. By induction,
x(t−n)=(1−r)n−1x(0)+n∑k=1(1−r)n−kI((g∘π−)k−1(s01,s02)). |
Thus the condition x(t−n)>0 is equivalent to
x(0)>−n∑k=1(1−r)−(k−1)I((g∘π−)k−1(s01,s02)). |
We define X(s01,s02) to be the least value so that if x(0)>X(s01,s02) then (s1(t−∗),s2(t−∗))∈Γ−A for some t∗>0. Hence,
X(s01,s02)=−(min1≤n≤N0(s0,s1)n∑k=1(1−r)−(k−1)I((g∘π−)k−1(s01,s02))). | (26) |
In particular,
X(s01,s02)=−I(s01,s02)if N0(s01,s02)=1. |
The following proposition extends Lemma 4.10.
Proposition 4.12. Assume (sin1,sin2)∈Ω1 and μ(r)>0. Let (s1(t),s2(t),x(t)) be a solution of (3) with (s1(0),s2(0))∈Ω1 and x(0)>0.
(i) If x(0)≤X(s1(0),s2(0)), then there are at most N0(s01,s02)−1 impulses.
(ii) If x(0)>X(s1(0),s2(0)), then the solution converges to the unique periodic orbit given by Theorem 4.6.
Proof. (ⅰ) Suppose x(0)≤X(s1(0),s2(0)) and the solution has at least N0=N0(s01,s02) impulses. Denote the first N0 impulse times by t1<t2<⋯<tN0. Then, by Lemma 3.2 and the definition of X(s1,s2), for some k∈{1,⋯,N0},
x(t−k)=x(0)−X(s1(0),s2(0))≤0, |
contradicting the positivity of the solution.
(ⅱ) If x(0)>X(s1(0),s2(0)), then the solution has at least N0 impulses. Denote the N0th impulse time by tN0. Then,
(s1(t+N0),s2(t+N0))=(g∘π−)N0(s1(0),s2(0))∈Ω1A. |
Since (s1(t+N0),s2(t+N0))∈Γ+, by the definition of Ω1A, I(s1(t+N0),s2(t+N0))>0. Hence, the result follows from Lemma 4.10.
Example 1. Consider (1) with the Monod functional responses fi(si)=misiai+si, i=1,2, and parameters (m1,m2,a1,a2)=(2,2,1.9,0.3), (Y1,Y2,D)=(4,1.9,0.5), and (ˉs1,ˉs2,sin1,sin2,r)=(0.6,0.5,1,1,0.4).
We compute the following quantities using their definition.
(ˉs+1,ˉs+2)=(0.76,0.7),˜s2≈0.36,ˆs2≈0.16,μ(r)≈0.03,V−≈−0.39. |
Taking the initial values (s01,s02)=(0.23,0.6), we have V(s01,s02)=−2.32. Then,
N0(s01,s02)=⌈ln(V(s01,s02)/V−)−ln(1−r)⌉=⌈3.4908⌉=4. |
The approximated values of I((g∘π−)n(s01,s02)), 1≤n<N0, are as follows.
n0123(1−r)−nI((g∘π−)n(s01,s02))−0.2970−0.1785−0.04410.0846 |
Therefore X(s01,s02)≈0.2970+0.1785+0.0441=0.5196.
In Figures 3 and 4, the initial data satisfies x(0)=0.5<X(s01,s02) and x(0)=0.53>X(s01,s02), respectively. By Proposition 4.12 the fermentation succeeds only in the latter case.
If (sin1,sin2)∈Ω1 and μ(r)≤0, then, by Theorem 4.6, system (1) has no periodic solution. The following proposition asserts that the fermentation fails in this case.
Proposition 4.13. Assume (sin1,sin2)∈Ω1.
(i) If μ(r)<0, then for every solution of (1) with positive initial conditions, only finitely many impulses occur.
(ii) If μ(r)=0, then for every solution of (1) with positive initial conditions, either only finitely many impulses occur, or the time between impulses tends to infinity.
Proof. Let (s1(t),s2(t),x(t)) be a solution of (1) with positive initial conditions. Suppose the solution has infinitely many impulses. Denote the impulse times by t1<t2<⋯. Then by Lemma 4.4, (s1(t−k),s2(t−k))→(ˉs1,ˆs2) and (s1(t+k),s2(t+k))→(ˉs+1,ˆs+2) as k→∞.
(ⅰ) In the case μ(r)<0, by Lemma 3.2 and the relation I(π+(ˉs1,ˆs2))=μ(r),
limk→∞(x(t−k+1)−x(t+k))=μ(r). |
(ⅱ) On the other hand, the impulsive map in (1) gives x(t+k)=(1−r)x(t−k). This implies limk→∞x(t−k)=1rμ(r)<0, contradicting to the positivity of the solution.
In the case μ(r)=0, by Lemma 3.2 and the relation I(π+(ˉs1,ˆs2))=μ(r)=0,
limk→∞(x(t−k+1)−x(t+k))=0. |
By the relation x(t+k)=(1−r)x(t−k), it follows that limk→∞x(t±k)=0. Hence, the trajectory of (s1(t),s2(t),x(t)), t∈(tk,tk+1), approaches the heteroclinic orbit of (3) from (ˉs+1,ˆs+2,0) to (ˉs1,ˆs2,0). This implies limk→∞(tk+1−tk)=∞.
By (25), the function N0(s1,s2) of (s1,s2)∈Ω1 has an upper bound
ˉN=max{N0(0,ˉs2),N0(ˉs1,0)}. |
We summarize our results as follows.
Theorem 4.14. Consider system (1).
(i) If (sin1,sin2)∈Ω0, then every solution has at most finitely many impulses.
(ii) If (sin1,sin2)∈Ω1 and μ(r)≤0, then the fermentation fails in the sense that for every solution with positive initial conditions, either only finitely many impulses occur, or the time between impulses tends to infinity.
(iii) If (sin1,sin2)∈Ω1 and μ(r)>0, then there is a unique periodic orbit. Moreover, for any solution (s1(t),s2(t),x(t)), with positive initial conditions, the number of impulse times is either infinite or is less than ˉN. The case with infinitely many impulses occurs if and only if
(s1(0),s2(0))∈Ω1andx(0)>X(s1(0),s2(0)). |
Proof. The theorem follows from Lemmas 4.3 and 4.5, and Propositions 4.12 and 4.13.
In the following Corollary, we consider a case in which we are guaranteed that X(s01,s02)≤0. In this case, by Proposition 4.12, it follows that any solution of (1) with (s1(0),s2(0))=Ω1 and x(0)>0 converges to the periodic orbit.
If ˉs1>λ1 and ˉs2>λ2, we define Ωλ to be the region in the s1-s2 plane that lies between the lines s2=ˉs2+R21(s1−λ1) and s2=λ2+R21(s1−ˉs1) and above or to the left of Γ−, i.e.,
Ωλ={(s1,s2):s1>ˉs1ors2>ˉs2,andV(λ1,ˉs2)≤V(s1,s2)≤V(ˉs1,λ2)}. | (27) |
For every (s1,s2)∈Ωλ, we have min{f1(s1),f2(s2)}>D, and so growth of x is always positive in this region.
Corollary 4.15. Assume (sin1,sin2)∈Ω1 and μ(r)>0. If ˉs1>λ1 and ˜s2≥λ2, then any solution (s1(t),s2(t),x(t)) with (s1(0),s2(0))∈Ωλ and x(0)>0 converges to the unique periodic orbit of (1).
Proof. First note that, since ˜s2≥λ2, the line through (ˉs+1,ˉs+2) is in Ωλ, and Ωλ∪Ω1A is connected. For any (s01,s02)∈Ωλ, we have I(s01,s02)>0 and (g∘π−)(s01,s02)∈Ωλ∪Ω1A. By (26), X(s01,s02)<0<x(0), and by Theorem 4.14, (s1(t),s2(t),x(t)) converges to the periodic orbit.
Example 2. Consider (1) with the Monod functional responses fi(si)=misiai+si, i=1,2, and parameters (m1,m2,a1,a2)=(2,2,1.4,0.6), (Y1,Y2,D)=(2,0.7,0.5), and (ˉs1,ˉs2,sin1,sin2,r)=(0.7,0.6,1,1,0.4). Then,
(ˉs+1,ˉs+2)=(0.82,0.76),˜s2≈0.42,ˆs2≈0.14,μ(r)≈0.04,V−≈−0.19. |
The equation fi(si)=D, i=1,2, yields λi=aiDD−mi, which gives λ1≈0.4677 and λ2=0.2. Since λ1<ˉs1 and λ2<˜s2, Ωλ∪Ω1A is a connected set, and the hypotheses in Corollary 4.15 are satisfied.
We take the initial value (s01,s02)=(0.6,0.7). Then, V(s01,s02)=−0.59. A direct calculation gives V(λ1,ˉs2)≈−0.79, so that V(λ1,ˉs2)<V(s01,s02)<0. This implies that (s01,s02)∈Ωλ. By Corollary 4.15, the fermentation succeeds for every initial value x(0)>0. An illustration is shown in Figure 5.
Remark 4.16. If λ2>˜s2 then the set Ω1A∪Ωλ is not connected. Then for some (s01,s02)∈Ωλ, (g∘π−)(s01,s02) is in the gap between Ωλ and Ω1A. We are unable to rule out the possibility that the net growth in this gap is negative, and so it is conceivable that I((g∘π−)(s01,s02))<0. In particular, we can choose (s01,s02)∈Ωλ such that I(s01,s02)<−I((g∘π−)(s01,s02)). Therefore,
X(s01,s02)≥−I(s01,s02)−(1−r)−1I((g∘π−)(s01,s02))>0. | (28) |
Therefore, for some positive initial concentrations of biomass, the reactor will fail, even though the initial conditions are in Ωλ.
In this section we regard r as a variable in the interval (r∗,1), where r∗ is the number given in Theorem 4.6.
For each r∈(r∗,1), there is a periodic orbit. In each period, there is exactly one impulse. As shown in the proof of Theorem 4.6, the left and right limits at an impulse are, respectively,
(ˉs1,ˆs2,x−(r))and(ˉs+1(r),ˆs+2(r),x+(r)), |
where
x−(r)=1rμ(r)andx+(r)=1−rrμ(r). | (29) |
The trajectory of the periodic orbit can be parametrized by
s2=ˆs2+R21(s1−ˉs1) |
and
x=X(s1;r)=x−(r)−Y1∫s1ˉs11−Dmin{f1(v),f2(ˆs2+R21(v−ˉs1))}dv | (30) |
=x+(r)+Y1∫ˉs+1s11−Dmin{f1(v),f2(ˆs2+R21(v−ˉs1))}dv | (31) |
with s1∈(ˉs1,ˉs+1) and ˉs+1=ˉs1+r(sin1−ˉs1).
Denote the minimal period of the periodic orbit by T(r). Then
T(r)=Y1∫ˉs+1(r)ˉs11min{f1(v),f2(ˆs2+R21(v−ˉs1))}X(v;r)dv. | (32) |
In the long run, the average amount of output divided by the total volume is
Q(r)=rT(r). |
Maximizing Q(r) for r∈(r∗,1) is equivalent to maximizing the output.
Lemma 5.1. The minimal period T(r), r∈(r∗,1) of the periodic orbit of (1) satisfies limr→1−T(r)=∞. Also limr→r∗T(r)=∞ if r∗>0.
Proof. As r→1, we have x+(r)→0. By (32) and (31), T(r)→∞.
If r∗>0, then by (29), x−(r)→0 as r→r∗. Along the periodic orbit, x≥x−(r)>0 and
min{f1(s1),f2(s2)}≥min{f1(ˉs1),f2(ˉs2)}>0. |
By (32) we conclude that T(r)→∞ as r→r∗.
Proposition 5.2. The function Q(r)=r/T(r), r∈(r∗,1), satisfies limr→1Q(r)=0. If r∗>0, then limr→r∗Q(r)=0. If r∗=0, ˉs1≥λ1, and ˆs2≥λ2, then limr→r∗Q(r)=min{f1(ˉs1),f2(ˆs2)}−D≥0.
Proof. By Lemma 5.1, limr→1Q(r)=1/(limr→1T(r))=0.
If r∗>0, then, by Lemma 5.1, limr→r∗Q(r)=r∗/(limr→r∗T(r))=0.
Next we assume r∗=0. By (32)
Q(r)=r/(Y1∫ˉs+1(r)ˉs11min{f1(v),f2(ˆs2+R21(v−ˉs1))}X(v;r)dv), | (33) |
where X is defined by (30). Since X(ˉs1;r)=x−(r)=μ(r)r, using L'Hôpital's rule and the definition of μ(r) in (16),
limr→0X(ˉs1;r)=Y1dˉs+1(r)dr|r=0(1−Dmin{f1(ˉs1),f2(ˆs2)}). |
Since ˉs+1=ˉs1+r(sin1−ˉs1), we have
dˉs+1(r)dr|r=0=sin1−ˉs1, | (34) |
and it follows that
limr→0X(ˉs1;r)=Y1(sin1−ˉs1)(1−Dmin{f1(ˉs1),f2(ˆs2)}). | (35) |
Furthermore, we have
T′(r)=Y1(sin1−ˉs1)min{f1(ˉs+1),f2(ˆs2+rR21(sin1−ˉs1))}X(ˉs+1;r)−Y1r2∫ˉs+1(r)ˉs1μ′(r)r−μ(r)min{f1(v),f2(ˆs2+R21(v−ˉs1))}X(v;r)2dv. | (36) |
Since ˉs+1(0)=ˉs and X(ˉs1;0)≠0, by (35) and (36) we obtain
T′(0)=Y1(sin1−ˉs1)min{f1(ˉs1),f2(ˆs2)}X(ˉs1;0)=1min{f1(ˉs1),f2(ˆs2)}−D. |
From the expression Q(r)=r/T(r), using L'Hôpital's rule, we conclude that limr→0Q(r)=1/T′(0)=min{f1(ˉs1),f2(ˆs2)}−D.
Assume r∗>0. Since limr→r∗Q(r)=limr→1Q(r)=0, by Proposition 5.2, Q(r) attains its maximum at some value of r in (r∗,1). Unfortunately, the analytical expression of the derivative ddrQ(r) is too complicated for finding a critical point of Q(r). We have only obtained the maximum value using a numerical simulation. An illustration of the maximal value of Q(r) is given in Figure 6.
We have modeled the self-cycling fermentation process assuming that there are two essential resources s1 and s2 that are growth limiting for a population of microorganisms, x, using a system of impulsive differential equations with state-dependent impulses. Assuming that the process is used for an application such as water purification, where the resources s1 and s2 are the pollutants, we assume that the threshold for emptying and refilling a fraction of the contents of the fermentor, resulting in the release of treated water, occurs when the concentrations of both pollutants reach an acceptable concentration set by some governmental agency. We called these thresholds, s1≤ˉs1 and s2≤ˉs2. We consider the process successful if once initiated, it proceeds indefinitely without a need for any subsequent interventions by the operator.
By solving the associated ODE system for s2 in terms of s1, we show that solutions, when projected onto the s1-s2 plane, are lines with slope given by the ratio of the growth yield constants. In order to derive necessary conditions for successful operation of the fermentor, we first divide the s1-s2 plane into two regions: Ω0 and Ω1. The model predicts that solutions of the associated system of ODEs with initial conditions in Ω0 approach the axes without ever reaching the thresholds for emptying and refilling and the reactor fails, independent of the initial concentration of microorganisms. Solutions of the associated system of ODEs with initial conditions in Ω1 have the potential to reach the threshold for emptying and refilling, but in this case, successful operation can also depend on the initial concentration of the population of microorganisms.
In most cases, at startup the input concentration of the pollutant would be the concentration of the pollutant in the environment, which we are assuming is constant, i.e., (s1(0),s2(0))=(sin1,sin2). If, for any solution starting at these input concentrations of the resources, (sin1,sin2), and positive concentration of biomass, x(0)>0, the threshold for emptying and refilling, s1≤ˉs1 and s2≤ˉs2, is reached with net positive growth of the biomass, the model analysis predicts that we can choose an emptying/refilling fraction, r, so that the system cycles indefinitely. In this case the solution approaches a periodic solution with one impulse per period.
If the system has a periodic solution, the (s1,s2) components of the periodic orbit lie along the line with slope given by the ratio of the growth yield constants joining (sin1,sin2) and the point in the s1-s2 plane where both thresholds are reached. The net change in the biomass on the periodic orbit, that we denote μ(r), must then also be positive.
For other initial conditions in Ω1, in order for the process to operate successfully, it is not enough that μ(r)>0. There is also a minimum concentration of biomass, X, that depends on the initial concentration of the resources, that is required for the reactor to be successful. If the initial concentration of biomass is larger than X, then our analysis predicts that the reactor will cycle indefinitely and solutions will converge to the periodic orbit. If the initial concentration of biomass is less than or equal to X, then the reactor will cycle a finite number of times and then fail. If there is no periodic orbit, then the reactor will either cycle a finite number of times and then fail, or will cycle indefinitely, but the time between cycles will approach infinity.
Besides depending on the initial concentration of the resources at start up, the minimum concentration of biomass at startup required for successful operation depends on the emptying/refilling fraction in an interesting way. The closer r is to one, the smaller the number of impulses that are required for solutions to get to the periodic orbit. However, the time spent in a region of negative growth could be larger, and so X would be larger. The closer r is to zero, results in less time spent in regions with negative growth, but more impulses are then required to get close to the periodic orbit. Each impulse removes biomass from the reactor, and so X would also increase. This implies that there is an optimal value of r for which the reactor has the best potential for success. The values of the growth yield constants, Y1 and Y2, also play a role in the size of X. If their ratio is held constant, but each value is scaled by a constant c>0, then X is also scaled by the same constant c. Knowing this is important when selecting the population of microorganisms to use in the process.
If the choice of potential microorganisms for use in the process is restricted, then it might be easier to treat more highly polluted water than less polluted water, provided the microorganisms are not inhibited at high concentrations of the pollutant. We have shown that for successful operation, it is necessary that an r exists such that μ(r)>0. One way to increase μ(r) without changing anything else is to increase sin1 and sin2 in such a way that it still lies on the same line as before. It is also important to choose a population of microorganisms so that (sin1, sin2) lies in Ω1. It might only be possible to do this by increasing the concentration of one of the pollutants. However, another possibility might be to pre-process the input with a different population of microorganisms that moves (sin1, sin2) into an acceptable position so that a second population can then treat the water effectively.
We also make what might appear to be other surprising observations. Although the break-even concentrations play a role, it is not necessary for both break-even concentrations to be below their respective thresholds for emptying and refilling for the process to be successful (see Figure 2). Also, the process can still fail when both break-even concentrations are below their respective thresholds (see Figure 2).
For growth on a single, non-inhibitory, limiting resource in the self-cycling fermentation process, it has been shown that when the system has a periodic orbit, every solution either converges to the periodic orbit, or converges to an equilibrium without a single impulse [19]. If the resource is inhibitory at high concentrations, it has been shown that solutions may also converge to an equilibrium after a single impulse, but if there are at least two impulses then the solution is destined to converge to the periodic orbit [7]. In contrast, if there are two limiting essential resources, we have shown that there may be many impulses before the system converges to an equilibrium, even when the system has a periodic orbit. The example in Figure 3 demonstrates failure after two impulses.
An important issue when setting up the self cycling fermentation process is the choice of the emptying and refilling fraction, r. In the application we considered we were interested in optimizing the total amount of output. In the example, shown in section 5, we demonstrated that the optimal value of the emptying/refilling fraction is r≈0.64. This result is consistent with what was shown in the single resource cases [7,19], where the optimal choice of r was also not 1/2. Another reason for implementing a self-cycling fermentation process instead of a continuous input process is to maximize the concentration of some microorganism in the output over some time period. For example, one recent 'proof of concept' study [23] investigated using the self-cycling fermentation process to improve the production of cellulosic ethanol production. In their investigation, and many other applications of self-cycling fermentation the emptying/refilling fraction r is set to one half. While this is convenient for experiments and measurements, our results indicate that this might not be the optimal choice of r.
All authors declare no conflicts of interest in this paper.
[1] | D. D. Bainov and P. S. Simeonov, Systems with Impulse Effect: Stability, Theory and Applications, Ellis Horwood Series: Mathematics and its Applications, Ellis Horwood Ltd., Chichester; Halsted Press [John Wiley & Sons, Inc.], New York, 1989. |
[2] | D. D. Bainov and P. S. Simeonov, Impulsive Differential Equations: Periodic Solutions and Applications, vol. 66 of Pitman Monographs and Surveys in Pure and Applied Mathematics, Longman Scientific & Technical, Harlow; copublished in the United States with John Wiley & Sons, Inc., New York, 1993. |
[3] | D. D. Bainov and P. S. Simeonov, Impulsive Differential Equations: Asymptotic Properties of the Solutions, vol. 28 of Series on Advances in Mathematics for Applied Sciences, World Scientific Publishing Co., Inc., River Edge, NJ, 1995, Translated from the Bulgarian manuscript by V. Covachev [V. Khr. Kovachev]. |
[4] | G. Butler and G. S. K. Wolkowicz, Exploitative competition in a chemostat for two complementary, and possibly inhibitory, resources, Math. Biosci., 83 (1987), 1–48. |
[5] | A. Cinar, S. J. Parulekar, C. Ündey and G. Birol, Batch Fermentors: Modeling, Monitoring, and Control, Marcel Dekker, New York, 2003. |
[6] | F. Córdova-Lepe, R. D. Valle and G. Robledo, Stability analysis of a self-cycling fermentation model with state-dependent impulses, Math. Methods Appl. Sci., 37 (2014), 1460–1475. |
[7] | G. Fan and G. S. K. Wolkowicz, Analysis of a model of nutrient driven self-cycling fermentation allowing unimodal response functions, Discrete Contin. Dyn. Syst., 8 (2007), 801–831. |
[8] | J. P. Grover, Resource Competition, Population and community biology series 19, Chapman and Hall, New York, 1997. |
[9] | A. Halanay and D. Wexler, Teoria Calitativa a Sistemelor du Impulsuri, Ed. Acad. RSR, Bucharedt, 1968. |
[10] | S. Hughes and D. Cooper, Biodegradation of phenol using the self-cycling fermentation process, Biotechnol. Bioeng., 51 (1996), 112–119. |
[11] | J. A. León and D. B. Tumpson, Competition between two species for two complementary or substitutable resources, J. Theor. Biol., 50 (1975), 185–2012. |
[12] | B. Li, G. S. K. Wolkowicz and Y. Kuang, Global asymptotic behavior of a chemostat model with two perfectly complementary resources with distributed delay, SIAM J. Appl. Math., 60 (2000), 2058–2086. |
[13] | MATLAB, version 9.3.0 (R2017b), The MathWorks Inc., Natick, Massachusetts, 2017. |
[14] | L. Perko, Differential Equations and Dynamical Systems, Springer-Verlag, 1996. |
[15] | A. Samoilenko and N. Perestyuk, Impulsive Differential Equations, World Scientific, Singapore, 1995. |
[16] | B. E. Sarkas and D. G. Cooper, Biodegradation of aromatic compounds in a self-cycling fermenter, Can. J. Chem. Eng., 72 (1994), 874–880. |
[17] | D. Sauvageau, Z. Storms and D. G. Cooper, Sychronized populations of Escherichia coli using simplified self-cycling fermentation, J. Biotechnol., 149 (2010), 67–73. |
[18] | H. L. Smith and P. Waltman, The Theory of the Chemostat: Dynamics of Microbial Competition., vol. 13 of Cambridge Studies in Mathematical Biology, Cambridge University Press, Cambridge, 1995. |
[19] | R. J. Smith and G. S. K. Wolkowicz, Analysis of a model of the nutrient driven self-cycling fermentation process, Dyn. Contin. Discrete Impuls. Syst. Ser. B Appl. Algorithms, 11 (2004), 239–265. |
[20] | Z. J. Storms, T. Brown, D. Sauvageau and D. G. Cooper, Self-cycling operation increases productivity of recombinent protein in Escherichia Coli, Biotechnol. Bioeng., 109 (2012), 2262–2270. |
[21] | D. Tilman, Resource Competition and Community Structure, Princeton University Press, New Jersey, 1982. |
[22] | J. Von Liebig, Die Organische Chemie in Ihrer Anwendung auf Agrikultur und Physiologie, Friedrich Vieweg, Braunschweig, 1840. |
[23] | J. Wang, M. Chae, D. Sauvageau and D. C. Bressler, Improving ethanol productivity through selfcycling fermentation of yeast: a proof of concept, Biotechnol. Biofuels, 10 (2017), 193. |
[24] | J. Wu, H. Nie and G. S. K. Wolkowicz, A mathematical model of competition for two essential resources in the unstirred chemostat, SIAM J. Appl. Math., 65 (2004), 209–229. |
1. | Tyler Meadows, Gail S.K. Wolkowicz, Growth on multiple interactive-essential resources in a self-cycling fermentor: An impulsive differential equations approach, 2020, 56, 14681218, 103157, 10.1016/j.nonrwa.2020.103157 | |
2. | Stacey R. Smith, Tyler Meadows, Gail S.K. Wolkowicz, Competition in the nutrient-driven self-cycling fermentation process, 2024, 54, 1751570X, 101519, 10.1016/j.nahs.2024.101519 |