## 1 Introduction

*et al.*, 2017; Kok Kiong

*et al.*, 1999; Yordanova and Ichtev, 2017). Our main aim is to propose a new feedback control algorithm for advanced bioreactors by using nonlocal formulations of the control functionals. This method enables automatic adaptation of rates of produced drugs to the treatment procedures specified by medical doctors. Such a technology gives a very convenient, flexible and robust tool for patients. The analysis is based on virtual simulation of real physical processes and demonstrates a potential of virtual mathematical modelling technique in biomedicine applications.

*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 defines the flux of

*P*at $x=0$. Thus a nonclassical combination of boundary conditions is used. The inverse problem is formulated to find the equivalent boundary condition for the substrate function $S(X,t)$.

## 2 Mathematical Models

##### (1)

\[\begin{array}{l}\displaystyle \frac{\partial S}{\partial t}={D_{S}}\frac{{\partial ^{2}}S}{\partial {x^{2}}}-\frac{VS}{{K_{M}}+S},\hspace{1em}(x,t)\in D=\{0<x<X,\hspace{2.5pt}0<t\leqslant T\},\\ {} \displaystyle \frac{\partial P}{\partial t}={D_{P}}\frac{{\partial ^{2}}P}{\partial {x^{2}}}+\frac{VS}{{K_{M}}+S},\end{array}\]*t*and

*x*are time and space variables, $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. This type of bioreactors is interesting for medical applications, since the enzymatic reaction converts a substrate of the enzyme

*S*(which is a prodrug material) into the active drug

*P*. Such technology can be considered as a smart technology producing the drug on demand. In the presented model the reaction conversation is described as the most simple Michaelis-Menten process. An extended review of models on nonlinear reactions in bioreactors is given in Murray (2003), Čiegis and Bugajev (2012), a very good practical user guide on such models is presented in Hillen and Painter (2009). We note that this model is also considered in paper Ivanauskas

*et al.*(2017). The enzyme is uniformly distributed in the reactor, and the substrate S diffuses in the bioreactor with the diffusion coefficient ${D_{S}}$. The product of the reaction

*P*diffuses inside the bioreactor with the diffusion coefficient ${D_{P}}$.

##### (3)

\[\begin{array}{l}\displaystyle P(0,t)=0,\hspace{2em}{D_{P}}\frac{\partial P}{\partial x}(X,t)=0,\hspace{1em}t>0,\\ {} \displaystyle {D_{S}}\frac{\partial S}{\partial x}(0,t)=0.\end{array}\]*P*at the boundary $x=0$: where $Q(t)$ defines the flux of the drug prescribed by a doctor in accordance with the therapeutic protocol.

*S*we get a well-posed non-stationary boundary value problem.

*et al.*, 2012). A more flexible for applications mathematical model is obtained if for $s(t)$ we consider the variational problem (Tikhonov and Arsenin, 1977): find $s(t)$ such that

##### (6)

\[ \bigg|{D_{P}}\frac{\partial {P_{s}}}{\partial x}(0)-Q\bigg|=\underset{\tilde{s}\in W}{\min }\bigg|{D_{P}}\frac{\partial {P_{\tilde{s}}}}{\partial x}(0)-Q\bigg|,\]*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 next section we investigate the stability of the obtained inverse problem and show that it can be treated as a well-posed model with a bounded stiffness constant.

*S*, if the function $q(t)$ is given. For the system of equations (1)–(4) this function should be obtained by solving the inverse problem.

## 3 The Delayed Feedback Control Loop

*et al.*, 1999). Instead of solving directly the inverse problem for the boundary condition (5) (or (7)) we consider the approach based on the nonlocal delayed feedback control method. Our aim is to select some efficient manipulated variable and to formulate an equivalent well-posed boundary value problem in order to produce the required flux of drugs at the boundary $x=0$. Thus we are interested to develop a dynamic control system based on proportional delayed feedback controllers.

*et al.*, 2017). The authors have attempted to minimize the error over the drug production by adjustment of a control variable $S(X,t)$ to a new value determined by a weighted sum: where $e(t)=Q(t)-{Q_{R}}(t)$ is the difference between a desired amount of produced drugs $Q(t)$ and a measured process variable ${Q_{R}}(t)$. ${K_{p}}$, ${K_{i}}$ and ${K_{d}}$ denote non-negative coefficients for the proportional, integral, and derivative terms of the error. The selection of $S(X,d)$ as a control variable seems quite questionable in this model. The stability analysis of the proposed PID controller is not presented in Ivanauskas

*et al.*(2017) and optimal values of coefficients ${K_{p}}$, ${K_{i}}$ and ${K_{d}}$ are selected experimentally.

*P*at $x=0$ depends monotonically on the boundary value

*s*of the substrate concentration

*S*.

*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 $\tilde{S}(x)$ and $\tilde{P}(x)$ is considered:

### 3.1 The Control Algorithm Based on the Boundary Value of Substrate *S*

*t*we can apply one iteration of the boundary value error correction and also include the information on a change of $Q(t)$ values over time:

*τ*defines a time step. The obtained control algorithm can be considered as a representative of delayed feedback control algorithms (Pyragas, 2006; Novičenko, 2015). Such algorithms are often used in practical computations, but in order to guarantee the stability of the control technique the parameter

*μ*should be adapted to the behaviour of the system and it is not sufficient to consider the solution of a stationary nonlinear system.

##### (13)

