Fractional Lotka-Volterra model with time-delay and delayed controller for a bioreactor

In this subsection a fractional Lotka-Volterra model with time-dealy for a bioreactor based on the Lotka–Volterra equations is introduced. One reason for choosing this well known predator-prey mathematical model is the appearance of periodic solutions, a phenomena observed in the experimental data of a bioreactor. Evidently, the Lotka-Volterra equations alone are not able to capture the rich variety of dynamics involved in a bioreactor. For this purpose, two context specific modifications are proposed: By following an approach gaining prominence in dealing with these complexities, first fractional exponents in the variables have been introduced [36, 37], and secondly, a time-delay term has been added in one of the differential equations [7]. In doing so it turned out the the mathematical model adequately fits the data provided.

For modeling a bioreactor, the following fractional Lotka-Volterra equations with time-delay and delayed control is introduced:

where $a,\ b,\ c,\ d,\ e\in\mathbb{R}^{+}$ are constant positive parameters of the bioreactor, $\alpha,\,\beta\in\mathbb{R}^{+}$ are positive fractional-exponents, $s_{0}\in\mathbb{R}^{+}$ is the positive initial supply of substrate, $\tau\in\mathbb{R}^{+}$ is a non-negative constant delay, $s(t)$ is the amount of substrate present in the bioreactor at time $t$, $x(t)$ is the amount of such produced biomass at time $t$, and $u(t)\in[0,1]$ is the input signal to manipulate the inflow of substrate into the bioreactor.

When setting $x(t)=0=u(t)$, the first equation in the model (1) reduces to $\dot{s}(t)=-as(t)^{\beta}$. That is to say, in the absence of biomass, the rate of change of the substrate is proportional to a scaling of the substrate by an exponent $\beta$. The negative sign reflects that the substrate in the bioreactor actually decays in the absence of biomass, in contrast to a classical Lotka-Volterra approach where intrinsic rate of prey population increase. It has been observed experimentally that the value of $\beta$ satisfy the inequality $0<\beta<1$, in particular the substrat, in absence of biomass, is subject to a non-exponential decay.

On the other hand, for $s(t)=0$, the second equation in the model (1) becomes $\dot{x}(t)=-c{x(t)}^{\alpha}+ex(t)$. That is to say, in the absence of substrate, the biomass follows a decay proposed in the Norton-Simons-Massagué (NSM) model. While the the NSM model actually describes the growth of tumors, the sign of the terms on the right hand side of the NSM model have been reversed to describe instead the decay of an organism under energy constrains. This approach has been successfully taken in very recent works on the analysis of cancer cells [37] and epidemiological models [36]. Here lies another of the proposed modification of the classical Lotka-Volterra model, where the rate of change of predators in absence of prey is assumed to be proportional to their population.

The Lotka-Volterra equations are complemented with the terms containing the scaled products which represent the rates of decay and growth provided by the biochemical iterations between substrate and biomass. That is to say, the terms $bx(t)^{\alpha}s(t)^{\beta}$ and $dx(t)^{\alpha}s(t)^{\beta}$ replace $bx(t)s(t)$ and $dx(t)s(t)$ used in the original model, respectively. The populations $x(t)$ and $s(t)$ are considered to be scaled by the same exponents $\alpha$ and $\beta$ described above.

Finally, the delay in the variable $x(t)$ in the first equation of the system (1) is justified as follows: the effect of the biomass in the rate of consumption of the substrate is affected by a dead time produced either by a gestation time-delay or by a time delayed measurement. Mathematically, this delay will turn out to produce the possibility for a stable equilibrium point of the fractional Lotka-Volterra model (1), in sharp contrast with the non-stability of equilibrium points of a classical Lotka-Volterra model. Thus the introduction of this time-delay favors the modeling of our bioreactor, where stability of fixed points is observed experimentally. Furthermore, biomass production $x(t)$ can only be manipulated indirectly by supplying substrate $s(t)$, so that the only controllable variable is the substrate $s(t)$ through a control action $u(t)\in[0,1]$ that regulates the flow of the substrate through the bioreactor. Here $u(t)=0$ represents zero flow (valve fully closed) and $u(t)=1$ represents total flow (valve fully open), see Figure 1. Therefore, the first equation in (1) has control action $u(t)$ while the second equation has no control action.

