Fractional-Order Closed-Loop Model Reference Adaptive Control for Anesthesia

Abstract

The design of a fractional-order closed-loop model reference adaptive control (FOCMRAC) for anesthesia based on a fractional-order model (FOM) is proposed in the paper. This proposed model gets around many difficulties, namely, unknown parameters, lack of state measurement, inter and intra-patient variability, and variable time-delay, encountered in controller designs based on the PK/PD model commonly used for control of anesthesia, and allows to design a simple adaptive controller based on the Lyapunov analysis. Simulations illustrate the effectiveness and robustness of the proposed control.

1. Introduction

In medical practice, the application of general anesthesia plays a significant role in a patient’s well-being, through the administration of a combination of drugs to provide adequate hypnosis (unconsciousness and amnesia to avoid traumatic recalls), paralysis, or muscle relaxation (to attain immobility, an absence of reflexes, and proper operating conditions), and analgesia (pain relief). This process is accomplished by an anesthesiologist who must continuously observe and adjust the rates and overall amounts of anesthetic agents delivered to the patient, preserving the stability of the autonomic, cardiovascular, respiratory, and thermoregulatory systems [1].

Even though been a subject of intense research in previous decades, the process of anesthesia is a complicated process and is still not well understood, resulting in a challenging control problem [2,3,4]. Moreover, in drug delivery systems, the controller has to tackle issues, such as inter- and intra-patient variability, multivariable characteristics, variable time delays, dynamics dependent on the hypnotic agent, model variability, and stability issues [5,6,7]. The current state of the art in the understanding of consciousness and the mechanisms of anesthetic-induced loss of consciousness is limited. Consciousness is very subjective and ethereal, so it is still very difficult to model. At present, the models available are the mean field models of drug action [5], which describe the phenomena presented in electroencephalograms (EEG) associated with different brain states [7].

There have been many attempts to automatize this process, aimed at the application of closed-loop control to drug delivery systems to assist anesthesiologists to improve the safety of the patient by avoiding excessive over-dosages and under-dosages in their patients, minimizing side effects and the risk of awareness during anesthesia [8]. Optimizing the delivery of anesthetics could lead to personalized healthcare, where the individual patient characteristics are taken into account for optimal and flexible drug administration. The first step towards an automated anesthesia process is to derive a mathematical model that adequately describes the system (or the experimental data). For control purposes, we need avoid very complex models with many parameters that cannot be determined or estimated independently mainly due to the lack of measurements and adequate sensors.

Commonly, the mathematical model used to study the depth of anesthesia is a PK/PD (pharmacokinetic/pharmacodynamic) model, with a third-order linear PK model and a PD model consisting of a first-order linear transfer function (which represents the time-lag between the drug infusion and the observed response) and a static nonlinearity, i.e., a Wiener structure [6]. Despite its plausibility being accepted by biomedical and control community, being a physiologically based empirical model, it presents many difficulties for controller designs, namely due to its large number of uncertain parameters resulting from significant variability among different individuals, time delay, lack of state measurements, and nonlinearity. The limited amount of real-time data and the poor excitation properties of the input signals constitute further challenges, making the identification of an online individualized model very difficult. Many attempts to simplify this model have been made, including PK/PD model structures with some fixed parameters [9,10], first-order plus time-delay models with an output nonlinearity [11,12], piecewise linear models [13], and low-complexity, control-oriented models [11,14].

One notable different approach was the recent introduction of fractional PK models [15,16]. In reference [17], the authors presented the relationship between the diffusion process and fractional-order models. Also, fractional-order pharmacokinetic models [15,16,18,19] that represent the experimental data in a more precise way thanks to the $t^{-\alpha }$ decay of the fractional operators, were introduced, and with this, a new line of investigation on the area of drug delivery systems has opened.

The problem of control of drug administration in anesthesia has been studied since the 1950s [20]. It is clear that the control of anesthesia has many challenges, such as multivariate features [21], different dynamics that are dependent on the drug and place of administration [22,23], stability problems [24], and performance of the control algorithm [25,26]. Given the complex nature and uncertainty of the process, it is not surprising that reliable models for the control of drug administration are not available.

The level of system uncertainty of fixed and robust gain controllers can unnecessarily sacrifice system performance, while adaptive controls can tolerate much higher levels of uncertainty and improve performance [27].

The control of anesthesia has been an active subject of research in past decades, and many control schemes have been developed, for example, fixed-gain controllers, mostly PID , were considered in references [25,28,29,30,31,32]. Adaptive controllers were developed in [33,34,35]. To deal with delays in the system, predictive controllers were used in reference [36,37,38]. Also, sliding mode control has been used in [39] and control based on neural networks in [40,41]. There have been a few published studies using a fractional control [2,42,43].

The implementation of such controllers is based on the assumption that the state measurements are available and that we have a clear and measurable output with little noise influence. However, in reality, the measured output may be noisy, and the system measurements do not produce this information directly. Instead, the state information needs to be inferred from the available output measurements. All of these challenges bring us the need to use estimation techniques that can estimate the state of each patient and adjust them based on the dynamics of each patient to deal with the system constraints [44].

This paper considers the problem of designing a fractional-order closed-loop model reference adaptive control (FOCMRAC) for anesthesia. To this end, three simple fractional-order model (FOMs) are proposed, which are shown to be able to follow the anesthesia dynamics of the PK/PD model. Then, based on these models, a simple FOCMRAC is designed. The stability of the closed-loop system is established using Lyapunov analysis. The main contributions of this paper are the proposal of the FOCMRAC scheme for fractional-order systems and its application to the control of anesthesia. Also, according to the best knowledge of the authors, the application of fractional-order controllers in control of anesthesia is very limited, so a novel and direct approach is proposed using fractional calculus tools.

Section 2 presents the mathematical tools used for the main results. Section 3 introduces the principles, the PK/PD model, and the challenges of the anesthesia process and presents the proposed FOMs. Section 4 gives the design of the FOCMRAC using the Lyapunov analysis. Section 5 gives the design of an identification scheme for fractional-order systems. Section 5 shows the simulation results using 30 virtual patients, and Section 6 presents the conclusions.

2. Preliminaries

In this section, the mathematical tools employed to develop the main results are briefly summarized.

The notion of fractional-order integral of order $\alpha >0\in \mathbb{R}$ is a natural consequence of Cauchy’s formula for repeated integrals, which reduces the computation of the primitive corresponding to the n-fold integral of a function $f\left(t\right)$ to a simple convolution [45]. The most common definition of a fractional-order integral is the Riemman–Liouville integral.

Definition 1.

The Riemann-Liouville integral is defined as

$${\displaystyle {}_{t_0}{D}_t^{-\alpha }f\left(t\right){=}_{t_0}{I}_t^{\alpha }f\left(t\right)=\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_{t_0}^t\frac{f\left(\tau \right)}{{(t-\tau )}^{1-\alpha }}d\tau ,}$$

(1)

where $\alpha \in {\mathbb{R}}_{+}$ and

$$\mathsf{\Gamma }\left(\alpha \right)={\int }_0^{\infty }{t}^{\alpha -1}{e}^{-t}dt.$$

(2)

is the Euler’s gamma function.

Among the existing definitions of the fractional derivative, there are three main definitions: the Grünwald–Letnikov, Riemman–Liouville and Caputo [46]. A review of the existing definitions of the fractional operator can be consulted in [47].

One of the advantages of Caputo derivative, which is a modification of the Riemann–Liouville definition, is that the initial conditions for fractional differential equations take the same form as those for integer-order differentiation which has a well-understood physical meaning. This definition is the one used in this paper.

Definition 2.

The Caputo derivative is defined as

