Numerical algorithms for solving the optimal control problem of simple bioreactors

Abstract. The modified nonlocal feedback controller is used to control the production of drugs in a simple bioreactor. This bioreactor is based on the enzymatic conversion of substrate into the required product. The dynamics of this device is described by a system of two nonstationary nonlinear diffusion–convection–reaction equations. The analysis of the influence of the convection transport is one the aims of this paper. The control loop is defined using the relation, which shows how the amount of the drug produced in the bioreactor and delivered into a human body depends on the substrate concentration specified on the external boundary of the bioreactor. The system of PDEs is solved by using the finite volume and finite difference methods, the control loop parameters are defined from the analysis of stationary linearized equations. The second aim of this paper is to solve the inverse problem and to determine optimal boundary conditions. These results enable us to estimate the potential accuracy of the proposed devices.


Introduction
In this paper, we continue the analysis of bioreactors based on the enzymatic conversion of substrate into the required product. In [3], the control loop was defined by using the relation how the amount of the drug produced in the bioreactor and delivered into a human body depends on the substrate concentration specified on the external boundary of the bioreactor. Only diffusion transport and nonlinear reaction processes were taken in account. The control loop parameters were defined from the analysis of stationary linearized equations. In this paper, we add the convection transport possibility. Thus the dynamics of the new device is described by a system of two nonstationary nonlinear diffusion-convection-reaction equations.
We note that the feedback control mechanism is used in many technological applications [11,12]. The recent developments of this technique for optimal control of processes in smart bioreactors is one of the interesting new theoretical and computational challenges [3,7,13].
The first aim of this paper is to investigate the influence of the convection transport to the efficiency of the bioreactor device. The second aim is to solve the inverse problem and to determine optimal boundary conditions for the control loop. These results enable us to estimate the accuracy of the proposed control loop technology.
We note that smart biochemical devices enable automatic adaptation of rates of produced drugs to the treatment procedures specified by medical doctors. Such devices are convenient, flexible and robust tools for patients.
The rest of this paper is organized as follows. In Section 2, a system of two nonstationary nonlinear diffusion-convection-reaction equations is formulated. This model describes the dynamics of the substrate S (prodrug) and the product P (drug). The three classical boundary conditions (two for P and one for S) are specified at the boundary of domain. The last fourth boundary condition should be determined from the objective function of the optimization problem, since the required flux of the product P is specified at the boundary x = 0. In order to define a well-posed boundary value problem, we formulate the inverse problem to find the equivalent boundary condition for the substrate function S(X, t). Here X defines the length of the bioreactor.
In Section 3, the algorithm of the proportional nonlocal controller is formulated. It is based on a modification of the algorithm proposed in [3]. The parameters of the control function are defined by solving the stationary (limit) system of equations when the nonlinear interaction term is linearized around some constant value of S.
The finite volume method is used to approximate the diffusion and convection processes in Section 4. The time derivatives and reaction terms are approximated by the symmetrical Euler method. The predictor-corrector method is applied to linearize the obtained system of nonlinear discrete equations.
The inverse problem to find optimal reconstructed boundary conditions S(X, t) is solved in Section 5. A simple derivative-free heuristic algorithm is proposed to solve the given problem.
In Section 6, results of computational experiments are presented and analysed. Final conclusions are presented in Section 7.

