AnesthesiaGUIDE: a MATLAB tool to control the anesthesia

The goals of this paper are: (a) to investigate adaptive and fractional-order adaptive control algorithms for an automatic anesthesia process, using a closed-loop system, and (b) to develop an easy-to-use tool for MATLAB/Simulink to facilitate simulations for users with less knowledge about anesthesia and adaptive control. A model reference adaptive control structure was chosen for the entire system. First of all, to control the patient’s state during the surgery process, the patient mathematical model is useful, or even required for simulation studies. The pharmacokinetic/pharmacodynamics (PK/PD) model was determined using MATLAB’s SimBiology tool, starting from a previously available block diagram, and validated through simulation. Then, to achieve the desired control performances, two controllers are designed: a PI adaptive controller and a PIλ (PI-fractional) adaptive controller, using the MIT algorithm. The time response during anesthetic drug infusion for each patient can be plotted with the AnesthesiaGUIDE tool, which is also designed in MATLAB/Simulink. The tool was tested on data from 12 patients, subjected to general anesthesia, with successful results. Through this tool, the article provides a good opportunity for any user to experience with adaptive control for the anesthesia process.


Introduction
In medicine, general anesthesia plays a very important role during surgical interventions. General anesthesia is a drug-induced state and is achieved through the delivery of different substances which induce hypnosis, muscles relaxation and analgesia [1]. The process is considered a complex one and can be done manually or automatically. Even today, when the process is mostly done automatically, the anesthesiologist must continuously observe the patient's state in order to react to not desirable effects caused by over or under-dosing [2].
In the case of automatic anesthesia, a closed-loop structure is used to control the depth of hypnosis (DOH) in realtime. The drug effect is measured, and its value is compared with a DOH reference value. The difference between these two, namely the control error, is taken into account by the controller so that the drug infusion is adjusted in such a way that the DOH reference value is timely achieved and the patient is safe during the surgery [2][3][4][5].
In the literature, many closed-loop anesthesia systems are proposed. Some relevant examples may be mentioned here: Ilyas proposed a sliding-mode control algorithm [6] [7], Khaqan developed a nonlinear strategy [8], Navarro designed a fractional-order closed-loop model reference adaptive control based on a fractional-order model [9], Mendez used fuzzy predictive control techniques [10], Caiado investigated H-inf robust control [11], Reznanian developed a predictive controller [12], Usharani developed a PID controller strategy [16], Pasin performed a metaanalysis comparing closed-loop with manual drug delivery [17], Sopasakis built a fractional pharmacokinetics model [18], and Gonzalez-Cava proposed a study where an adaptive model deal simultaneously with the variabilities in the clinical response of the patients [19]. Actual trends and challenges are presented by Ghita in [20] and Zaouter in [21]. All the approaches guarantee a steady and safe hypnosis level.
In this paper, a model reference adaptive control (MRAC) system for anesthesia is developed based on two adaptive controllers: proportional-integrative (PI) and fractional proportional-integrative (PI λ ). The control performances achieved by the system are evaluated through simulations, and the efficiency of the proposed methods is proved. Also, a MATLAB/Simulink software tool named "AnesthesiaGU-IDE" was developed in order to support researchers and students with less knowledge in MATLAB or adaptive control theory to experience and verify these methods.
The current paper is structured in four sections: i) the first section is the current introduction and the appropriate context setup, ii) in section II a new pharmacokinetic/ pharmacodynamics (PK/PD) model is developed using MATLAB, iii) section III presents the control strategy with which the desired system's performances are achieved, and finally iv) the conclusions are drawn in section IV.