$${\displaystyle {}_{t_0}{D}_t^{\alpha }f\left(t\right)=\frac{1}{\mathsf{\Gamma }(n-\alpha )}{\int }_{t_0}^t\frac{f^{\left(n\right)}\left(\tau \right)}{{(t-\tau )}^{\alpha -n+1}}d\tau ,}$$

(3)

where $n-1<\alpha <n\in {\mathbb{Z}}_{+}$.

The following property is established using the of Riemann–Liouville integral and Caputo derivative:

Property 1.

Newton–Leibniz formula generalization [48]

$${\displaystyle {}_{t_0}{I}_t^{\alpha }{}_{t_0}{D}_t^{\alpha }f\left(x\right)=f\left(t\right)-f\left(t_0\right).}$$

(4)

All modeling applications of the fractional calculus give rise to the notion of fractional-order models, namely, models described by an integrodifferential equation involving fractional order derivatives. Based on the Caputo derivative, a fractional-order system can be described as

$${}_{t_0}{D}_t^{\alpha }x\left(t\right)=f(x,t),$$

(5)

with initial conditions $x\left(t_0\right)$, where $\alpha \in (0,1)$, $f:\mathsf{\Omega }\times {\mathbb{R}}_{+}\to {\mathbb{R}}^n$ is piecewise continuous in t and locally Lipschitz in x on $\mathsf{\Omega }\times {\mathbb{R}}_{+}$, and $\mathsf{\Omega }\in {\mathbb{R}}^n$ is a domain that contains the origin, $x=0$.

Definition 3.

[49] The constant $x_0$ is an equilibrium point of the Caputo fractional dynamical system (5), if and only if $f(t,x_0)=0$.

When $\alpha \in (0,1)$, the fractional-order system (5) has the same equilibrium points as the integer-order system $\dot{x}=f(t,x)$ [50].

In nonlinear systems, the Lyapunov direct method provides an important tool to analyze the stability of a system without explicitly solving the differential equations.

The following method generalizes the idea for fractional-order systems which shows that a system is stable if there is a Lyapunov function for the system.

Theorem 1.

(Fractional-order extension of the Lyapunov direct method [50]) Let $x=0$ be an equilibrium point for the fractional-order system (5) and $\mathbb{D}\subset {\mathbb{R}}^n$ be a domain containing the origin. Let $V(t,x(t\left)\right):[0,\infty )\times \mathbb{D}\to \mathbb{R}$ be a continuously differentiable function and locally Lipschitz with respect to x, such that

$${\gamma }_1\parallel x\parallel \le V(t,x\left(t\right))\le {\gamma }_2{\parallel x\parallel }^{ab}$$

(6)

$${}_{t_0}{D}_t^{\alpha }V(t,x\left(t\right))\le -{\gamma }_3{\parallel x\parallel }^{ab}$$

(7)

where $t\ge 0$, $x\in \mathbb{D}$, $\alpha \in (0,1)$, ${\gamma }_1$, ${\gamma }_2$, ${\gamma }_3$, a and b are arbitrary positive constants. Then, $x=0$ is Mittag–Leffler stable. If the assumptions hold globally on ${\mathbb{R}}^n$, then $x=0$ is globally Mittag–Leffler stable.

Similar to integer-order systems, the idea of this fractional-order extension of the Lyapunov theorem is that the stability condition is derived by constructing a positive definite function, V, and calculating the fractional derivative of the V function.

Next is a version of Theorem 1 that is useful for the analysis of adaptive systems, especially when the ${}_{t_0}{D}_t^{\alpha }V(x,t)$ negative semidefinite is present.

Theorem 2.

[51] Let $x=0$ be an equilibrium point for the non-autonomous, fractional-order system (5). Let us assume that a continuous function exists, $V(x,t)$, such that

• $V(x,t)$ is positive definite.
• ${}_{t_0}{D}_t^{\alpha }(x,t)$, with $\alpha \in (0,1]$, is negative semidefinite.

Then, the origin of the system (5) is Lyapunov stable. Furthermore, if $V\left(x\right(t),t)$ is decrescent, then the origin of system (5) is Lyapunov uniformly stable.

One of the most used Lyapunov candidate functions to analyze the stability of integer-order system is the quadratic functions. However, in fractional cases, the use of these functions is not immediate, since evaluating the fractional derivative of the Lyapunov candidate function, in general, involves the evaluation of infinite sums, which include higher order integrals and derivatives of the fractional states of the system, which is not an easy task [52].

Next, some inequalities that facilitate the use of quadratic Lyapunov candidate functions in the analysis of stability of fractional-order system using Lyapunov’s direct method are presented.

Lemma 1.

[51] Let $x\left(t\right)\in {\mathbb{R}}^n$ be a vector of differentiable functions. Then, for any time instant $t\ge t_0$, the following relationship holds:

$$\frac{1}{2}{}_{t_0}{D}_t^{\alpha }\left(x^T\left(t\right)Px\left(t\right)\right)\le x^T\left(t\right)P_{t_0}{D}_t^{\alpha }x\left(t\right)$$

(8)

where $P\in {\mathbb{R}}^{nxn}$ is a constant, square, symmetric and positive definite matrix.

Lemma 2.

[51] Let $A\left(t\right)\in {\mathbb{R}}^n$ be a time-varying differentiable matrix. Then, for any time instant, $t\ge t_0$, the following relationship holds:

$${}_{t_0}{D}_t^{\alpha }\left[tr\left(A^T\left(t\right)A\left(t\right)\right)\right]⩽2tr\left(A^T{\left(t\right)}_{t_0}{D}_t^{\alpha }A\left(t\right)\right),\forall \alpha \in (0,1].$$

(9)

Barbalat’s lemma is a well-known and useful tool to deduce the asymptotic stability of integer-order nonlinear uncertain and time-varying systems (like adaptive systems). Due to some different properties between fractional-order integrals and integer-order integrals, it is not easy (but is possible) to use Barbalat’s lemma in fractional-order systems. In this paper, we will use the extension of this lemma for fractional-order systems that was presented in reference [43].

Lemma 3.

(Extension of Barbalat’s lemma [43]) Let $f:\mathbb{R}\to {\mathbb{R}}_{+}$ be a function that is uniformly continuous, and ${}_{t_0}{I}_t^{\alpha }$ is given by the Riemann–Liouville integral with $\alpha \in (0,1)$.

If ${\displaystyle \underset{t\to \infty }{lim}{}_{t_0}{I}_t^{\alpha }f\left(t\right)}$ exists and is finite, then $f\left(t\right)\to 0$, as $t\to \infty .$

This extension of the Barbalat lemma is very useful for concluding the convergence of the error in fractional-order adaptive systems.

Other important concepts are the positive realness (PR) and strictly positive realness (SPR). According to [53,54], these concepts for integer-order systems are valid for fractional-order commensurate systems with $\alpha \in (0,1)$.

Definition 4.

The $m\times m$ transfer function matrix $G\left(s\right)$ is called strictly positive real (SPR) if

1. 

All elements of $G\left(s\right)$ are analytic in $\mathcal{R}\ge 0$.

2. 

$G\left(s\right)$ is real for real s.

3. 

$G\left(s\right)+G^{T\ast }\left(s\right)$ for $\mathcal{R}\ge 0$ and a finite s.

A more extense review on fractional-order Lyapunov theory can be consulted in references [55,56].

3. Anesthesia

The effects of drugs on patients in the operating room vary with drug dosage, from patient to patient, and with time. Different doses of drugs result in different concentrations in various tissues, producing a range of therapeutic and sometimes undesirable responses. Responses depend on drug pharmacokinetics (time course of drug concentration in the body) and drug pharmacodynamics (the relationship between drug concentration and drug effect). These processes may be influenced by factors including pre-existing disease, age, and genetic variability. Patient responses to drugs may also dynamically altered by factors such as temperature, pH, circulating ion and protein concentration, levels of endogenous signaling molecules, and coadministration of the drug in the operating room environment [1].