A summary of the parameters in the fractional Lotka-Volterra model (1) is shown next.

• •

  The parameter $a$ is the intrinsic rate of the amount of scaled substrate [Units: 1/h].

• •

  The parameter $\beta$ is the fractional-exponent scaling the amount of substrate [Units: Dimensionless].

• •

  The parameter $b$ is the rate of consumption of the scaled substrate by the scaled biomass [Units: 1/h].

• •

  The parameter $\alpha$ is the fractional-exponent of growth of the biomass [Units: Dimensionless].

• •

  The parameter $c\not=0$ denotes the net rate of growth the scaled biomass in response to the size of the scaled substrate [Units: 1/h].

• •

  The parameter $d$ measures the output of the scaled biomass from the bioreactor [Units: 1/h].

• •

  The parameter $e$ quantifies anabolism (growth rate) of the biomass.

To obtain a parametric identification of the model proposed in (1), we begin by assuming that the input signal $u(t)=D$ is constant and positive. This amounts to open the valve for the substrate to a constant input flow. Our purpose is to achieve the tuning of the parameters so that our model resembles the experimental data of a bioreactor. In addition, we will give conditions for obtaining limit cycles and Hopf bifurcations using the time-delay as a bifurcation parameter.

Thus, the mathematical model (1) is rewritten as

We remark that in our setting, equation (3) does admit a positive solution such that $s^{\ast}>0$. For completeness of our analysis, next we provide some algebraic conditions for such solution to exist, namely, conditions for the left hand side of equation (3) to change of sign for a couple of values of $x$. Indeed, if $x=\left(\frac{d}{e}\right)^{\frac{1}{1-\alpha}}$ then $d-ex^{1-\alpha}=0$ and so the left hand side of equation (3) reduces to $Ds_{0}$ which is positive. Similarly if $x=\left(\frac{d-cs_{0}^{\beta}}{e}\right)^{\frac{1}{1-\alpha}}$ then $\frac{d-ex^{1-\alpha}}{c}=s_{0}^{\beta}$ and so the left hand side of equation (3) reduces to $(-a-bx^{\alpha})s_{0}^{\beta}$. Thus, a condition for the left hand side of equation (3) to be negative and hence for equation (3) to admit a nonzero solution is $-a-bx^{\alpha}<0$. In in addition one requires the inequality $x^{\ast}<\left(\frac{d}{e}\right)^{\frac{1}{1-\alpha}}$ one also guarantees that $s^{\ast}>0$. In what follows, we will assume that the equilibrium points $x^{\ast}$ and $s^{\ast}$ obtained from Proposition 1 are both positive.

The previous result provides conditions of the delay $\tau_{0}$ for which the quasi-polynomial (4) has roots on the imaginary axis $\lambda_{0}=\pm i\omega_{0}$. If these roots are dominant, then it is natural to determine when the roots move to the right half-plane under infinitesimal variations of this delay $\tau_{0}$. The above is know as direction of the crossing [45] and it is determined using the following result.

Note that the sign given in (7) is independent of $\tau_{0}$, therefore if the linearization of system (2) is stable for $\tau=0$, then it remains stable for all $\tau\in(0,\tau_{0})$, where $\tau_{0}$ is the first value of (6). Thus the following result is evident.

In this section, we carry out a stability analysis of the fractional mathematical model (1) in closed-loop with a delayed controller of the type (24). As well as tuning the controller gains.

Using Proposition 5, the linearization of model (1) is of the form (23), where

Next, we propose the following delayed controller to stabilize the model

where $k_{r}\in\mathbb{R}$ is the controller gain and $h\in\mathbb{R}^{+}$ is the controller delay. Note that in (9) the proposed controller is only acting on the biomass $x$ and not on the substrate $s$. We base our proposal in the fact that, during the experiments, we only have access to the measurement of $x$ and not of $s$. In addition, the time-delay in the controller is justified by the time-delay in between the recordings of the biomass output.

Let $A_{0},\ A_{1}\in\mathbb{R}^{2\times 2}$ and $B\in\mathbb{R}^{2}$ given in (8), then the quasi-polynomial of closed-loop systems (23) with (9) is