The patient model
To derive the control algorithm for the patient's state during a surgery process and to simulate the control system, a PK/ PD patient's mathematical model is necessary. In this work, a four-compartment model of the patient hypnotic state under propofol (hypnotic drug) infusion rate was developed using MATLAB's SimBiology toolbox [13]. It is based on a block diagram available in many related papers from the literature [9][10][11][12][13][14].
The obtained differential equations describing the PK/PD model are similar to those presented by Schneider in [9][10][11][12][13][14]: The first three equations represent the PK part of the model, while the last one is the PD part. The meaning of each parameter in (1) is detailed in Table 1.
Because of the negligible volume of the effect-site abstract compartment, the constants k 1e and k e0 are equal [9,14]. The PK parameters depend on the patient's characteristics only such as age, gender, height, and weight [9,10,14,15].
To measure the DOH, the bispectral index (BIS)-a signal derived from the electroencephalogram (EEG)-is used. The BIS value range is from 0 to 100. If the value of BIS is equal to 0, on the EEG a flat line can be seen (no cerebral activity), and if the BIS value is equal to 100, it means that the patient is fully conscious. The ranges between 60-70 and 40-60 represent high and moderate hypnotic conditions, respectively. During the surgery, the target value for the BIS is 50 and it is recommended to be kept between 40 and 60 [6][7][8][9].
The values of the parameters EC 50 and γ are patient specific.
Applying the Laplace transform on Eq. 1, the transfer function of the PK/PD model, which represents the complete patient model, is obtained: x e (t) + EC 50 Table 1 Nomenclature for the PK/PD model [9,14] The BIS value can be calculated from the concentration of the effect-site compartment by using the nonlinear Hill function: Symbol and unit Description Amount of drug in the peripheral compartments Propofol infusion rate k 12 , k 13 , k 21 , k 31 [min −1 ] Flow rates between compartments (positive constants) k 10 [min −1 ] Drug elimination rate k 1e [min −1 ] Rate constant of the effect-site compartment k e0 [min −1 ] Elimination rate constant from the effect-site compartment Drug concentration of the effect-site compartment E 0 Awake state; typically is 100 E max The maximum effect achieved after drug infusion EC 50 Patient's sensibility to the drug Sets the degree of function nonlinearity where the coefficients are: b 1 = k 21 + k 31 , b 2 = k 21 k 31 , a 1 = k 10 + k 12 + k 13 + k 1e + k 21 + k 31 , a 2 = k 10 k 21 + k 10 k 31 + k 12 k 31 + k 13 k 21 + k 1e k 31 + k 12 k 31 , a 2 = k 10 k 21 k 31 + k 21 k 31 k 1e . Based on the newly introduced parameter k 1e by the SimBiology toolbox in the PK model, a sensitivity analysis is performed to determine which parameter has a significant influence on the overall behavior of the model. This is examined by checking the magnitude of the computed sensitivities integrated over time (Fig. 1). It can be seen that the PK model is very sensitive to the k 1e parameter.
To validate the model presented in Eq. 3, the time responses for several patients to a constant input are plotted and analyzed (Fig. 2). The patients' clinical parameters (i.e., age, weight, height, and gender) are taken from [14], and the values for the pharmacokinetic parameters are calculated considering the equations proposed by Schneider and presented in many papers [8,9].
The patients' responses are different for the same propofol drug dose (50 mg/min). In the absence of a control loop [14], their time response is normal to tend to zero. The same behavior is obtained with the new model as well (Fig. 2), so the model can be validated.
Some patients have a higher sensitivity to the drug (the response decreases faster, e.g., patients 4 and 5), while others have less sensitivity (the response decreases slower, e.g., patients 1 and 2). Such behaviors are strongly influenced by the EC 50 EC 50 and parameters. The new parameter introduced in the PK model can adjust the patients' sensitivity to the drug administration, which might be considered an improvement (e.g., some responses do not decrease so faster).

MRAC-PI and MRAC-PI λ controllers design
To control the anesthesia, a model reference adaptive control structure is proposed (Fig. 3). This structure is preferred for two reasons. First, there are two loops that can assure the system's control performances imposed by the reference model (chosen by the user) and the stability of the closed-loop system. Second, the patient's behavior under drug infusion can be analyzed without recalculating the controller's parameters for each patient; the values will be adapted automatically for a new patient by the adjustment mechanism in the system.
To achieve this objective, i.e., to maintain the BIS level at 50, two adaptive control laws were proposed: an adaptive PI and an adaptive PIλ (PI fractional) control law, respectively. The two PI controllers are described by the following transfer functions: where k p and k i are the controller's parameters and α is the fractional order [22]. In both cases, the controller's parameters must be adapted by an adjustment mechanism so that the system's performances are fulfilled. For the adjustment mechanism, the MIT rule (gradient method) is going to be used [22]: where θ is a vector containing the controller's parameters, = [k p , k i ] ; γ is the adaptation gain; J( ) = (1∕2)e 2 is the minimization criterion, chosen so that the controller's parameters are adjusted in the direction of a negative gradient; e is the error between the system's output and reference model output (very small).
Further on, the adjustment mechanism proposed in (6) is modified for the PI λ controller's parameters (by introducing the fractional-order calculus [22]): Considering the patient's model, by its transfer function presented in (3), and the controller's transfer functions presented in (4)-(5), the closed-loop transfer functions are: Based on the Eqs. 8 -9 and considering the required system's performances (i.e., the overshot as small as possible and the settling time around 180 s), the transfer function of the reference model is chosen as follows: Basically, the performances are imposed by the term (60 s + 1) of the denominator. All the other terms are chosen to have little or no influence over the system behavior. (The poles and zeros are far away from the dominant pole 1/60.) Having Eq. 9 and considering the adjustment mechanism from Eq. 6, the variation of the controller's parameters k p and k i is [15]: If the PI λ controller is used, the adjustment mechanism for k p and k i parameters is similar to the one presented in Eq. 11. Hence, G 0 (s) = B m (s) A m (s) = b m0 s 3 + b m1 s 2 + b m2 s + b m3 s 5 + a m1 s 4 + a m2 s 3 + a m3 s 2 + a m4 s + a m5 = (s + 1)(2s + 1)(3s + 1) (60s + 1)(10s + 1)(11s + 1)(12s + 1)(13s + 1) Based on the reference model action and together with the proposed controllers, the system's performances and the stability in closed loop are guaranteed, i.e., the BIS value to be around 50 and the error to lead to zero.