General anesthesia consists of providing the patient with a reversible state of loss of consciousness (hypnosis), analgesia, and muscle relaxation to allow the patient to be operated on without pain by administering anesthetic drugs intravenously or through inhalation, providing maximum safety, comfort, and vigilance during the surgical act. The description of the general process with its variables is shown in Figure 1. As can be seen, the variables that can be manipulated are the anesthetics, relaxants, and serums. The disturbances in the system are signals that can occur at any time, such as surgical stimulation, blood loss, etc. The output variables are divided into measurable and non-measurable, and the main interest in the control of anesthesia is the non-measurable variables: hypnosis, analgesia, and muscle relaxation.

In practice, the anesthesiologist has to observe and control a large number of hemodynamic and respiratory variables as well as clinical signs that indicate the state of hypnosis and analgesia.

Most of the drugs used do not only operate on the desired effect but alter other aspects, for example, the effects of the anesthetic drug propofol not only affect the level of hypnosis but also increase the level of analgesia. The drug remifentanil, whose main objective is to increase the level of analgesia, has the same behavior, but as a side effect also increases the level of hypnosis. Due to this cross effect between these drugs, the anesthesiologist must adjust the desired levels of hypnosis and analgesia with different amounts of both drugs. From the point of view of control engineering, this problem is a problem of multiple inputs and multiple outputs.

In the daily work routine, the anesthesiologist calculates the amount of the necessary drug with the help of dose regimes given by the supplier of the drug, which, in most cases, are based on the patient’s body weight.

For a more detailed description of the process of general anesthesia, the reader is referred to references [1,6,57].

3.1. Pharmacokinetic/Pharmacodynamic Model

The PK/PD models most frequently used for propofol (hypnotic drug) are the fourth-order compartmental models (Figure 2) [58,59]. These models have been developed, tested, and validated on a wide range of real patient data and have often used in the literature for the control of anesthesia [60].

In this paper, we use the model presented in [58] given by

$${\displaystyle \dot{x}\left(t\right)=\left[\begin{array}{ccc}-(a_{11}+a_{21}+a_{31})& a_{12}& a_{13}\\ a_{21}& -a_{12}& 0\\ a_{31}& 0& -a_{13}\end{array}\right]x\left(t\right)+\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]u\left(t\right),}$$

(10)

where $x_1$ represents the concentration of the drug in the central compartment (intravascular blood); $x_2$ y $x_3$ represent the concentration of the drug in the peripheral compartments; $a_{12}$, $a_{13}$, $a_{21}$, $a_{31}$ are positive constants representing the flow between compartments; $a_{11}$ is the elimination rate of the drug through the metabolism; and $u\left(t\right)[mg/min]$ is the rate of infusion of anesthetic (propofol) in the central compartment. $a_{ij}$ are in $\left[mi{n}^{-1}\right]$ and $x_i$ are in $\left[mg\right]$.

An additional compartment, namely, the effect compartment, is introduced to represent the time delay between the observed effect and the plasma concentration. The effect compartment model links the plasma concentration (concentration in the central compartment) to the effect concentration with a first-order differential equation:

$$\dot{y}\left(t\right)=a_{eff}(x_1\left(t\right)-y\left(t\right)),y\left(0\right)=x_1\left(0\right),t\ge 0,$$

(11)

where $a_{eff}$ is the time constant; $x_1\left(t\right)$ is the concentration in the central compartment defined in (10); and $y\left(t\right)$ is the concentration of the effect compartment.

The pharmacokinetic parameters can be obtained through the following equations [58]:

$$\begin{array}{ccc}\hfill V_1& =& 4.27\left[l\right]\hfill \\ \hfill V_2& =& 18.9-0.391(age-53)\left[l\right]\hfill \\ \hfill V_3& =& 2.38\left[l\right]\hfill \\ \hfill C_{11}& =& 1.89+0-0456(weight-77)-0.06681(ibm-59)\hfill \\ & & +0.0264(height-177)[l/\min ]\hfill \\ \hfill C_{12}& =& 1.29-0.024(age-53)[l/\min ]\hfill \\ \hfill C_{13}& =& 0.836[l/\min ]\hfill \\ \hfill a_{11}& =& \frac{C_{11}}{V_1}\left[{\min }^{-1}\right];a_{12}=\frac{C_{12}}{V_1}\left[{\min }^{-1}\right];a_{13}=\frac{C_{13}}{V_1}\left[{\min }^{-1}\right]\hfill \\ \hfill a_{21}& =& \frac{C_{21}}{V_2}\left[{\min }^{-1}\right];a_{31}=\frac{C_{13}}{V_3}\left[{\min }^{-1}\right]\hfill \\ \hfill a_{eff}& =& 0.456\left[{\min }^{-1}\right]\hfill \\ \hfill lb{m}_m& =& 1.1·weight-128·\frac{{weight}^2}{{height}^2}\hfill \\ \hfill lb{m}_f& =& 1.07·weight-148·\frac{{weight}^2}{{height}^2}\hfill \end{array}$$

As we can observe in the previous equations, the PK parameters depend on the biometrical characteristics of the patient.

The bispectral index (BIS) is a signal derived from the electroencephalogram (EEG) used to assess the level of consciousness in anesthesia. A BIS value of 0 equals a flat line in the EEG, while a BIS value of 100 is the expected value of a fully conscious adult patient. The 60–70 and 40–60 ranges represent light and moderate hypnotic conditions, respectively. The target value during surgery is 50, giving us a gap between 40 and 60 to guarantee adequate sedation [6].

The BIS can be related to the concentration of the effect compartment by the nonlinear static function, termed the Hill equation [6]:

$$BIS\left(y\right)=E_0-E_{max}\frac{y^{\gamma }\left(t\right)}{y^{\gamma }\left(t\right)+E{C}_{50}^{\gamma }},$$

(12)

where $E_0$ denotes the base value (awaken state) and by convention, typically is given a value of 100; $E_{max}$ is the maximum effect achieved by drug infusion; $E{C}_{50}$ is the drug concentration to half maximal effect and represents the patient’s sensibility to the drug; and $\gamma$ determines the degree of nonlinearity of the function.

One of the significant challenges in the control of anesthesia is the variability among patients. This variability can occur as a result of patient physiology (age, gender, disease), variations in the PK processes (rates of absorption, distribution, metabolism, and elimination), and differences in PD (sensitivity of receptor, etc.) [61]. Also, in the medical practice, no state of this model is available for measurement, only the output (BIS) is measurable for feedback.

General anesthesia, as mentioned before, is comprised of three components: hypnosis, neuromuscular blockade, and analgesia. The three components can be modeled with a PK/PD model, but only hypnosis and neuromuscular blockade have clear physical measurements. In practice, the anesthesiologist measures the blood pressure, heart rate, sweating, and lacrimation to infer the level of analgesia in the patient. Therefore, there is not a clear index to measure the level of analgesia, but recently, some advances have been made [62]. In this paper we will concentrate on the hypnotic component of general anesthesia, noting that the control designed also can be applied in the control of neuromuscular blockade.

3.2. Challenges

Drug responses in humans are the result of integrated effects, including signal amplification and dampening mechanisms at cellular, tissue, and physiological system levels, and the pharmacodynamic responses can be therapeutic, toxic, or lethal.

In anesthesiology, it is easy to observe that the response to the same drug dose varies widely among patients. Part of the process of delivering anesthesia is titrating a drug dose to provide optimal therapy for a specific patient, particularly when the drug has a significant level of toxicity. At the same time, the anesthesiologist must know the dosing ranges that are appropriate for broad populations of patients to provide dosing guidelines.