Mathematical models
In this paper, we consider an one-dimensional mathematical model used to simulate dynamics of simple bioreactors [5]. It is based on a system of two equations where D = {0 < x < X, 0 < t T }, S(x, t) and P (x, t) are real valued functions, S defines the concentration of the substrate of the enzyme and P is the concentration of a product. During the enzymatic reaction, the substrate S (which is a prodrug material) is converted into the active drug P . Such a technology is an example of smart technologies targeted to produce the drug on demand. The reaction conversation is described as the most simple Michaelis-Menten process. It is interesting also to consider more complicated bioreaction processes [5].
Diffusion coefficients D S , D P and convection coefficient α are assumed to be constants. The enzyme is uniformly distributed in the reactor, and the substrate S and the product of the reaction P are transported by the diffusion and convection processes inside the reactor.
The well-written review of mathematical models on nonlinear bioreactions is given in [2,10], a practical user guide on such models is presented in [5]. Similar diffusionreaction models are considered in [3,7].
In order to define a full mathematical model, we formulate initial conditions and three boundary conditions: The last condition specifies the flux of P at the boundary x = 0: where Q(t) defines the flux of the drug prescribed by a medical doctor in accordance with the therapeutic protocol (the objective function in the optimization task). Equation (4) is written by using the full flux of the product P (x, t) at the boundary x = 0. This form is consistent with the conservation law defining the physical process, and it should be used for the approximation in order to guarantee the discrete conservation property. Still in the definition of the full mathematical model, condition (4) can be written in a more simple form It is clear that such a combination of boundary conditions (3)-(4) is not used in any classical well-posed boundary value problem. In order to apply such bioreactors in real life applications, it is proposed in [3] to find the equivalent boundary condition for the substrate function S(X, t) = s(t), where s(t) is unknown function. When this function is specified (or obtained by solving some auxiliary problem), then we get a well-posed nonstationary boundary value problem to find S.
In general, the inverse problems are ill posed [1]. As a more flexible approach, we consider the variational problem [14]: find s(t) such that where W defines a feasible set of boundary conditions, and P s defines a product function when s(t) is used as a boundary condition in (5). In the following sections, we will solve this inverse problem for some typical cases of the objective function Q.

The feedback control loop
In this section, we briefly define a modification of the feedback control loop technology from [3], where it was used to achieve the desired regime of drug production for the diffusion-reaction model. We generalize it for the diffusion-convection-reaction models.
The main idea is instead of solving directly the inverse problem (6) to define a simple control algorithm based on some efficiently manipulated variable. In this way, an equivalent well-posed boundary value problem is formulated. It should guarantee that a close approximation of the required flux of drugs at the boundary x = 0 is produced. Thus, in general, we are interested to develop a dynamic control system based on proportional feedback controllers.
Our approach is different from a classical proportional-integral-derivative controller (PID controller) technique, which was used in paper [7]. It is well known that the main challenge in application of PID controllers deals with the selection of optimal values of coefficients K p , K i , K d . This task requires to do big scale computational experiments and the stability analysis of the obtained controllers is very nontrivial.
Our approach to construct the proportional controller is based on the idea to control a so-called steady-state error [3]. The asymptotic analysis of stationary (limit) system of equations is done and the nonlinear interaction term is linearized around some constant value of S. Due to the maximum principle, it is recommended to linearize this term around a zero value of S.
The following system of two linear differential equations for functions S(x) and P (x) is considered: The solution of problem (7) is given by Substituting it into (8) and integrating, we get the full flux of P at x = 0: The control algorithm based on the boundary value of substrate S It is proposed in [3] to use relation (9) in order to define the proportional controller algorithm. The supply of the substrate into the bioreactor is specified as where Q(t) defines the required flux of the product at the boundary x = 0. The coefficient 1/µ(α) defines the so-called control gain parameter [12]. The dependence of 1/µ(α) on the convection velocity α is presented in Fig. 1. The remaining parameters of the mathematical model are defined in Section 6.
Remark 1. It follows from the presented results that, due to the convection process, the required amount of substrate supply is decreasing for the increased velocity α, thus the convection reduces the costs of the bioreactor.
The smart bioreactors have a possibility to perform the electrochemical monitoring of the enzymatic reaction. It is assumed that the concentration of the produced drug flux Q R (t) can be measured. Then we can use the feedback control loop to regulate the total amount of the produced product. Following [3], we add the correction term in the definition of the modified objective function Q(t): Then Q(t) is used in the definition of the control boundary condition (10). In this algorithm, the surplus/deficit of the produced drug is distributed uniformly over the remaining time interval, and the surplus of the product is compensated dynamically. In this section, we consider the discrete approximation of problem (1)-(4). Let Ω t be a t-grid Ω t = t n : t n = t n−1 + τ, n = 1, . . . , where τ is the discretization step size. Also we introduce a uniform spatial grid We consider numerical approximation U n j to the exact solution values of function U (x j , t n ) at the grid points (x j , t n ).
For functions defined on the grid Ω x × Ω t , we introduce the backward difference quotient and the averaging operator with respect to t and three difference operators with respect to x: We use the finite volume method to approximate the diffusion and convection processes in each space control volume [x j−1/2 , x j+1/2 ]. The obtained fluxes are approximated by the standard finite difference operators. This method leads to the conservative discretization in space. The finite difference method is used to approximate the time dependent operators by the symmetrical Euler scheme. Thus the differential problem (1)-(4) is approximated by the discrete scheme where Q Rh defines a measured value of the product The convection process can be approximated by using the monotone upwind approximation The proposed discrete model includes the dynamic control condition (13). The nonlinear boundary value problem (12) is linearized by using the predictor-corrector technique. The approximation error is of order two with respect to time and space if the central difference scheme is used to approximate the convection process, i.e. it is bounded by C(τ 2 + h 2 ). If the convection is approximated by the upwind scheme, then the approximation error is bounded by C(τ 2 + h).
For a fixed proportional type control algorithm, the boundary condition S(X, t) is defined as a scaled function Q(t). Thus we have a classical set of boundary conditions, and this modified initial-boundary value problem is well posed. In this case, the known theoretical results [4,6] can be applied to prove that the solution of the discrete scheme (12)- (14) converges in the L ∞ -norm with the order equal to the accuracy of approximation. It is obvious that we are getting only a control based solution of the modified problem. The main aim is to select efficient control strategies enabling us to approximate accurately the objective function Q(t). (13) is nonlocal. Since the nonlocal and nonlinear terms are approximated on (n − 1)th level, the standard factorization algorithm is used to solve the obtained systems of linear equations. More details on discrete approximations of problems with nonlocal boundary conditions are given in [8].