where

As mentioned before, for parametric identification, in this section we assume that $u(t)=D$ in the mathematical model (1). In addition, experimental data is taken from a biochemical process as reference points to identification, see Table 1. For the simulations, we consider initial conditions of $s_{0}=10$ g/L for the substrate and $x_{0}=0.1$ g/L for the biomass, as well as a constant opening of the substrate flow valve equal $D=0.15$ 1/h. These initial parameters are the same as those implemented in the biochemical process described above. To obtain the values of the other parameters nonlinear regression was used based on the Levenberg-Marquardt least squares minimization algorithm. In consequence, the values obtained from our model come closer to the experimental data, see Figure 2 and 3. They are

Whereby, $\kappa_{1}=0.37985$, $\kappa_{2}=0.02011$ and $\kappa_{3}=0.09791$. The parameter values obtained in the present study fall within the range of those reported in the literature, due to the different operating conditions used in each case, i.e., different carbon source, continuous or batch operation, temperature, pH, among others ([4, 27, 30, 52]).

It is recommended, therefore, to use ($\epsilon_{1}$) in lieu of correlation-based measures to provide a relative assessment of model performance. The statistics used absolute values rather than squared differences. Parameter efficiency ($\epsilon_{1}$) for all variables using the L-M method are: Biomass (0.850) and Substrate (0.834) and, Interpretation of correlation-based measures, 0.85 indicate that the model explains 85.0% of the variability in the observed data. Therefore, the proposed model was able to predict experimental data.

Now, using Proposition 1 and parameters given in (16) an equilibrium point of the mathematical model (2) is

Therefore, the linearization of (2) using (16) and (17) is

and by (4) its quasi-polynomial is

Thus, the polynomial (5) gives

Finally, using Proposition 5 and the solution $\omega_{0}=0.23876$ of the above polynomial, we get some candidate values of delay $\tau=\tau_{0}$ which may be critical values of $\tau_{0}$, where the mathematical model (2) has a Hopf bifurcation. Some calculated critical values are

These values were obtained using the equation (6) for $n=0,1,2,3,4$. To corroborate that the points given in (18) are critical points, in Figure 4 the phase diagrams of model (2) are presented. Note that for initial conditions $(s_{0}=1.85,\,x_{0}=4.9)$ near the equilibrium point (17), even more the system response (2) remains oscillating around the equilibrium point, thus verifying the existence of a limit cycle on $\tau_{0}=4.9608$.

While in Figure 5 the stability of the model (2) for $\tau=1,\dots,4$ is shown, using phase diagrams. Note that for the same initial condition $(s_{0}=1.85,\,x_{0}=4.9)$, but different values of $\tau$, the system response converges to the equilibrium point (17), corroborating the stability of (2) for all $\tau\in(0,\tau_{0})$, as postulated in the Proposition 4.

On the other hand, by Proposition 3 we have that (7) is

Thus, in Figure 6, the system response and a phase diagram of (2) are presented, when $\tau=7$. Clearly, the solution $(s(t),\,x(t))$ diverges from the equilibrium point (17), showing via simulation that the model (2) is unstable, when $\tau>\tau_{0}$ and the initial condition is nearby the equilibrium point and corroborating the provisions of the Proposition 4.

Now, we stabilize the mathematical model (1) using delayed controllers of the form (9).

Consider the parameters given in (16), $\tau=7$. By Proposition 5 and $x^{*}=4.77631$ given, the equilibrium point of the model (1) is

Thus, the linearization of (1) using (8) is

and the quasi-polynomial (3.3) is

Using Proposition 6, the $\sigma$-stability regions in the parametric plane $h-k_{r}$ are given in Figure 7 for $\sigma>0$. Notice that these regions are concentric and their size decreases as $\sigma$ increases. Furthermore, a collapse of these regions to the point $(7.38,0.031)$ marked with a red asterisk $({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\ast})$ can be observed when $\sigma=\sigma^{*}$, which implies that $\sigma^{*}=0.24$ is the maximum achievable exponential decay of the linear systems (20) or quasi-polynomial (21).