The classical PK/PD models are coarse abstractions of a real distribution and elimination process. Still, these models describe the measured and observed concentrations. Another more complex type of model exists, called physiologically-based models [63]. With these models, it is possible to have more information available, for example, the influence of changes in cardiac output and organ blood flow on the time course of the concentration in the blood and various organs and tissues. Unfortunately, the development of such accurate models is expensive and can only be performed on animals. Furthermore, these accurate models are complex and have too many variables and parameters to be useful in the development of controllers.

The potential benefits of a closed-loop anesthesia delivery are more consistent drug administration, less inter- and intra-patient variability, less over- and under-dosages, faster control action to unexpected arousal (perturbation rejection), smaller quantities of drug usage, faster recovery of the patient, better hemodynamic control, and less hypotension during the induction of anesthesia. These factors are what keeps the interest in the control community to design reliable control schemes. Moreover, this problem is still considered an open problem.

The current state of the art understanding of the unconsciousness and the mechanisms of drug-induced unconsciousness is limited. Therefore, is very challenging to translate these little-known mechanisms into an accurate mathematical model.

One of the most critical challenges in the design of high-performance controllers for anesthetic drug delivery is the lack of accurate models.

Another challenge is the identification of a model from clinical data, as it has been shown [64] that the information available in the operating room (infusion rate of the drug and BIS index usually) is insufficient to identify a full-order PK/PD model. The excitation in the input signal is not sufficiently frequency rich and is not able to excite all modes in the model, because this input cannot be selected freely. Additionally, this model has a Wiener structure. Therefore, the task of identifying an individualized model online is very challenging.

Linear and nonlinear reduced-order model structures have been proposed to improve the identifiability and the control synthesis of the anesthesia process.

In [13], the authors presented a piecewise linear model, with which they used different LTI models to represent different phases of the process. For example, one model was used for the induction phase and another for the maintenance phase. Then, they synthesized a controller for each phase, similar to gain scheduling schemes.

In reference [65], a standard PK/PD model was proposed in which the authors assumed that only the PD parameters were responsible for the inter-patient variability, and an extended PK model was used and the model was linearized with the PD parameters identified via a Kalman filter.

A first-order plus time-delay was proposed in reference [10]. Here, the PK parameters are fixed, the Hill curve is linearized around an operating point, and the PD parameters are calculated using a standard least square algorithm.

In [66], a reduced-order model obtained using model reduction technics was obtained. This model is an affine model with only four parameters and takes advantage of the redundancy shown in the PK model, that is, the adjacent poles and zeros.

In [14] a MISO Wiener model was proposed for the pharmacokinetics and pharmacodynamics of propofol and remifentanil, this model uses a PK part with a reduced number of parameters (with a combination of three fixed parameters and one unknown), for the hypnotic drug and other for the analgesic. Also, a combine Hill equation (with a reduced number of parameters) that combine the effect of both drug in the level of unconsciousness.

Every proposed model has his particular advantages and disadvantage (which can be seen in detail in his respective references), but all these studies show us the necessity to develop a simple model for identification and control synthesis to circumvent the difficulties of using the PK/PD model with a Wiener structure.

Moreover, the majority of the control schemes based on the proposed models relied on the inversion of the non-linearity and supposed that the states are available for measurement, However, in practice the parameters of the non-linearity are unknown, and the states are not measurable online, thus adding more uncertainty in the control schemes presented.

3.3. New Modeling Paradigm: Fractional Calculus

Some researchers had proposed the necessity for a fractal view of physiology that explicitly takes into account the complexity of the living matter and its dynamics. Complexity in this context incorporates recent advances in physiology concerned with the application of concepts from fractal geometry, fractal statistics, and nonlinear dynamics, to the formation of a new kind of understanding within the life sciences.

The complexity of the human body and the characterization of that complexity through fractal measures and their dynamics involve the use of fractional calculus. Not only anatomical structures are fractal [67], such as the convoluted surface of the brain, the lining of the bowel, neural networks, and the placenta, but the output of dynamical physiologic networks are fractal as well [68].

The time series for the inter-beat intervals of the heart, inter-breath intervals, and inter-stride intervals have all been shown to be fractal or multifractal statistical phenomena. Consequently, the fractal dimension is a significantly better indicator of organismic functions in health and disease than the traditional average measures. The observation that human physiology is primarily fractal was first made in the 1980s [68]. Also, in reference [69], it was suggested that from the point of view of fractal physiology, blood flow and ventilation are delivered in a fractal manner in both space and time in a healthy body.

The control of physiologic variables is one of the goals of medicine, in particular, understanding and controlling physiological networks in order to ensure their proper operation.

Therefore, it seems reasonable that one novel strategy for modeling the dynamics and control of complex physiologic phenomena is through the application of fractional calculus [69].

A fundamental mechanism in the absorption of a drug in the human body is the process of diffusion. In reference [17], the authors present the relationship between the diffusion process and fractional-order models. Also, the introduction of fractional-order pharmacokinetic models are presented in references [15,16,18,19] which represent the experimental data in a more precise way, thanks to the $t^{-\alpha }$ decay of the fractional operators, and with this, a new line of investigation was opened ion the area of drug delivery systems.

In [15,70,71], it was suggested that biological systems (like in pharmacology and bioengineering) could be represented with a fractional-order model with a more simple structure compared with the integer-order counterpart, simplifying the control design by using a less complex model.

Other interesting applications of fractional calculus are shown in reference [72], where a fractional-order impedance model was used to model the respiratory mechanics. In reference [62] a fractional-order integro-derivative function to model the pain perception in analgesia is presented, and in reference [73] a review of the recent results reported on the model of biological phenomena using fractional calculus is presented.

Therefore, we can see that fractional calculus can offer us a new point of view to understand certain physical phenomena, especially those which, from the point of view of integer-order systems, seem too complex.

3.4. Proposed Model

Under a process where a considerable amount of uncertainty exists between individuals (inter-patient variability) and within the same individual (intra-patient variability), to have a valid deterministic model is challenging (in practice, it is almost impossible or too complex to be useful). So, if we could have a model for a single patient, this model is only valid for a period because of the changes in the physiological variables during surgery, which imply a change in the parameters of the model, hence the need to update the model online. Therefore, it would be convenient to have a generic model that has the ability of capture a wide range of dynamics (in this case, a range characterizing inter-patient variability) and adapted online to individualize the model for each patient.

Based on these facts, simple FOMs were considered. Three commensurate fractional-order models to represent the input–output behavior of the Wiener system (10–12) are proposed as follows:

$$G_1\left(\lambda \right)=\frac{b_0}{\lambda +a_0},$$

(13)

$$G_2\left(\lambda \right)=\frac{b_1\lambda +b_0}{{\lambda }^2+a_1\lambda +a_0},$$

(14)

$$G_3\left(\lambda \right)=\frac{b_2{\lambda }^2+b_1\lambda +b_0}{{\lambda }^3+a_2{\lambda }^2+a_1\lambda +a_0},$$

(15)

where $\lambda =s^{\alpha }$, s is the complex variable, and $\alpha$ is the commensurate order, with $0<\alpha <1$, $a_0$, $a_1$, $a_2$ and $b_0$, $b_1$, $b_2$ as the model parameters.

These control-oriented models should be evaluated based on the performance in the closed-loop, rather than from their prediction capabilities, because the latter is not the objective of the model, but ideally, a good prediction capability is also desired, as well as, not instead of, a closed-loop performance.

The proposed fractional models have a generic structure, and it can be shown that the input–output response of the patient models given by the PK/PD model can be captured by the fractional models proposed, given that the S-shape response of the patient model is part of all the possible responses of the proposed fractional models.