Inverse reconstruction of the boundary condition
In this section, we solve the discrete version of the variational problem (6). Our aim is to do an inverse reconstruction of the boundary condition (5). Let ω t be a set of time points Thus we formulate the following discrete variational problem: Next, we define the algorithm to find an approximate solution of problem (16). In order to get an initial approximation of the solution, we use an appriori information defined by the control formula (10). We also mention general derivative-free DIRECT type algorithms, which are widely used to solve global optimization problems for a black box type applications [9].
The proposed optimization algorithm consists of two steps.
Step 1: The local optimization step. Let us assume that, for each j, three different values of s j are selected: Then applying the full search technique we find the optimal set of parameters . . , M. This step requires to solve the discrete problem (12)- (14) for 3 M different boundary conditions S M . Remark 3. All 3 M problems are independent, thus they are solved by using a parallel master and slave version of the algorithm. For selected sizes of discrete meshes, the computation time of one job is much larger than the communication time between the master and slave processes. The tasks are distributed dynamically, and a deterministic predictor-corrector algorithm is used to solve nonlinear systems of equations. Thus the efficiency of the parallel algorithm is close to one.
Step 2: A set of parameters is updated. After finishing Step 1, the feasible set of parameters (17) is updated in the following way. If k j = 2, then If k j = 1, then If k j = 3, then The iteration stopping condition is defined as s j3 − s j1 ε, j = 1, . . . , M.
Our aim is to test the robustness of the proposed control scheme in cases when the convection starts to be a dominating transport process. Also we compare this simple proportional control algorithm with more sophisticated control schemes, which include the ideas of time-delayed feedback control [12] or the inverse determination of the control function [14].
First, we remind the important fact, that the drug production process reacts with a fixed delay to the changes of the substrate boundary condition (5) [3]. A brief theoretical justification of this effect in the case of pure diffusion transport process is given in [3].
Here we present results of numerical experiments when the convection transport process is also included into the model.
First, we investigate the case when α = 1. In this example, the convection transport is still weak. In Fig. 2, the dynamics of product P molar flow rate at x = 0 is shown (a scale factor 10 −6 should be added) when the objective function Q is a stepwise function: where 1/µ(α) is the gain control coefficient. Such scaling of P function is used in all figures.
The three control algorithms are compared. Our first conclusion is that the simple proportional control algorithm (10) is quite robust and the produced flux of the product follows the requite regime of the objective function. Still we see two weaknesses of this  algorithm. First, it produces the amount of product smaller than required. To address this topic, the corrected control algorithm (11) is applied. This technique dynamically compensates the deficit of the produced drug. Second, the standard proportional control algorithm generates the flux of the product, which follows the dynamics of the objective function with some time delay. In order to show that the latter issue can be resolved within the framework of the mathematical model (1)-(4), we present results when the control boundary condition is obtained solving the inverse variational problem (6). The spline function S M (t) (15) is defined on M = 7 time points. It follows from the obtained results that state of the art modelling techniques can help to develop very accurate smart medical devices.
Computations of the optimization algorithm were performed using resources at the High Performance Computing Centre "HPC Saulėtekis" in Vilnius University Faculty of Physics. The master and slaves template was used to implement the given search algorithm.
Next, we investigate the dependence of the solutions of problem (1)-(4) on the convection velocity α. The corrected proportional algorithm (11) is used in all experiments. The results are presented in Fig. 3. It is clearly seen that, for larger values of α, the convection transport process starts to dominate, and the substrate (S) is moved to the left side of the device faster than the enzymatic conversion reaction is going on. Thus too much of the unused substrate is accumulated near the boundary x = 0. Still we want to note that, for larger values of α, the required amount of substrate is reduced in comparison with pure diffusion transport case (see Fig. 1).
We propose one more modification of the control algorithm, which is based on the well-known time-delayed feedback control method [12]. In our algorithm, we apply the forward looking control boundary condition S(X, t) = 1 µ(α) Q(t + τ ),  where τ defines the time delay in the changes of the product flow with respect to perturbations of the control boundary conditions. This parameter should be fitted to experimental results. The results of this modified control algorithm are shown in Fig. 3 (the brown line). Next, we use the proposed algorithms to reconstruct a flux function typical in realworld applications The corrected proportional control algorithm (11) is used to define the required boundary conditions. The results of experiments are presented in Fig. 4 for α = 0, 5, 10. Again, we see that the proposed control algorithm is robust for both the diffusion and convection dominated transport cases. The modified time-delayed feedback control scheme with τ = 0.1 reduces the delay effect.
It is important to note that the variational approach to define the optimal control boundary conditions gives a very good quality of the flux dynamics. For the convection parameter α = 10, the spline function S 9 (t) (15) is defined by using M = 9 time points. The results of computations are presented in Fig. 4 (see the violet line).

Conclusions
In this paper, we have modified the feedback control algorithm from [3] for a mathematical model, which includes also the convection transport process and describes the drug delivery system. The diffusion-convection-reaction system simulates the enzymecontaining bioreactor, and the prodrug is converted into an active drug during the reaction.
The model is simulated numerically, and the finite volume method is used to approximate the given nonstationary equations. The system of partial differential equations is approximated with the second-order accuracy in space and time. Results of computational experiments show that the proposed control algorithm is accurate and robust. Thus it can be used in medical practices.
The potential accuracy of the proposed control algorithms is analysed by comparing the obtained boundary conditions with the optimal boundary conditions. Such optimal conditions are computed by solving the inverse variational problem. The parallel masterslaves method is used to implement the search algorithm. The efficiency of the obtained parallel algorithm is close to one.