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

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


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

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

Preliminaries. Let us consider a general model with impulsive effects:
where (, ) and (, ) are indefinitely differentiable with respect to (, ) ∈ Ω ⊂ R 2 , and , , and  are linearly dependent on  and , that is,   ,   ,   ,   ,   , and   are constant.Denote Then  IMP is referred to as the control set, which contains all the states at which the control strategy is taken on and  and  describe the effects of the control strategy.When the state variables (, ) arrive at the control set  IMP , the impulsive control measures are taken; then the state variables (, ) jump from  IMP to the set  PHA .Let x() = ((), ()) be any solution of system (5).Let   (x, x 0 ) = x(), given the initial value x(0) = x 0 .Then, the (forward) orbit is the set of all values that this trajectory obtains (x, Let  ⊆ R 2 = (−∞, +∞) 2 be an arbitrary set, let  ∈ R 2 be an arbitrary point, and  > 0. Then the distance between  Discrete Dynamics in Nature and Society and  and the noncentral −neighborhood of  are denoted by Definition 1 (periodic orbit [36,37]).An orbit (x) of system ( 5) is said to be periodic if there exists a positive integer  ⩾ 1 such that x  = x 0 .Denote  = min{ | x  = x 0 }.Then the orbit (x) is said to be an order- periodic orbit.
Definition 2 (-close [38]).Two orbits (x) and (x) are close if there is a reparameterization of time (a smooth, monotonic function) Definition 3 (orbitally stable [36][37][38]).An orbit (x) is said to be orbitally stable if, for any  > 0, there is a neighborhood  of x so that, for all x in , (x) and (x) are -close, as illustrated in Figure 2.
Definition 5 (successor function [37]).Suppose that the control set  IMP and the phase set  PHA are both lines, as illustrated in Figure 4.The coordinate on the phase set  PHA is defined as follows: the origin is defined as the intersection point between  PHA and the -axis, denoted as   , and the coordinate of any point  on  PHA is defined as the distance between  and   .Let  3 denote the intersection point between the trajectory starting from

Main Results
It is visible that  IMP = {(, ) |  =  CL ,  ≥ 0} and Integrating both sides of (10) from  0 to  yields Now it is assumed that (( 0 ), ( 0 )) ∈  PHA ; that is, ( 0 ) =  UP .Let  1 be the first time for the trajectory to reach the control set  IMP ; that is, ( 1 ) =  CL .Then Denote  1 = ( 1 ).Due to the effects of the impulsive control, one has Thus Theorem 7. System (4) admits a unique order-1 periodic solution.
Proof.For the order-1 periodic solution  = (),  = (), by (11), one has Meanwhile, by ( 4), one has Substituting ( 19) into (20) yields Integrating both sides of ( 21) from  0 to  0 +  yields Besides the qualitative analysis of model (4) presented above, the stability of the controlled system is also an important issue to be considered.As far as the stability is concerned, most of the works in literature rely on the Analogue of Poincaré Criterion [36].Recently, a geometric method is used to verify the stability of the order-1 periodic solution for a general semicontinuous dynamic system [13].For the convenience of readers, a detailed overview is presented.

Stability Criterion Representation.
To test the stability of the order-1 periodic solution  = (),  = (), the behaviour of the orbit of system (5) starting from    () ∩  PHA for given small  > 0 should be characterized.Let  be a point on the order-1 periodic orbit with rectangular coordinates ((), ()), where 0 ≤  ≤  is the arc length between  and .To begin with, assume  = ; that is,  = 0.For a given point  ∈    () ∩  PHA , here assume that  locates above , where the rectangular coordinates of  are ((), ()).The trajectory starting from  will intersect the normal of  and the intersection point is denoted by   .When  increases, this trajectory intersects the normal of  and the control set  IMP at   and , respectively.Define If  Γ < 1, then the order-1 periodic solution  = (),  = () is orbitally asymptotically stable.
Theorem 10 (stability criterion [13]).The order-1 -periodic solution  = (),  = () of system (5) Denote Then which implies that the phase set  PHA is also a straight line with slope   PHA , where On the other hand, by the effect of the impulsive controls, one has x O s Figure 5: Illustration of the disturbance nearby the periodic orbit of model (5).
Let  denote the angle between the normal at  and the phase line and let  denote the angle between the normal at  and the impulse line.When  is sufficiently close to , one has (,   ) ≈ (, ) cos(), (,   ) ≈ (, ) cos(), which implies that Denote  1 ( 1 ) as the angle between  PHA ( IMP ) and the -axis and  2 ( 2 ) as the angle between the tangent at  () and the -axis.Then there are tan( To characterize the ratio (,   )/(,   ), it needs to introduce the curvilinear coordinates (, ); the description is as follows:  represents length of the arc starting from  and its increasing direction is consistent with that of time ;  represents the normal length and its positive direction is to the right side when traveling along the periodic orbit.
The relationship between the rectangular coordinates (, ) and the curvilinear coordinates (, ) on point  is  = () −   (),  = () +   (), where 0 and  0 denote the values of  and  at the order-1 periodic orbit x = ((), ()); that is,  0 = ((), ()) and  0 = ((), ()).It can be inferred that Since  and  are periodic functions of , then (, ) in ( 36) is also a periodic coefficient nonlinear equation with  as the variable and  as the undetermined function.By (35) it can be obtained that  = 0 is a solution of (36) corresponding to the order-1 periodic orbit.Since  and  are continuously differentiable, (, ) also has continuous partial derivatives with respect to , and (36) can be rewritten as In order to calculate    (, )| =0 , it is noted that (  ()) 2 + (  ()) 2 = 1, which implies that     +     = 0; that is,  0   +  0   = 0. Thus, Therefore, where () represents the curvature of the orthogonal trajectory of ( 5) at point .The first-order approximation of ( 37) is / = (), with the solution Note that which results in the fact that Combined with (34), this yields the conclusion.

Verifications and Biomass Productivity Optimisation
The microbial cultivation process for pulse dosage supply of substrates and removal of bioreactor medium has been   model parameters.Furthermore, from Figure 6(b), it can be concluded that the biomass concentration is controlled below a certain level through substrate regulation.Figure 7 illustrates the order-1 periodic solution with period  = 0.6 [h] = 36 [min].The phase diagrams for different initial biomass concentrations on  PHA are illustrated in Figure 8, from which the stability of the order-1 periodic solution can be observed.

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

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

Figure 1 :
Figure 1: Schematic view of the impulsive bioprocess.
analyzed theoretically.In order to complement the theoretical results, simulations of system (4)'s dynamics are presented firstly.The model parameters used in simulations are as follows: M = 0.5 [1/h], ] = 0.5,  ] = 1 [g/L],   = 16 [g/L],  = 0.5,  = 0.025,  = 0.8,  = 0.4, and   = 15% = 0.15.The low substrate control level is a decision parameter, which should be set according to an actual demand.Firstly, in order to verify the theoretical results,  CL is set to 10%   = 1.6 [g/L].Next, an optimal  * CL is determined by a biomass productivity minimization problem.4.1.Verification of the Main Results.The change of substrate concentration (), biomass concentration (), and phase diagram (, ) starting from ( 0 ,  0 ) = (12, 0.5) are shown in Figure6.It can be observed that the trajectories tend to be periodic, and this phenomenon is not dependent on the

Figure 8 :Figure 9 :
Figure 8: Verification of Theorem 7. The phase diagrams for different initial biomass concentrations on  PHA .

sFigure 10 :
Figure 10: The dependence of BP  on the substrate control level  CL for different   .

Figure 11 :
Figure 11: The dependence of BP  on   for different substrate control level  CL .

Figure 12 :Figure 13 :
Figure 12: The dependence of BP  on the model parameters  and  for the given  CL = 6 [g/L] and   = 0.1.

Figure 14 :
Figure 14: The dependence of BP  on the model parameters  for the given  CL = 6 [g/L],   = 0.1.
1 and the control set  IMP , and  4 denote the phase point of  3 after control measures.Let  2 denote the intersection point if the trajectory starting from  1 intersects the phase set  PHA .Then  4 is called as the successor point of  1 , and the mapping from  1 to  4 constructs the successor function.If  2 exists, the successor function is called type-II; that is,  II SOR ( 1 ) = ( 4 ,   ) − ( 1 ,   ); otherwise, it is called type-I; that is,  I SOR ( 1 ) = ( 4 ,   ) − ( 1 ,   ).Moreover, the successor function is continuous on  PHA .If there exists a point (  ,   ) ∈  PHA such that +       + [(1 +   )   −     ] +   + [(1 +   )   −     ] As illustrated in Figure 5, the slope of the control set is   IMP = −   /   (in case of    = 0,   IMP is equal to ∞).Due to the specific forms of the impulsive control, for any (  ,   ) ∈  IMP , if    ,    ,    , and    are constants, one has and tan( 2 ) =    .Thus 2 (,   )  (,   ) .