These models we can see them as phenomenological models (or some kind black-box model), namely, models that can capture the input–output dynamics of the patient model, but with the disadvantage of the loss of physical meaning of the model parameters.

It was shown in reference [74] that a set of parameters exists (depending on the respective structure) with which the fractional-order model can capture the response of a particular patient’s model. However, with the use of adaptive control, we aim to control a large set of patients with one controller, so that we do not need to identify a specific set of parameters for the FOM for a given patient—we only need a model structure capable of capture the overall response.

It is known that the particular response of a PK/PD model has an S-shape response, and reference [75,76] showed that simple structures like those proposed can capture this type of response.

Also, it is worth noting that one interesting application of the fractional operators and fractional differential equations is that of the study of the representation of an integer-order system by a fractional-order system. Fractional systems make it possible to carry out an efficient reduction of high-order integer-order systems, so we can represent those systems by a fractional system characterized by lower fractional order (compared with the high-order of the integer system) and a more simple structure [77,78].

4. Fractional-Order Closed-Loop Model Reference Adaptive Control

One recent approach of adaptive control that improves transient behavior by using a closed-loop architecture for reference models, named closed-loop model reference adaptive control (CMRAC), has been proposed [79]. In this approach, the focus is to have an adaptive system with output feedback. It has been shown that such closed-loop reference models can lead to a separation principle-based adaptive controller which is simpler to implement than classical ones based on observers or filtered signals. The simplification comes with the use of the reference model states in the construction of the regressor and not the classic approach where the regressor is constructed from filtered plant inputs and outputs.

Next, the generalizations of this scheme (Figure 3), denoted the fractional-order closed-loop model reference adaptive control (FOCMRAC), are presented. A preliminary version of this extension was presented in reference [80].

The commensurate fractional-order system, which is a state-space representation of the models (13–15) is given by

$$\begin{array}{cc}\hfill {}_0{D}_t^{\alpha }x\left(t\right)=& Ax\left(t\right)-B\Lambda u\left(t\right)\hfill \\ \hfill z\left(t\right)=& C^{T}x\left(t\right)\hfill \end{array}$$

(16)

where $x\in R^n$, $n\le 3$, $u\in R$, and $z\in R$. A and $\Lambda$ are unknown, B y C are known, and only z is available for measurement.

The control objective is to design a control signal, u, such that x follows the state, $x_m$, of the reference model given by

$$\begin{array}{cc}\hfill {}_0{D}_t^{\alpha }{x}_m\left(t\right)=& A_m{x}_m\left(t\right)-Br\left(t\right)-L(z-z_m)\hfill \\ \hfill z_m\left(t\right)=& C^T{x}_m\left(t\right)\hfill \end{array}$$

(17)

where $r\in R$ is the reference signal, and L is a feedback gain to be designed.

The following assumptions are made:

• The product, $C^{T}B$, is full rank.
• The pair, $A_m,C^T$, is observable.
• The system in (16) is at the minimum phase.
• ${\Theta }^{\ast }\in R^{nx1}$ exists such that $A+B\Lambda {\Theta }^{\ast T}=A_m$ and $K^{\ast }\in R$ such that $\Lambda K^{\ast T}=I$ (matching conditions).
• $\Lambda$ is diagonal with positive elements.
• The uncertain matching parameter, ${\Theta }^{\ast }$, and the input uncertainty matrix, $\Lambda$, have a priori upper bounds.

  $${\overline{\theta }}^{\ast }\triangleq sup\parallel {\Theta }^{\ast }\parallel and\overline{\lambda }\triangleq sup\parallel \Lambda \parallel .$$

  (18)

Next, some results needed for the design and analysis of the control scheme are presented. The proofs can be found in reference [79].

Lemma 4.

For the SISO case, in system (16), satisfying the suppositions 1–3, a $L_s$ exists, such that

$$C^T{(sI-A_m-L_s{C}^T)}^{-1}B=\frac{a}{s+\rho }$$

(19)

where $\rho >0$ is an arbitrary parameter and $a=C^{T}B$.

Lemma 5.

If $L_s$ is chosen as (19) and $M\triangleq C^{T}B$, the SISO transfer function , $M^T{C}^T{(sI-A_m-L_s{C}^T)}^{-1}B$, is SPR. Therefore, $P=P^T>0$ and $Q_s=Q_s^T>0$ exist, such that

$$\begin{array}{cc}\hfill {(A_m+L_s{C}^T)}^{T}P+P(A_m+L_s{C}^T)=& -Q_s\hfill \\ \hfill PB=& CM.\hfill \end{array}$$

(20)

Lemma 6.

Choosing $L=L_s-\rho B{M}^T$ where $L_s$ is defined by (19) andis arbitrary, $\rho >0$, the transfer function $M^T{C}^T{(sI-A_m-L{C}^T)}^{-1}B$ is SPR and satisfies

$$\begin{array}{cc}\hfill {(A_m+L{C}^T)}^{T}P+P(A_m+L{C}^T)=& -Q\hfill \\ \hfill Q\triangleq & Q_s+2\rho CM{M}^T{C}^T\hfill \end{array}$$

(21)

where P and $Q_s$ are defined in (20) and $M=C^{T}B$.

Assuming that L is chosen using the above lemmas, in the following theorem, the adaptive scheme is formulated.

Theorem 3.

Consider the fractional-order system given by (17) that satisfies assumptions 1–6 and the closed-loop model reference given by (17) and the control law

$$u\left(t\right)={\Theta }^T\left(t\right)x_m+K^T\left(t\right)r\left(t\right)$$

(22)

with adaptive laws

$$\begin{array}{cc}\hfill {}_0{D}_t^{\alpha }\Theta =& -\mathsf{\Gamma }{x}_m{e}_y^T\hfill \\ \hfill {}_0{D}_t^{\alpha }K=& -\mathsf{\Gamma }r{e}_y^T\hfill \end{array}$$

(23)

where $e_y=C^{T}e$, $\mathsf{\Gamma }>0$. Then, all the signals in the closed-loop system given by (16), (17), (22), and (23) are bound, $\forall t\ge 0$. Futhermore, the tracking error is $e\to 0$ when $t\to \infty$.

Proof. 

As shown in (16) , (17), and (22), the dynamic of the error, $e=x-x_m$, is given by

$$\begin{array}{cc}\hfill {}_0{D}_t^{\alpha }e=& (A_m+L{C}^T)e+B\Lambda ({\tilde{\Theta }}^T{x}_m+{\tilde{K}}^{T}r)\hfill \\ \hfill e_y=& C^{T}e.\hfill \end{array}$$

(24)

Consider the Lyapunov candidate function

$$V=\frac{1}{2}{e}^{T}Pe+\frac{1}{2}tr\left(\frac{{\tilde{\Theta }}^{T}P\tilde{\Theta }}{\mathsf{\Gamma }}\right)+\frac{1}{2}tr\left(\frac{{\tilde{K}}^{T}P\tilde{K}}{\mathsf{\Gamma }}\right),$$

(25)

taking the Caputo derivative and applying lemmas 1 and 2

$${}_0{D}_t^{\alpha }V\le e^{T}P{}_0{D}_t^{\alpha }e+tr\left(\frac{{\tilde{\Theta }}^{T}P{}_0{D}_t^{\alpha }\tilde{\Theta }}{\mathsf{\Gamma }}\right)+tr\left(\frac{{\tilde{K}}^{T}P{}_0{D}_t^{\alpha }\tilde{K}}{\mathsf{\Gamma }}\right).$$

(26)

Substitute (24), (20), and (21) to form