AnesthesiaGUIDE -a MATLAB tool
To test the behavior of each patient under drug infusion, a MATLAB tool, named "AnesthesiaGUIDE", was developed. By using it, any user can experience the methods proposed in this paper and run simulations on the presented control scheme, less knowledge on MATLAB/Simulink and adaptive control theory being required.
The AnesthesiaGUIDE tool opens when the command "AnesthesiaGUIDE" is entered in the MATLAB's command window. The main window of the application (Fig. 4) is divided into five sections: -the first one is named "Patient parameters. " In this section, the user must fill the values of a few patient specific parameters (age, weight, height, and gender); after pressing the "Calculate" button, the values from section three, "Coefficients" are displayed. These values represent the PK/PD model, and they are read-only, as the section panel becomes inactive. During the simulation, they cannot be modified; -the second one, called "BIS, " is for the BIS coefficient; -section four, "Controller, " is used to control the system; the two types of controllers are presented; if the fractional PI controller is used, the value of must be set; -the system's output and the error can be seen for one patient after the button "START" is pressed, in section five, "Simulations," the button becomes active after the values from the previous panels are completed and the simulations will run in real time until the button "STOP" is pressed. Note: to not alter the values used, during the simulation the user interface is inac- tive, excepting the aforementioned buttons; at the end of the simulation, the patient state and the error can be seen in the graphs.
The application is available for download at https:// ati. unitbv. ro/ anest hesia guide/, where a short description and tutorial are also available.

Simulations and results
By using this tool, the state evolution for several patients under anesthesia is analyzed. The tool was tested on the clinical parameters of 12 patients from [14], but for better visibility, the results are presented only for six patients. Figure 5a shows the response of each patient under drug infusion ( p = 5 , i = 0.00001 ). For all patients, the induction phase achieves the value of 50 for the BIS index, in nearly 50 s. On some of the patients, the BIS time response has an oscillatory behavior (except patient 2). Then, the maintenance phase (steady-state regime) begins. The evolution of the error is depicted in Fig. 5b. The error leads to zero, so the stability in a closed-loop is guaranteed.
In these simulations, the PI adaptive controller was used. The results show the controller's capability to achieve and maintain the desired hypnosis level. The duration of the transient responses (induction phase), 50 s, is lesser than in other studies from the literature [7,9,12].
The results obtained when the PI λ adaptive controller is used are presented in Figs. 6 and 7, respectively. In Fig. 6a, all patients achieve the hypnosis level without any oscillations, which demonstrates the improved efficiency of the fractional-order controllers (with λ = 0.5 in these case). However, by increasing the value of λ (see Figs. 6b, c), the control performances are maintained while having small oscillations or overshoot.
By comparing the results obtained with the proposed controllers, we may conclude that better results are obtained when a PI λ controller is used. The IAE (Integral of Absolute Error) index is calculated to support this conclusion ( Table 1). As a remark, better performance (small overshoot), means a smaller value of IAE. Other basic criteria to evaluate the system's performances are the rise time and settling time. So, the ISE (Integral Square Error) and ITAE (Integral Time Absolute Error) indexes were also calculated ( Table 2). In some cases (e.g., patient 1), checking the ISE and ITAE values, a good rise time but poor settling time was obtained with the PI controller and vice versa with the PI λ controller (λ = 0.75).

Conclusions
The main contributions of this paper are to develop a new PK/PD model of the anesthesia process and to design a strategy to control it. Both are implemented in a tool called "AnesthesiaGUIDE", which is made for use with MATLAB/ Simulink. This tool offers the user an easy way to understand and simulate the control system, despite having less knowledge of control theory or anesthesia. Using these two methods, the hypnosis level is achieved and maintained. The tool can be improved by adding a performance analysis based on IAE, ISE, and ITAE criteria. So far, the tool cannot run simulations for several patients at once and some errors caused by setting wrong parameters are detectable while others not (e.g., if in the field 'Age' some letters are introduced, the user will be notified through an error message, but for the value 400 or -10 the simulation starts).
By deriving the PK/PD model from the block diagram in the literature using SimBiology, a new model was obtained that can adjust the patients' sensitivity to the drug administration.
Having the differential equations of the PK/PD model, two adaptive strategies to control the system are proposed. The results show that the approach presented to control the anesthesia offers good control performances. Town