\[ \tilde{Q}(t)=Q(t)+\bigg({\int _{0}^{t-\tau }}\hspace{-0.1667em}\hspace{-0.1667em}Q(s)\hspace{0.1667em}ds-{\int _{0}^{t-\tau }}\hspace{-0.1667em}\hspace{-0.1667em}\hspace{-0.1667em}{Q_{R}}(s)\hspace{0.1667em}ds\bigg)\Big/(T+{T_{0}}-t),\hspace{1em}t-\tau <t\leqslant T.\]*τ*defines a time step of the numerical integration algorithms. It should not be interpreted as a time delay of the system reaction to the changes of the boundary condition (5). In the following sections we investigate such a delay in more details.

### 3.2 The Control Algorithm Based on the Flux of Substrate *S*

*t*is defined as

### 3.3 Approximation of the Flux of *P*

*P*at point $x=0$:

*P*at $x=0$ can be approximated by the integral term:

##### (17)

\[ {D_{P}}\frac{\partial P}{\partial x}(0,t)\approx \frac{2{D_{P}}}{{b^{2}}-{a^{2}}}{\int _{a}^{b}}P(x,t)\hspace{0.1667em}dx.\]## 4 Finite Volume Scheme

*t*-grid

*τ*is the discretization step. Also we introduce a uniform spatial grid

*t*and two difference operators with respect to

*x*:

##### (18)

\[ {\partial _{\bar{t}}}{S_{j}^{n}}={D_{S}}{A^{h}}{S_{j}^{n-1/2}}-\frac{V{S_{j}^{n-1/2}}}{{K_{M}}+{S_{j}^{n-1/2}}},\hspace{1em}{x_{j}}\in {\Omega _{x}},\]##### (19)

\[\begin{array}{l}\displaystyle -{D_{S}}{\partial _{x}}{S_{1}^{n-1/2}}+\frac{h}{2}\bigg({\partial _{\bar{t}}}{S_{0}^{n}}+\frac{V{S_{0}^{n-1/2}}}{{K_{M}}+{S_{0}^{n-1/2}}}\bigg)=0,\\ {} \displaystyle {S_{J}^{n}}=\frac{1}{\mu }\bigg[Q({t^{n}})+\bigg({\int _{0}^{{t^{n-1}}}}Q(s)ds-{\sum \limits_{m=1}^{n-1}}{Q_{Rh}^{m}}\tau \bigg)\Big/\big(T+{T_{0}}-{t^{n-1}}\big)\bigg],\end{array}\]##### (20)

\[\begin{array}{l}\displaystyle {\partial _{\bar{t}}}{P_{j}^{n}}={D_{P}}{A^{h}}{P_{j}^{n-1/2}}+\frac{V{S_{j}^{n-1/2}}}{{K_{M}}+{S_{j}^{n-1/2}}},\hspace{1em}{x_{j}}\in {\Omega _{x}},\\ {} \displaystyle {P_{0}^{n}}=0,\hspace{1em}{D_{P}}{\partial _{x}}{P_{J}^{n-1/2}}+\frac{h}{2}\bigg({\partial _{\bar{t}}}{P_{J}^{n}}-\frac{V{S_{J}^{n-1/2}}}{{K_{M}}+{S_{J}^{n-1/2}}}\bigg)=0,\end{array}\]*et al.*(2016) for more details on discrete approximations of problems with nonlocal boundary conditions). The stability analysis of the dynamical process will be considered in the next section.

*h*. The error due to approximation of the flux by the integral formula is bounded by $Cb$. This error will make a small perturbation of the real flux of the produced drug, since the control technique depends on the approximate flux value ${S_{Rh}}$.

## 5 Results and Discussion

*et al.*(2017):

*P*molar flow rate at $x=0$ is shown (a scale factor ${10^{8}}$ is used) when the boundary condition for the substrate

*S*is a stepwise function (a scale factor ${10^{3}}$ is used):

##### (23)

\[ S(X,t)=\left\{\begin{array}{l@{\hskip4.0pt}l}0,\hspace{1em}& t<0.5;\\ {} 5,\hspace{1em}& t\geqslant 0.5.\end{array}\right.\]*S*and

*P*functions is used in presentation of results in all computational examples.

##### Fig. 1

*P*) molar flow rate at $x=0$, when for the boundary condition (23) is applied for the substrate (

*S*).

##### Remark 1.

### 5.1 Inverse Reconstruction of the Boundary Condition

*S*. The first one defines a piecewise linear function The second test boundary condition is defined as

##### Fig. 3

#### 5.1.1 Sensitivity of the Inverse Reconstruction Procedure to Perturbations of $Q(t)$

*et al.*, 2012). Yet, for the given mathematical model the obtained above theoretical stability estimates of the stationary solutions show that well-posedness of the proposed feedback control algorithm can be expected. We also note that the obtained system of two nonlinear PDEs with the prescribed four classical boundary conditions defines a stable mathematical model. Thus the discrete scheme is stable with respect to approximation errors.

##### Fig. 4

#### 5.1.2 Stability Analysis of the Control Scheme

### 5.2 Feedback Control of Different Treatment Protocols

##### Fig. 6

#### 5.2.1 Robustness of the Control Method

*μ*and different parameters of the model which are distributed within some compact set. Using the maximum principle which is valid for the given mathematical model we propose to use the maximum value of the control parameter

*μ*computed for the given set of parameters. In computational experiments we have fixed the the control parameter $1/\mu =181.88$ and used it for different values of $V\in [0.0005,0.0011]$ and ${K_{M}}\in [0.2,0.3]$. The obtained results have proved that in all cases the produced drug rate ${Q_{R}}(t)$ was close to the required treatment protocol $Q(t)$.