$$\begin{array}{cc}\hfill {}_0{D}_t^{\alpha }V\le & -e^{T}Qe+e^{T}PB{\tilde{\Theta }}^T{x}_m+e^{T}PB{\tilde{K}}^{T}r+\hfill \\ & +tr\left(\frac{{\tilde{\Theta }}^{T}P{}_0{D}_t^{\alpha }\tilde{\Theta }}{\mathsf{\Gamma }}\right)+tr\left(\frac{{\tilde{K}}^{T}P{}_0{D}_t^{\alpha }\tilde{K}}{\mathsf{\Gamma }}\right).\hfill \end{array}$$

(27)

Using the properties of the operator, $tr(\ast )$, and the fact that $e_y=C^{T}x$ and taking the adaptation laws as (23), we have

$${}_0{D}_t^{\alpha }V\le -e^{T}Qe.$$

(28)

Given that $\gamma >0$ and ${}_0{D}_t^{\alpha }V$ is negative semidefinite, from Theorem 2 the stability of the closed-loop system can be concluded. Applying the fractional integral and property 1 to (28), we have

$${}_0{I}_t^{\alpha }{e}^{T}Qe\le V\left(0\right).$$

(29)

Then, the fractional integral (29) exits, and in accordance with Lemma 3, it is concluded that the tracking error $e\to 0$ where $t\to \infty$. □

5. Fractional-Order Parameter Identifier

The next identification scheme is an extension for fractional-order systems of the scheme presented in [27]. This identification scheme does not need a state measurement and only uses information about the input and output to construct the identifier. This scheme will be used to show that the fractional model proposed in the next section captures the behavior of the PK/PD model.

Consider the commensurate fractional-order SISO system:

$$\begin{array}{cc}\hfill {}_0{D}_t^{\alpha }x=& Ax+Bu,\hfill \\ \hfill z=& C^{T}x,\hfill \end{array}$$

(30)

where $x\in R^n$, $n\le 3$, and only $y,u$ are available for measurement. System (30) can be written as

$$z=\frac{B\left(\lambda \right)}{A\left(\lambda \right)}u=C^T{(\lambda I-A)}^{-1}Bu,$$

(31)

where $\lambda =s^{\alpha }$ and $\alpha$ is the commensurate order with $0<\alpha <1$, and $A\left(\lambda \right),B\left(\lambda \right)$ are in the forms

$$\begin{array}{cc}\hfill A\left(\lambda \right)=& {\lambda }^n+a_{n-1}{\lambda }^{n-1}+\cdots +a_1\lambda +a_0,\hfill \\ \hfill B\left(\lambda \right)=& b_m{\lambda }^m+b_{m-1}{\lambda }^{m-1}+\dots b_1\lambda b_0,\hfill \end{array}$$

(32)

where the constants $a_i$ and $b_i$ for $i=0,1,2,\dots ,n-1$ are the system parameters, $n\le 2$ and $m\le 1$. Now, consider the linear parameterization of (31):

$$z=W\left(\lambda \right){\theta }^{\ast T}\psi ,$$

(33)

where ${\theta }^{\ast T}$ is the parameter vector, $\psi$ is the regressor that contains the filtered measurable signals, $u,y$, and $W\left(\lambda \right)$ is a proper stable transfer function.

$$\begin{array}{ccc}\hfill {\theta }^{\ast }& =& {[b_m,b_{m-1},\dots ,b_0,a_{n-1},a_{n-2},\dots ,a_0]}^T,\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill z& =& \frac{1}{\Lambda \left(\lambda \right)}y=\frac{{\lambda }^n}{\Lambda \left(\lambda \right)}y,\hfill \end{array}$$

(35)

$$\begin{array}{ccc}\hfill \varphi & =& \left[\frac{{\beta }_{n-1}^T\left(\lambda \right)}{\Lambda \left(\lambda \right)}u,-\frac{{\beta }_{n-1}^T\left(\lambda \right)}{\Lambda \left(\lambda \right)}y\right],\hfill \end{array}$$

(36)

$$\begin{array}{ccc}\hfill \Lambda \left(\lambda \right)& =& {\lambda }^n+{\eta }_{n-1}{\lambda }^{n-1}+{\eta }_{n-2}{\lambda }^{n-2}+\cdots +{\eta }_0\hfill \end{array}$$

(37)

$$\begin{array}{cc}=& {\lambda }^n+{\eta }^T{\beta }_{n-1}\left(\lambda \right).\hfill \end{array}$$

(38)

With the parameterization (33), the signals z and $\psi$ can be generated only with the information of u and y.

Because ${\theta }^{\ast }$ is a constant vector, we can write (33) in the form $z=W\left(\lambda \right)L\left(\lambda \right){\theta }^{\ast T}\varphi$, where $\varphi =L^{-1}\left(\lambda \right)\psi$.

$L\left(s\right)$ is chosen so that $L^{-1}\left(\lambda \right)$ is a proper stable transfer function, and $W\left(\lambda \right)L\left(\lambda \right)$ is a proper SPR transfer function.

A state-space representation of the parameterization (33) is given by

$$\begin{array}{c}{}_0{D}_t^{\alpha }{\varphi }_1={\Lambda }_c{\varphi }_1+lu,\hfill \\ {}_0{D}_t^{\alpha }{\varphi }_2={\Lambda }_c{\varphi }_2+ly,\hfill \\ z=y+{\beta }_{n-1}^T{\varphi }_2={\theta }^{\ast T}\varphi ,\hfill \end{array}$$

(39)

where

$${\Lambda }_c=\left[\begin{array}{cccc}-{\eta }_{n-1}& -{\eta }_{n-2}& \cdots & -{\eta }_0\\ 1& 0& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & 1& 0\end{array}\right],l=\left[\begin{array}{c}1\\ 0\\ ⋮\\ 0\end{array}\right].$$

(40)

Because $\Lambda \left(\lambda \right)=det(\lambda I-{\Lambda }_c)$ and $\Lambda \left(\lambda \right)$ is stable, it follows that ${\Lambda }_c$ is a stable matrix.

The state-space model (39) has the same input–output response as (30) and (31), provided that all state initial condition are $x_0=0$, ${\varphi }_1={\varphi }_2=0$.

Now the identification scheme is formulated in the following theorem.

Theorem 4.

Consider that the system (30), the linear parameterization (33) and the matrices $A_c$, $B_c$ and $C_c$ are associated with the state-space system , $W\left(\lambda \right)L\left(\lambda \right)=C_c^T{(\lambda I-A_c)}^{-1}{B}_c$, and the adaptive law

$${}_0{D}_t^{\alpha }\widehat{\theta }=\mathsf{\Gamma }ϵ\varphi ,$$

(41)

where $\mathsf{\Gamma }>0$ is the adaptation gain, and $\widehat{\theta }$ are the estimates ${\theta }^{\ast }$ in (33), and the identification error is given by $e=z-\widehat{z}$ and $ϵ=C_{c}e$. Then, the identification error is $ϵ\to 0$ as $t\to 0$.

Proof. 

The dynamics of the identification error are given by

$${}_0{D}_t^{\alpha }e=A_{c}e-B_c\tilde{\theta }\varphi ,$$

(42)

where the parameter error is defined as $\tilde{\theta }=\widehat{\theta }-{\theta }^{\ast }$.

Consider the following Lyapunov candidate function

$$V(\tilde{\theta },e)=\frac{1}{2}{e}^T{P}_{c}e+\frac{1}{2}{\tilde{\theta }}^T{\mathsf{\Gamma }}^{-1}\tilde{\theta },$$

(43)

where $\mathsf{\Gamma }={\mathsf{\Gamma }}^T>0$ is a constant matrix, and $P_c=P_c^T>0$ satisfies the algebraic equation

$$\begin{array}{c}\hfill P_c{A}_c+A_c^T{P}_c=-Q,\\ \hfill P_c{B}_c=C_c.\end{array}$$

(44)

By applying the Caputo derivative to (43), we get

$$\begin{array}{cc}\hfill {}_0{D}_t^{\alpha }V\le & e^T{P}_c{}_0{D}_t^{\alpha }e+{\tilde{\theta }}^T{\mathsf{\Gamma }}^{-1}{}_0{D}_t^{\alpha }\tilde{\theta }\hfill \\ \hfill \le & e^T[P_c{A}_c+A_c^T{P}_c]e-e^T{P}_c{B}_c\tilde{\theta }\varphi +{\tilde{\theta }}^T{\mathsf{\Gamma }}^{-1}{}_0{D}_t^{\alpha }\tilde{\theta }\hfill \\ \hfill \le & -e^{T}Qe-e^T{P}_c{B}_c\tilde{\theta }\varphi +{\tilde{\theta }}^T{\mathsf{\Gamma }}^{-1}{}_0{D}_t^{\alpha }\tilde{\theta }.\hfill \end{array}$$

(45)

From (45), we need

$$-e^T{P}_c{B}_c\tilde{\theta }\varphi +{\tilde{\theta }}^T{\mathsf{\Gamma }}^{-1}{}_0{D}_t^{\alpha }\tilde{\theta }=0.$$

(46)

From (44), we obtain $P_c{B}_c=C_c$, which implies that $e^T{P}_c{B}_c=ϵ$. Then, (46) can be written as

$$-{\tilde{\theta }}^T\varphi ϵ+{\tilde{\theta }}^T{\mathsf{\Gamma }}^{-1}{}_0{D}_t^{\alpha }\tilde{\theta }=0$$

(47)

which leads to (41).

Let the adaptation law be given in (41), and then

$${}_0{D}_t^{\alpha }V\le -e^{T}Qe.$$

(48)

Since (48) is negative semidefinite from Theorem 2, the stability of the closed-loop system can be concluded.

By applying the Riemann–Liouville integral in both sides of the inequality (48) and property 1, we have

$$I^{\alpha }{e}^{T}Qe\le V\left(0\right).$$

(49)

Then, the integral (49) exists and in accordance with Lemma 3, it is concluded that the identification error $e\to 0$ when $t\to \infty$. □

6. Simulations

In this section, numerical simulations are presented. The FOCMRAC scheme designed in the previous section was implemented in 30 virtual patients. Simulations tested the robustness of the control schemes against the intra-patient variability, inter-patient variability, disturbance, noise, and time delay. Also, this FOCMRAC design was compared with the scalar FOMRAC presented in [43].

To perform the simulations, we used the Matlab toolbox for fractional control NINTEGER v.2.3 developed in [81]. In this toolbox, various numerical approximations are implemented for the fractional-order operators that are explained in detail in reference [81]. The quality of the simulations is related to the approximation used and the associated parameters. In this work, we used the Crone approximation with the MacLaurin expansion, n = 10 and a bandwidth $[0.0011000]$. This approximation was chosen to make a trade-off between the quality of the simulation and the computing time.

6.1. Identification

To assess the variability among patients, 30 patient models (taken from different studies: patients 1–10 [82], patients 11–20 [36], patients 21–30 [83]) were used to emphasize the variability among the population. Figure 4 shows the responses of these models to a step input.

Table 1. Patients’ pharmacodynamic parameters and biometric features.

Patient	Age	Height	Weight	Gender	${\mathbf{EC}}_{50}$	$\mathit{\gamma }$
1	56	160	88	F	13.94	2.0321
2	48	158	52	F	13.88	1.0133
3	51	165	55	F	20	2.0196
4	56	160	65	F	20	1.8930
5	64	146	60	F	14.85	1.0702
6	59	159	110	F	20	2.6169
7	45	155	58	F	3.35	0.9172
8	51	163	55	F	12.17	1.8645
9	32	172	56	F	16.91	1.4517
10	68	160	64	F	15.52	0.9334
11	40	163	54	F	6.33	2.24
12	36	163	50	F	6.76	4.29
13	28	164	60	M	4.93	2.46
14	43	163	59	F	12.10	2.42
15	37	187	75	M	8.02	2.10
16	38	174	80	F	6.56	4.12
17	41	170	70	F	6.15	6.89
18	37	167	58	F	13.70	1.65
19	42	179	78	M	4.82	1.85
20	34	172	58	F	4.95	1.84
21	15	180.5	71	M	3.95	1.74
22	7	132	25.1	M	4.24	1.90
23	10	139	41.1	F	3.83	2.17
24	8	128	22	F	5.77	1.56
25	10	138	33.6	M	3.88	1.89
26	16	154.9	52.5	F	8.80	1.49
27	8	130	25.3	M	5.44	1.52
28	15	169	48	M	3.85	1.88
29	13	151	65	M	3.45	1.58
30	7	121	24	M	3.64	1.59

shows the pharmacodynamic and the biometric characteristics of the 30 patients.

As can be observed from Table 1, the parameters of the PK/PD models showed significant variation, depending on age, weight, height, and gender, causing substantial variation in the step response.

To verify the ability of the proposed FOMs to capture the dynamics of the PK/PD model, a simulation was carried out with the three structures proposed.

We used a nominal patient to carry out this simulation. Figure 5 shows the identification scheme, and Figure 6 shows the model’s outputs and the identification error. We can observe that the three fractional-order structures proposed can capture the step response of the PK/PD model. In Figure 7, the evolution of the model’s parameters is shown.

It can be seen that the input–output behavior of the PK/PD model was well captured by the proposed model structures.

It is worth noting that the objective of this study was the design of an adaptive control scheme, not to model a particular system, and, for adaptive control, in general, only a model structure is needed. Therefore, with the identification scheme and the simulation presented, it is shown that the proposed models were able to capture the input–output behavior of the presented example. Moreover, it was shown that the simple fractional-order models (13)–(15) are able to capture the responses of the patients, as illustrated in Figure 4.

6.2. Control

In this section, the effectiveness of the adaptive schemes designed is illustrated via simulations (Figure 8). It is worth noticing that, in the simulations, the plants (patients) are represented by the Wiener system (PK/PD model), and the proposed fractional-order models are only used to design and analyze the control schemes.

Table 2. Tuning parameters.

Control with 30 Patients
1st order FOMRAC	1st order FOCMRAC	2nd order FOCMRAC	3rd order FOCMRAC
$\alpha =0.005$	$\alpha =0.05$	$\alpha =0.005$	$\alpha =0.1$
$\gamma =0.0008$	$\gamma =0.005$	$\gamma =0.05$	$\gamma =0.05$
	$L=-30$	$L_1=-10$	$L_1=-10$
		$L_2=-2$	$L_2=-0.5$
			$L_3=-1$
Control with Perturbations
1st order FOMRAC	1st order FOCMRAC	2nd order FOCMRAC	3rd order FOCMRAC
$\alpha =0.005$	$\alpha =0.05$	$\alpha =0.005$	$\alpha =0.005$
$\gamma =0.0008$	$\gamma =0.0008$	$\gamma =0.001$	$\gamma =0.0008$
	$L=-0.1$	$L_1=-0.5$	$L_1=0.5$
		$L_2=-0.1$	$L_2=-0.5$
			$L_3=-0.05$
Time-Delay Robustness
1st order FOMRAC	1st order FOCMRAC	2nd order FOCMRAC	3rd order FOCMRAC
$\alpha =0.005$	$\alpha =0.05$	$\alpha =0.005$	$\alpha =0.1$
$\gamma =0.0008$	$\gamma =0.005$	$\gamma =0.01$	$\gamma =0.01$
	$L=-25$	$L_1=-10$	$L_1=-10$
		$L_2=-2$	$L_2=-0.5$
			$L_3=-1$

shows the values of the design parameters of the implemented control schemes. These values were chosen to make a trade-off between the speed of convergence and the transient performance.

The simulations are performed with the PK/PD model of anesthesia given by (10)–(12). The control objective was to take the patient to the level of $BIS=50$.

6.3. Inter-Patient Robustness

The inter-patient variability denotes the variation in the mathematical models among individuals. Every patient has a specific model.

The simulations were performed using 30 virtual patients applied to four different control schemes. The FOMRAC presented in reference [43] and the proposed FOCMRAC given by Theorem 3 were tested.

In the case of the FOCMRAC scheme, three controllers were implemented, each one based on the models (13)–(15) respectively, and applied to the 30 virtual patients.

Figure 9 shows the response of the 30 virtual patients with the four control schemes. We can observe that all schemes met the control objective–taking all of the patients to $BIS=50$. It can be seen that the controllers based on the second-order (14) and third-order models (15) had better performances than the controllers based on the first-order structure (13).

Figure 10 shows the control input of the four implemented controllers and it a similar behavior can be seen between the controllers.

In Figure 11, the tracking errors are shown. The errors of the FOMRAC scheme caused more oscillations in the induction phase and a higher convergence time than in the FOCMRAC schemes.

In Figure 12, Figure 13, Figure 14 and Figure 15 are shown the evolution of the controller parameters.

To evaluate controller performance, performance measures were calculated during the maintenance phase of anesthesia, in accordance with [31].

The performance error (PE) was calculated as the difference between the actual and the target values:

$$PE=100\left(\frac{BI{S}_{measured}-BI{S}_{target}}{BI{S}_{target}}\right).$$

(50)

The bias or median performance error (MDPE) described whether the measured values were either above or below the target values and thus represented the direction (under-shoot or over-shoot) of the PE:

$$MDPE=median\left(PE\right).$$

(51)

The inaccuracy or median absolute performance error (MDAPE) described the size of the errors:

$$MDAPE=median\left(\right|PE\left|\right).$$

(52)

The Wobble measured the intra-individual variability in the PE:

$$Wobble=median\left(\right|PE-MDPE\left|\right).$$

(53)

The global score (GS) was calculated according to the following equation:

$$GS=\frac{MDAPE+Wobble}{\%timeBISvaluebetween40and60}.$$

(54)

The overall control performances of the FOMRAC scheme [43] and the proposed FOCMRAC scheme are shown in

Table 3. Overall performance characteristics for the fractional-order adaptive controllers applied to 30 virtual patients. The values are reported as mean values and the minimum and maximum values within the 30 patients are presented.

Controller	1st Order FOMRAC	1st Order FOCMRAC	2nd Order FOCMRAC	3rd Order FOCMRAC
PE(%) [min, max]	2.0702 [−3.99, 6.38]	0.0803 [−1.74, 0.40]	0.1826 [0.09, 0.26]	0.1456 [0.06, 0.25]
MDPE(%) [min, max]	2.2505 [−1.94, 7.39]	0.2282 [−0.33, 0.12]	0.0986 [0.05, 0.21]	0.0759 [0.04, 0.15]
MDAPE(%) [min, max]	2.8269 [0.82, 7.34]	0.3478 [0.03, 3.24]	0.0986 [0.05, 0.21]	0.0759 [0.04, 0.15]
Wobble(%) [min, max]	0.5764 [−4$\times {10}^{-16}$, 4.92]	0.3196 [−6$\times {10}^{-16}$, 3.24]	0.0003 [−6$\times {10}^{-17}$, 0.005]	0.0028 [−6$\times {10}^{-18}$, 0.03]
Global score [min, max]	3.4033 [0.82, 12.26]	0.6674 [0.03, 6.48]	0.0986 [0.05, 0.22]	0.0787 [0.04, 0.16]

.

From Table 3, we can observe that the three FOCMRAC controllers led to an overall improved performance compared to the FOMRAC scheme. Moreover, the performance also improved depending on the model used, for example, the MDPE, MDAPE, Wobble, and GS indices decreased with an increase in the order of the controller, thus improving performance.

6.4. Robustness to Perturbations and Noise

During the maintenance phase, it is essential that the controller is capable of rejecting disturbances that occur during surgery.

The second simulation illustrated the robustness of the control scheme to perturbations and noisy measurements, specifically, those perturbations that affect the value of the BIS index in a patient. These perturbations occur because of, for example, intubation of the patient, painful stimuli, or blood loss. Figure 16 shows the artificial disturbance signal. Figure 17 shows the BIS response of patient 1 with the four control schemes. Note that the controllers are capable of compensating for these perturbations and noise, although the under-shoots in the responses are accentuated for the value of the adaptive gain and the lack of negative control. In Figure 18, the tracking error is shown, and Figure 19 shows the evolution of the controller’s parameters.

6.5. Robustness to Time Delay

The simulation illustrated the robustness of the control scheme to the time delay, represented by the parameter $a_{eff}$ in the effect compartment (11) of the PD part of the PK/PD model. In Figure 20, the step response of patient 1 is shown for different values of $a_{eff}$.

Figure 21, Figure 22, Figure 23 and Figure 24 shows the BIS output, the control input and the tracking error of patient 1 with different time delays, from approximately 0–8 min. It can be observed that despite the change in time delay, the four adaptive controllers were capable of compensating for the delay variation, thanks to the memory effect of the fractional operators.

6.6. Comparison between Fractional-Order and Integer-Order MRAC Schemes

In the last simulation, we made a comparison between the fractional-order MRAC schemes and integer-order MRAC schemes applied to the nominal patient (patient 1).

Figure 25 and Figure 26 shown the adaptive schemes based on the first-order model (13). We can see that the integer-order controller had an oscillatory response, and in particular, the MRAC scheme became unstable.

Figure 27 shows the CMRAC schemes based on the second-order model (14). We can observe that the integer-order controller had a constant oscillatory response around the reference level, and the fractional-order controller met the control objective without oscillations.

Figure 28 illustrates the response of the CMRAC schemes based on the third-order model (15). We can observe that the integer-order controller had a damped oscillatory response around the reference level but with a much larger control input, and the fractional-order controller met the control objective without oscillations.

These simulations show that a complex process like the control of anesthesia can be controlled and meet the control objective using simple fractional-order models, which is not possible with the same simple integer-order models and controllers.

7. Conclusions

In the case of control of anesthesia, we are dealing with a process that is not fully understood, and with the current integer-order model paradigm, dominated by the PK/PD approach, which, despite its plausibility and acceptance by the biomedical and control community, presents a significant challenge, not only due to the nature of the processes (unknown parameters, unknown time delay, states not available for measurement, positivity and poor excitation in the control input) but also due to the model structure (Wiener structure). Recent developments in biology and physiology with a focus on fractal dynamics, pharmacology, and anesthesia, using fractional calculus, have led to a new paradigm in the understanding and modeling of the physical phenomena in these areas. This new paradigm suggests that complex phenomena (including anesthesia) can be modeled more precisely and accurately and with a simple structure with fractional-order tools.

The main contributions of this paper are the extension of the CMRAC and the identification schemes for fractional-order systems and their application to the control of anesthesia.

It was proposed that simple fractional-order models could be used to represent the input–output behavior of the PK/PD model of anesthesia. It was shown that a simple, fractional-order structure can capture the response of the patient, which is represented by a nonlinear model (Wiener model). Based on the fractional models proposed, a FOCMRAC was designed to control the PK/PD model of anesthesia, showing, through simulations, that these controllers can meet the control objectives. Moreover, this scheme was shown to be robust against inter and intra-patient variability, time delay, parameter uncertainty, perturbation, and noise. These results represent a different and novel approach for attacking the problem of control of anesthesia, which still is an open problem and an active topic of research.