Fractional-Order Modeling of the Depth of Analgesia as Reference Model for Control Purposes

Abstract

Little research has been carried out in terms of modeling and control of analgesia. However, emerging new technology and recent prototypes paved the way for several ideas on pain modeling for control. Recently, such an idea has been proposed for measuring the Depth of Analgesia (DoA). In this paper, that solution is further exploited towards obtaining a novel fractional-order model and dedicated controller for DoA. First, clinical data from patients undergoing general anesthesia are used to determine a commensurate fractional-order model of the skin impedance at each sampling period. Second, we provide a proof of concept indicating that fractional order changes due to variations in the infused opioid drug (Remifentanil). Third, a fractional-order model for DoA is developed correlating the changes in the pain index (as the output signal) and the Remifentanil infusion rate (as the input signal). Standard optimization routines are used to estimate the parameters. A database of 19 real patients is used. Lastly, a preliminary fractional-order controller is designed and tested in simulation for the 19 patients. The closed-loop simulation results correspond to the expected clinical outcomes.

1. Introduction

General anesthesia is a core topic in current surgical procedures, with several research papers directed towards finding a simple and efficient solution for automatic control of drug dosing [1,2,3]. This is envisaged as the most suitable solution to avoid under- or over-dosage and to offer a dedicated support for decision making for the anesthesiologist [4,5,6]. General anesthesia consists of three main components: hypnosis, neuromuscular blockade and analgesia. Computer-based optimal drug dosage in general anesthesia is currently undergoing extensive research, mainly focused on the regulation of hypnosis and neuromuscular blockade. For a holistic approach, the regulation of analgesia is equally important.

Analgesia refers to the relief of nociception or pain (if the patient is awake) and is achieved by an interruption in the nervous system pathway between the sense organ and brain. The interruption can be performed using several types of analgesics, such as opioids, i.e., Remifentanil. The opioid drug infusion is currently decided by the anesthesiologist by monitoring several patient vital signs. Optimal drug dosing would be possible using computer-controlled analgesia. To achieve this, adequate nociception monitoring is crucial to give feedback information on the dose–effect response in patients. To date, this has received limited attention due to the lack of suitable devices for monitoring. Apart from the lack of feedback information, a model for analgesia is also required to properly design a computer algorithm for optimal drug dosing.

There is little research in terms of modeling and control of analgesia. This is mainly due to the limited choice of monitoring devices and signals to measure the Depth of Analgesia (DoA) as a relationship between the infusion rate of the opioid drug, Remifentanil, and a signal that evaluates the nociception level in patients. Although many studies on analgesia exist, as reported from clinical trials [7,8,9,10], there is no established mathematical model that relates Remifentanil to its effect upon a surgical or nociceptive stimulus. The purpose of this manuscript is to develop such a mathematical model, which, in turn, will be used to design a dedicated controller.

The emergence of devices that measure DoA or nociception levels [11,12] enables the development of preliminary mathematical models and of dedicated controllers. Two monitoring devices have emerged to be more suitable for estimating DoA based on indices such as Nociception Level (NOL) [13] or Skin Conductance (SC). These are the commercial devices Medasense (CE0344, Manufacturer: Medasense Biometrics Ltd., Ramat Gan, Israel)) and Medstorm (CE0413, Manufacturer: Med-Storm Innovation AS, Oslo, Norway) [13,14]. Both devices are built using standard electric circuitry and follow the classical laws of physics while delivering a signal (either NOL or SC) that is used to quantify a patient’s physiological pain response (nociception). Recent papers have compared these two commercial monitoring/measuring devices with a prototype AnspecPro monitor and have validated its utility in estimating nociception levels expressed in terms of various new indices [11,15]. Much like Medasense and Medtsorm, the AnpsecPro device is designed and follows the classical laws of physics, with the output being an electrical signal. For the AnpsecPro device, our previous manuscript has proposed several indexes, similar to the NOL and SC indexes, to quantify the nociception or DoA. At the same time, an analysis based on clinical data has been performed to establish a causality between the opioid (Remifentanil) infusion rate and the proposed nociception indexes. Simple mathematical models have been derived to link the input (Remifentanil) signal to the output (nociception index). These models are second-order integer transfer functions [15].

Nonetheless, fractional calculus, as a generalization of integer-order calculus, has been recently used more frequently to describe physiological phenomena in the human body, such as fractional-order models for tissues [16], cancer-spreading models [17,18,19,20,21], HIV dynamics [22,23], virus spreading [24,25], and drug diffusion [26], to mention just a few.

In this paper, we attempt to use fractional calculus to model DoA. Clinical data available and measured using the AnspecPro device are used as described in [12,15]. The output signal of the AnspecPro device is the measured skin impedance in a frequency range of [628.32, 9424.78] rad/s [27]. Skin impedance is measured over a certain period of time according to the clinical protocol [12,15] with a constant sampling period. First, various commensurate fractional-order (FO) models are determined for the skin impedance using frequency domain identification techniques and based on the mean frequency response. These are compared in terms of accuracy and complexity to establish the best structure of the impedance model. Next, at each sample, a commensurate FO model for the impedance is obtained. It is shown that the commensurate fractional order changes due to variations in the infused Remifentanil drug. The commensurate fractional order, denoted as q, is further considered a new index for estimating the DoA (output signal) and a correlation with the Remifentanil infusion rate (input signal) is established as a fractional-order transfer function. Validation with clinical available data is performed.

Once such a model for the DoA is developed, it can be used for the tuning of an optimal controller. Nowadays, most of the research for closed-loop control is centered around one of the key components in anesthesia, namely, depth of hypnosis. This is reflected in the Bispectral Index (BIS) signal, which is available via electroencephalography (EEG). Several closed-loop strategies have been developed and proposed for monitoring and maintaining the BIS signal within acceptable ranges. These range from simple proportional–integral–derivative (PID) control algorithms to more advanced solutions such as fractional-order controllers, event-based algorithms, model predictive control, to name just a few [1,3,4,5,6,28]. A second component of anesthesia that is matured in terms of closed-loop testing and analysis is the neuromuscular blockade [29]. In this case, the available literature has also introduced and validated several control possibilities, including, for example, PIDs and model predictive control [30,31]. Some high-quality review papers dealing with closed-loop control of anesthesia are also available [2,4,5,28]. For DoA, no control strategies have been proposed so far. In this manuscript, we propose a simple yet efficient and robust fractional-order controller for regulating the DoA. These fractional-order controllers are generalizations of the integer-order controllers [32,33,34] and they pose extra degrees of freedom, which increase the flexibility of the design [35,36,37] and allow for better closed-loop performance, especially in terms of robustness [32,33,34,38,39,40,41].

A feasibility analysis regarding closed-loop control is finally performed by designing this Proportional Integral–Fractional-Order Proportional Derivative (PI-FOPD) controller. The fractional-order transfer function with the commensurate fractional order q considered as the output signal and Remifentanil rate as the input signal is used as a model for analgesia in the tuning of the dedicated controller. In a practical setting, the AspecPro device measures the patient skin impedance in a certain frequency range, an electrical signal, which is used to estimate fractional-order models with commensurate order q(t), as indicated in Section 2.2 and Section 3.1. This occurs at each sampling interval. The signal q(t) is then considered by the controller as the output signal. Based on a reference value q*, considered to be equivalent to no pain, the fractional-order control algorithm generates a control signal that corresponds to an optimal infusion rate for Remifentanil. This control signal computed by the controller is an electrical signal fed to dedicated micropumps, which close the loop by supplying the required drug infusion rate.

Establishing a fractional-order model for the DoA is the main target of the current manuscript. The motivation behind this research consists of the growing trend towards a complete anesthesia control, which requires the control of not just the depth of hypnosis and neuromuscular blockade but that of analgesia as well.

The main innovation of the article consists of the following:

•

Development and validation of a novel commensurate fractional-order skin impedance model;

•

Development and validation of a novel ‘pain index’ and its correlation to the opioid (Remifentanil) infusion rate in the absence of surgical stimuli;

•

Development and validation of a fractional-order transfer function for analgesia suitable for closed-loop control;

•

Design of a simple preliminary fractional-order controller for analgesia to test its feasibility.

This manuscript is structured as follows. In the next section, a short overview of the patient database and the clinical trial is introduced, followed by the theoretical background for fractional-order modeling using frequency and time domain data. A glance at the theoretical concepts of the optimization algorithms used is also included in Section 2. The rationale behind using the proposed fractional-order models presented in this section is further explained in the next section, which also presents the results. Section 3 is organized in four main subsections. First, the results of the fractional-order impedance model are presented, followed by the development and validation of the fractional-order model for the depth of analgesia. An integer-order model is used for comparative purposes in the third subsection. Finally, the last subsection presents the necessary theoretical background regarding fractional-order control and tuning, along with the closed-loop simulation results for 19 patients in the database.

2. Materials and Methods

2.1. The Patient Database

Data from several patients under Total Intravenous Anesthesia were used in this research. The clinical protocol is described in detail in [12,15], where three non-invasive pain monitors (Medasense, Medstorm and AnspecPro) were evaluated in terms of their effectiveness to represent the nociception–anti-nociception sensitivity of the patients. Data referring to the DoA were collected with a Remifentanil target-controlled infusion (TCI model of Minto) to reach an adequate level of analgesia at the moment of intubation. A total of 70 patients were included in this study and, among these, the Medasense device was used to collect data from 26 patients, the Medstorm was used for 20 patients and the AnspecPro was used for 24 patients [12]. All patients in this study were included after their informed consent was received. Details regarding the inclusion/exclusion requirements are given in [12].

This trial was carried out in accordance with the defined protocol, which includes ISO 14155 [42] and is centered on abiding by applicable laws and good clinical practices. ISO 14155 is a worldwide standard for scientific and ethical excellence that regulates the design, conduct, record-keeping, and dissemination of research projects involving human subjects. Furthermore, this trial complied with the legal guidelines specified in European Regulation (EU) 2017/745 [43]. The Federal Agency for Medicines and Health Products of Belgium (FAGG) and the Ethics Committee of Ghent University Hospital approved this trial (EC/BC-08020, FAGG/80M0840, EudraCT: CIV-BE-20-07-0342442020, clinicaltrials.gov: NCT04986163, principal investigator: Martine Neckebroek). The reader is referred to [12,15] for more details on the clinical trial and how the data used in this research were collected.

A graphical representation of the protocol is given in Figure 1.

Useful data for developing a model for the DoA begin when the Remifentanil infusion (blue line) is administered to the patient and end with the start of the intubation. This corresponds to the analgesia region, depicted in Figure 1, which lasts on average for 145 s for the considered patients, until intubation. This is a major disturbance, which corrupts the data for identification and these cannot be used further. Data are collected with a sampling period of 5 s. As such, the quantity of data for each patient is limited and modeling becomes a difficult task. Methods to tackle this were proposed in [44].

Table 1. Biometric data for patients used in this work.

Patient No.	Sex	Age	Height (cm)	Weight (kg)	Data Until Intubation (Seconds)
3	F	42	168	70	86
10	M	73	181	83	157
14	X	19	168.7	59	250
21	F	62	168	88	221
24	F	36	168	63	95
25	M	58	179	94	240
33	F	46	170	82	115
35	F	41	167	65	100
38	F	18	163	64	70
39	F	64	158	58	134
41	F	30	156	44	90
48	M	31	170	85	142
54	F	21	166	63	140
56	M	63	168	82	135
59	F	71	167	83	86
63	X	26	170	50	170
64	F	33	172	94	165
66	F	70	165	62	155
70	X	20	155	49	206

shows the intubation start time relative to the beginning of the analgesic region considered as time zero.

In this research, only 19 of the 24 patients that were monitored using the AnspecPro device (set as MON2) were considered. Their biometric data are indicated in Table 1. Among these, 4 patients were male (M), 12 awee female (F) and 3 were identified as transgender (X). Patients ranged from 18 to 73 years old (with an average of 43.36), their height varied from 155 cm to 181 cm (with an average of height of 167.3 cm) and their weight varied from 44 to 94 kg (with the average of 70.4 kg).

2.2. Fractional-Order Model Identification Using Frequency Response Data

A generalized continuous time commensurate fractional-order model is represented by:

$$G\left(s\right)={\displaystyle \frac{b_n{s}^{nq}+b_{n-1}{s}^{\left(n-1\right)q}+\cdots {+b}_2{s}^{2q}+b_1{s}^q+b_0}{a_n{s}^{nq}+a_{n-1}{s}^{(n-1)q}+\cdots +a_2{s}^{2q}+a_1{s}^q+a_0}}$$

(1)

where n is the number of poles/zeros and q is maximum common multiple of powers of the Laplace operator s. Transfer functions of the form in (1) are defined as the commensurate fractional-order models. The transfer function in (1) will be used to model the skin impedance, as will be indicated in the next section. The choice of using (1) consists of its generalized form, which can be further altered to fit experimental data. Additionally, similar skin impedance models have been used, as will be further discussed in Section 3.

In this study, system identification based on the frequency response is performed. Given the known frequency behavior of (1) defined as G(j $\omega$), the fractional order q and the coefficients b_i and a_i, with i = 1,…, n are estimated by minimizing the square norm:

$$E=wG\left(j\omega \right)\left[a_n{\left(j\omega \right)}^{nq}+a_{n-1}{\left(j\omega \right)}^{(n-1)q}+\cdots +a_2{\left(j\omega \right)}^{2q}+a_1{\left(j\omega \right)}^q+a_0\right]-\left[b_n{\left(j\omega \right)}^{nq}+ b_{n-1}{\left(j\omega \right)}^{\left(n-1\right)q}+\cdots {+b}_2{\left(j\omega \right)}^{2q}+b_1{\left(j\omega \right)}^q+b_0\right]$$

(2)

where the weights are frequency dependent and computed according to:

$$w=\left\{\begin{array}{c}{\displaystyle \frac{{\omega }_2-{\omega }_1}{2{\omega }_1^2}}, \mathrm{i}\mathrm{f} i=1 \\ {\displaystyle \frac{{\omega }_{i+1}-{\omega }_{i-1}}{2{\omega }_i^2}}, \mathrm{i}\mathrm{f} 1<i<f\\ {\displaystyle \frac{{\omega }_f-{\omega }_{f-1}}{2{\omega }_f^2}}, \mathrm{i}\mathrm{f} i=f\end{array}\right.$$

(3)

assuming a set of available frequencies ${\omega }_i$, i = 1,…, f. To improve the quality of the approximation at low frequencies, the norm in (2) is a modification of the method of Levy with weights (3) [45]. To estimate the coefficients of (1), q and n need to be supplied.

2.3. Fractional-Order Model Identification Using Time Response Data

Input–output time domain data are also used in this research to determine fractional-order transfer functions in the following form:

$$G_{FO}\left(s\right)={\displaystyle \frac{k}{a_1{s}^{{\alpha }_1}+a_2{s}^{{\alpha }_2}+1}}$$

(4)

where k is the static input–output gain and ${\alpha }_1$ and ${\alpha }_2$ are fractional orders. The choice of model in (4) is explained in Section 3.2. The parameters in (4) are estimated by minimizing the error ${‖e\left(t\right)‖}_2^2$ given by $e\left(t\right)=y_{exp}\left(t\right)-y_{sim}\left(t\right)$, where $y_{exp}$ and $y_{sim}$ are the experimental output signal and the simulated one using (4). To simulate the fractional-order model in (4), the Grunwald–Letnikov method [46] is used. In order to obtain the best fit, the Trust Region Reflective algorithm is used. The algorithm requires an initial model guess. An integer-order model is supplied in this case, mathematically described as follows:

$$G_{IO}\left(s\right)={\displaystyle \frac{k}{c_1{s}^2+c_{2}s+1}}$$

(5)

where k, c₁ and c₂ are estimated using input–output data and the Matlab (R2019a) tfest function. Minimum and maximum values for the exponents ${\alpha }_1$ and ${\alpha }_2$ and the coefficients a₁ and a₂ are set as constraints for the optimization algorithm as follows: 1 < ${\alpha }_1<2$, 0 < ${\alpha }_2<1.5$, c₁/3 < a₁ < 3c₁ and c₂/3 < a₂ < 3c₂. In the Trust Region Reflective algorithm, which is used to estimate (4), the gain k is fixed and assumed to be that of (5) and only the remaining parameters ${\alpha }_1$, ${\alpha }_2$, a₁ and a₂ are estimated.

2.4. Optimization Algorithms

Two optimization routines were used to estimate the fractional-order transfer functions. The Nelder–Mead (NM) algorithm [47] was used to determine a best-fit fractional model satisfying the available frequency domain data, whereas the Trust Region Reflective (TRR) algorithm was employed to determine the best-fit fractional-order model satisfying the available time domain data.

The Nelder–Mead method is a Black-Box Optimization approach, where the optimization of an objective function is performed without any information regarding the analytic form and gradient of the objective [47,48]. The advantage of the NM algorithm is the quick convergence combined with a relatively small number of function evaluations. It achieves decent results compared to computationally expensive problems. With the Nelder–Mead simplex algorithm, an objective function of n design parameters can be minimized. The algorithm iteratively updates the worst design in order to converge on the smallest value of the objective function:

$$J={\sum }_{k=1}^n{\left(x_{\mathrm{e}\mathrm{x}\mathrm{p}}\right(k\left){-x}_{sim}\right(k\left)\right)}^2$$

(6)

where $x_{sim}\left(k\right)$ represents the frequency response of the fractional-order model for each frequency ${\omega }_k$ and $x_{exp}\left(k\right)$ is the measured frequency response. The algorithm sorts the solutions, reflects the worst one through the centroid of the remaining ones and checks convergence. Extension, contractions and shrinkage are also part of the iterative algorithm. The algorithm stops whenever convergence is achieved either in the parameters or in the objective. Convergence in the parameters is achieved whenever the largest differences between two adjacent solutions is less than a user-specified tolerance value ${ϵ}_x$ times the average of adjacent solutions. Convergence in the objective function is achieved whenever the difference between the worst and best objective function is less than a user specified tolerance ${ϵ}_f$ times the best objective function value [47]. In this study, the tolerance values are set to ${ϵ}_x={ϵ}_f={10}^{-4}$.

The Trust Region Reflective method [49] is one of the most popular and effective approaches for solving bounded constrained nonlinear minimizations. It uses non-linear least squares fitting. The fractional-order model in (4) that needs to be estimated is assumed to be mathematically defined as a function f(x) of its unknown parameters included in vector x. An objective function E(x) is defined as the error between the $y_{exp}$ and $y_{sim}$ signals, as mentioned in Section 2.3. The parameters in vector x need to be estimated such that the objective function is minimized. Maximum and minimum values for x are specified as constraints. In the TRR approach, the function f(x) is approximated with a quadratic function q(x) in a neighborhood N around the current point x and defined as the trust region. A trial step s is computed by minimizing the area N and evaluating the function f(x + s) with respect to f(x). Depending on the result, the trust region is kept unchanged or minimized in the search for the minimum [49] The algorithm stops when a minimum tolerance value is reached. In this study, the built in Matlab^® function lsqnonlin is used to implement the TRR algorithm for estimation of the parameters in (4).

3. Results

3.1. The Fractional-Order Impedance Mode

Previous research [50] has shown that the skin impedance, measured by the AnspecPro or the Medstorm devices, behaves as a dynamic system responding to noxious stimuli (acting as disturbances) or opioid drugs (such as Remifentanil, acting as the manipulated input). System identification requires persistent excitation of the dynamic process to be mathematically modeled, a feature built within the AnspecPro device. Unlike the Medstorm device, which uses data collected at a single frequency, in the case of the AnspecPro device, a multisinusoidal signal with an amplitude A, phase $\varphi$ and linearly distributed frequencies over $\omega =2\pi \left[100:50:1500\right]$ rad/s is used as an input [27]:

$$u\left(t\right)={\sum }_{k=1}^{29}{A}_k\sin \left({\omega }_{k}t\right)+{\varphi }_k$$

(7)

The amplitude at each of the k frequencies A_k and the corresponding phase ${\varphi }_k$ is chosen according to the protocol defined in [27,50] The multisine signal in (7) allows for an improved characterization of the frequency response of the skin impedance over a wider range of frequency and thus enables the use of system identification techniques. The AnspecPRO monitoring device used to collect data for analgesia modeling measures the skin impedance with a sampling period T_s = 5 s. The collected data consist of a vector of complex numbers denoting the skin impedance at each of the 29 frequencies in $\omega$ and each sampling period. These data are processed and it is possible to obtain in real-time the frequency response of the skin impedance Z( $j{\omega }_k$) at each frequency point and every sampling period. In this section, a mathematical model for the skin impedance is obtained based on the frequency response Z( $j{\omega }_k)$.

Several models for skin impedance have been proposed [27,50,51] in the absence of opioid drug and for healthy awake individuals. These models were derived based on physiological insight and are used as a starting point for those developed in this paper. Existing models [27,50,51] follow a recurrent pattern, a mathematical tool which can mimic the diffusion phenomena occurring in biological tissues [26]. The recurrent pattern can be mathematically represented as an integer-order transfer function with interlacing pole-zero pairs. Such a simplified model is represented as follows:

$$Z_{REC}\left(s\right)=K{\displaystyle \frac{\left(s^2+b_{1}s+b_2\right)\left(s^2+b_{3}s+b_4\right)\left(s^2+b_{5}s+6\right)}{\left(s^2+a_{1}s+a_2\right)\left(s^2+a_{3}s+a_4\right)\left(s^2+a_{5}s+a_6\right)}}$$

(8)

In previous research, a lumped fractional-order impedance model (FOIM) has been proposed for capturing the skin impedance [50]:

$$Z_{FOIM}\left(s\right)=R+{\displaystyle \frac{K}{1+{\left(\frac{s}{p}\right)}^{\beta }}}$$

(9)

where s is the Laplace variable, ${\left({\displaystyle \frac{1}{p}}\right)}^{\beta }$ is the pseudo-capacitance of the constant phase element and R and K are the gains at infinite frequencies and at low versus high frequencies, respectively. Fractional-order models in general and those in (9), in particular, are approximated in a certain frequency range as transfer functions of a higher order with interlacing poles and zeros. As such, they are similar to the recurrent models in (8). The justification of using fractional-order models is three-fold. First, by using fractional-order models as in (9), a more compact transfer function with fewer parameters is obtained compared to (8). Secondly, fractional-order models as in (9) have been proposed as a more suitable and accurate alternative to integer-order models [15,16,17,18,19,20,21,22,23,24,25,26,27,50,51] in biological phenomena. Third, research has shown that in drug diffusion problems, such as that of Remifentanil in the human body, the phenomena that occur are better captured with non-integer models [26]. All these justify the selection of a fractional-order (FO) transfer function to model the skin impedance.

In what follows, the measured experimental data for one patient are considered to be the frequency response Z( $j{\omega }_k$) at each of the 29 frequencies in the range $\omega =2\pi \left[100:50:1500\right]$ rad/s. These data were collected over a period of time, as indicated in Table 1, corresponding to M samples, with the sampling period T_s, thus resulting in a different $Z_i$( $j{\omega }_k$) at each sampling period i = 1,…,M. In this section, an averaged fractional-order impedance model was firstly derived. To estimate the averaged fractional-order impedance model denoted as $\overline{Z}$( $j{\omega }_k$), for each frequency, the measured frequency response was averaged over time: $\overline{Z}$  $(j{\omega }_k)={\displaystyle \frac{\sum _{i=1}^M{Z}_i\left(j{\omega }_k\right)}{M}}$. Various fractional orders were considered, and the results evaluated in terms of the cost function in (6) for all frequencies in the range of $2\pi \left[100:50:1500\right]$ rad/s.

To abide by low-complexity models, commensurate fractional-order models were considered in the following form:

$$Z_{FO}\left(s\right)={\displaystyle \frac{b_4{s}^{4q}+b_3{s}^{3q}+b_2{s}^{2q}+b_1{s}^q+b_0}{a_4{s}^{4q}+a_3{s}^{3q}+a_2{s}^{2q}+a_1{s}^q+a_0}}$$

(10)

where b_j and a_j, with j = 1, 2, 3, 4 are coefficients that need to be estimated and $q$ is the commensurate fractional order. A commensurate order larger than four was not considered suitable due to the limited availability of clinical data for validation in a frequency range of 29 values. At the same time, models exhibiting non-minimum phase dynamics were disregarded due to mismatch of physiological properties. Both pole-zero commensurate fractional-order models and all pole-only commensurate fractional-order models were considered.

The drawback of the NM algorithm consists of the dependence of the search performance upon the initial points. To achieve good optimization results, it is crucial to provide proper initialization. Several combinations of the initial q parameters were tested by trial and error. Using frequency domain identification tools, as described in Section 2, unstable transfer functions were obtained for commensurate fractional orders q < 1. Stable transfer functions were obtained for q in the range of [1, 2]. To estimate the fractional-order models, an initial starting value q* = 1.5 was used in the optimization method for all patients. For simplicity, only the results for patient 70 are presented here, although the same conclusions can be drawn by analyzing the results for different patients.

Table 2. Cost function and model complexity for impedance modeling.

Number of poles na and zeros nb of the FO models	na = 4 nb = 4	na = 3 nb = 3	na = 2 nb = 2	na = 4 nb = 0	na = 3 nb = 0	na = 2 nb = 0
Cost function J as defined in (6)	5.1544	6.2547	5.6648	8.0844	7.3194	6.4329

summarizes the cost functions obtained with the various fractional-order (FO) models.

Based on the results in Table 2, the minimum value of the cost function J = 5.1544 is obtained for fourth-order models with n_b = 4 and n_a = 4. To reduce the complexity of the impedance model, a simplified fractional-order transfer function was preferred with n_a = 2 and n_b = 0 and a corresponding J = 6.4329, 25% larger compared to the more complex fourth-order model. The final fractional-order impedance model has the following form:

$$Z_{FO}\left(s\right)={\displaystyle \frac{b_0}{a_2{s}^{2q}+a_1{s}^q+a_0}}$$

(11)

The resulting fractional-order impedance model for patient 70 is given in the next equation and the validation with clinical data is indicated in Figure 2a:

$$Z_{FO\_P70}\left(s\right)={\displaystyle \frac{6924}{2.87{·10}^{-12}{s}^{3.148}+7.36·{{10}^{-6}s}^{1.574}+1}}$$

(12)

Figure 2 portrays the validation of the FO impedance models, similar to (12), for three other patients included in the database. Notice that the magnitude approximation is more accurate. The phase is highly corrupted by noisy data, especially at high frequencies. To be able to apply the system identification algorithm from Section 2, the clinical (initial) data were smoothened and filtered using the Matlab^® programming tool. The plots in Figure 2, however, depict the raw data.

3.2. The Fractional-Order Analgesia Model

The ‘pain index’ (as a measure of nociception during anesthesia) is defined in what follows as a means to evaluate changes in the measured skin impedance due to variations in the progression rate of the opioid drug Remifentanil. Figure 3 shows the drug progression rate over a period of 80 s for patient 70. Skin impedance was measured over this period of time at every Ts = 5 s. To evaluate changes in the measured skin impedance due to variation in the Remifentanil input, a fractional-order model as in (11) is estimated at each sampling period. For the NM algorithm, the initial starting guess for the commensurate fractional order in (11) is considered to be q* = 1.5, similarly to Section 3.1. At each sample, the previous model is used as an initial guess.

Table 3. Commensurate fractional order at each sample for patient 70.

Sample	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16
q	1.63	1.62	1.56	1.53	1.53	1.54	1.55	1.54	1.52	1.54	1.56	1.58	1.58	1.57	1.55	1.54

includes the commensurate fractional orders of the impedance model obtained using the optimization routine at each sampling instant during the 80 s interval. After this point, the administered Remifentanil drug remains constant and the commensurate fractional order q does not exhibit any significant variations. From a control engineering perspective, it is assumed that the steady state has been reached.

A correlation between the Remifentanil input in Figure 3b and the commensurate fractional order q in Table 3 and Figure 3a shows that an increase in the drug rate leads to changes in q. Similar results are obtained for other patients under study, as indicated in Figure 4, Figure 5 and Figure 6. As such, the commensurate fractional order of the impedance q will be further considered the ‘pain index’ and the output q(t) of the analgesia model. The input to this model is the Remifentanil progression rate, denoted as Remi(t).

As shown in Figure 4, Figure 5 and Figure 6, there is a lack of data both in terms of quantity and quality, with poor information signals available from the patient. As such, pharmacokinetic–pharmacodynamic (PK-PD) compartmental models are difficult to be developed. The compartmental modeling approach assumes the drug is intravenously applied to the blood compartment and diffused in the muscle and fat compartments before reaching a hypothetical effect site compartment. This final compartment represents the effect of the drug on the nervous system [50]. Differential equations are used to model the pharmacokinetic compartments [51,52]:

$$\begin{array}{l}{\dot{x}}_1\left(t\right)=-\left(k_{10}+k_{12}+k_{13}\right)x_1\left(t\right)+k_{21}{x}_2\left(t\right)+k_{31}{x}_3\left(t\right)\\ {\dot{x}}_2\left(t\right)=k_{12}{x}_1\left(t\right)-k_{21}{x}_2\left(t\right)\\ {\dot{x}}_3\left(t\right)=k_{13}{x}_1\left(t\right)-k_{31}{x}_3\left(t\right)\end{array}$$

(13)

where x_i(t) are concentrations in the blood (central), muscle, and fat compartments and k_ij are the drug transfer rates from the i-th to the j-th compartment. From a control engineering perspective, the muscle compartment ${\dot{x}}_2\left(t\right)$ is a fast-acting compartment, which is characterized in a simplified approach as a small, parasitic time constant compared to the central compartment. ${\dot{x}}_1\left(t\right)$. In a controller design approach, this compartment can be neglected, as its effect is minimal. The central compartment is responsible for the rise time and rapid onset of the drug, whereas the fat compartment represented by ${\dot{x}}_3\left(t\right)$ is characterized by very slow dynamics, which lead to a sluggish effect close to the steady-state value. From this point of view, in control engineering, the muscle and fat compartments bring little information on the overall dynamics, with the main dynamics generated by the blood compartment. Thus, the integer-order model in (13) can be simplified to a first-order transfer function, with time constant ${\tau }_1$.

The effect site compartment is modeled by the following differential equation:

$${\dot{x}}_e\left(t\right)={-k}_{0e}{x}_e\left(t\right)-k_{1e}{x}_1\left(t\right)$$

(14)

where $k_{e0}$ is the drug metabolic rate and $k_{1e}$ is the drug transfer rate from the central compartment to the effect site, which can be written as a first-order transfer function:

$${\displaystyle \frac{x_e\left(t\right)}{x_1\left(t\right)}}=-{\displaystyle \frac{\frac{k_{1e}}{k_{0e}}}{\frac{1}{k_{0e}}s+1}}$$

(15)

The dynamics of the integer-order model in (15) are mainly characterized by the time constant ${\tau }_2={\displaystyle \frac{1}{k_{0e}}}$.

The pharmacodynamic part of the model is mathematically represented by a nonlinear Hill curve, which is ultimately reduced to a gain in steady state. This models the drug diffusion from effect site to their final effect. With all of the above considerations, the final simplified integer-order transfer function for DoA would be a damped second-order system of the form:

$$H_{DoA}\left(s\right)={\displaystyle \frac{k}{({\tau }_{1}s+1)({\tau }_{2}s+1)}}$$

(16)

As indicated in [26], for drug diffusion in the human body, a much better approach is a fractional-order model. In what follows, such a fractional-order model is proposed with two stable overdamped poles, no zeros and time delays:

$$H_{DoA}\left(s\right)={\displaystyle \frac{k}{a_1{s}^{{\alpha }_1}+a_2{s}^{{\alpha }_2}+1}}={\displaystyle \frac{q\left(s\right)}{Remi\left(s\right)}}$$

(17)

where k is the steady-state gain between the drug progression rate Remi(s) as the input and the nociception index q(s) as the output, with a₁ and a₂ as real coefficients, where ${\alpha }_1$ and ${\alpha }_2$ are the fractional orders. The fractional-order model in (17) considers the “pain index” (the fractional order of the impedance model q) as the output to be controlled, while the input is the Remifentanil progression rate. The TRR algorithm presented in Section 2 is used to estimate the parameters of (17), resulting in the following fractional-order model for analgesia for patient 70:

$$H_{DoA70}\left(s\right)={\displaystyle \frac{-0.24}{1087.9s^{1.84}+99.38s^{0.89}+1}}={\displaystyle \frac{q\left(s\right)}{Remi\left(s\right)}}$$

(18)

Figure 3a presents the ‘pain index’ q as indicated in Table 3, as well as the validation of the fractional-order model in (18). Similar results were obtained for the other patients included in this study, as indicated in

Table 4. Estimated parameters of the fractional-order analgesia models for all patients under study.

Patient No.	K	a1	a2	${\mathit{\alpha }}_1$	${\mathit{\alpha }}_2$	ISE	GOF [%]
3	−0.636	14,518	105	1.9	0.89	0.006	35.86
10	−0.33	21,699	753	2	1.16	0.049	26.81
14	−0.245	12,650	285	1.95	1.09	0.023	23.98
21	−0.84	8500	245	1.86	0.88	0.135	12.77
24	−0.525	22,445	335	1.88	0.94	0.012	53.52
25	−0.7	11,003	130	1.79	0.87	0.025	18.48
33	−0.45	7800	280	1.93	1.06	0.018	24.32
35	−0.715	3166.5	124.21	1.83	0.84	0.186	12.59
38	−0.8	3961.4	309.46	2	1.15	0.369	5.64
39	−0.11	1380.6	360.43	1.81	1.1	0.002	2.39
41	−0.33	12,361	785.42	2	1.26	0.192	4.86
48	−0.18	9500	470	1.88	0.88	0.011	23.20
54	−0.09	16,742	287.53	1.86	0.99	0.001	16.30
56	−0.2	9171.8	172.19	1.88	0.87	0.043	4.24
59	−0.35	4063	22.16	1.95	0.69	0.019	23.66
63	−1.38	5575.2	450	1.67	0.86	0.1323	25.66
64	−0.9	10,400	340	1.93	0.87	0.024	17.11
66	−0.104	7573	50.534	1.99	0.78	0.018	21.88
70	−0.24	1087.9	99.38	1.84	0.89	0.011	24.95

.

The performance measure used to evaluate the accuracy of the models is represented by the Goodness Of Fit (GOF), defined as follows:

$$GOF=100·\left(1-{\displaystyle \frac{\sqrt{\sum _{i=1}^n{(q_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(i\right)-q_{sim}(i\left)\right)}^2}}{\sqrt{\sum _{i=1}^n{(q_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(i\right)-q_{mean}(i\left)\right)}^2}}}\right) [\%]$$

(19)

where $q_{sim}$ is the simulated model output from (17), $q_{exp}$ is the commensurate fractional orders of the impedance model computed at each sample i and $q_{mean}$ is their mean value. Figure 4, Figure 5 and Figure 6 show the validation of the fractional-order models (17) for three other patients included in the study. In these cases, the best GOF was obtained, according to Table 4. The performance measure GOF in Table 4 is similar to that obtained in [12] for different analgesia models. The computed Integral of Squared Error (ISE) is also indicated in Table 4.

Remark No. 1: The fractional-order model in (17) is a measure of reference of analgesia related to the level of Remifentanil, considering a zero-nociception setting. The resulting model is validated for increasing levels of Remifentanil, according to the Minto model. As such, the proposed fractional-order model is a reference model for further investigation related to the level of nociception.

Remark No. 2: The signal can be influenced by the administered hypnotic (Propofol). However, during the period of time in which the data were used to derive the model, the level of hypnotics is stable and considered to be constant in every patient.

Remark No. 3: No fractional-order transfer function models for DoA have been proposed so far. Previous research reported in [15] mentions integer-order models for DoA of the form (16), where a different nociception index is considered the output signal. The settling times reported in [15] for the 19 patients are similar to those obtained in this manuscript for the developed fractional-order models. The gains, however, differ due to the different types of nociception index considered as a mean of evaluation the DoA.

3.3. Comparison with Integer-Order Analgesia Model

To validate the advantage of using fractional-order models, an integer-order one is estimated for patient 24. In this case, the best fit was obtained according to the GOF in Table 4. The obtained fractional-order model is given as follows:

$$H_{DoA24}\left(s\right)={\displaystyle \frac{-0.525}{22445s^{1.88}+335s^{0.94}+1}}$$

(20)

An integer-order model was estimated using the same TRR algorithm and the same constraints, leading to:

$$H_{DoA24\_IO}\left(s\right)={\displaystyle \frac{-0.525}{22200s^2+430s+1}}$$

(21)

The validation of the two models with respect to the computed commensurate order q is given in Figure 7. The GOF for the fractional-order model is computed as 53.52%, whereas the GOF for the integer-order model is 39.94%. This shows that the fractional-order model is more accurate, yielding an improvement of approximately 25% compared to the integer-order model. Considering the quality and quantity of the experimental data, such an improvement is desirable for a subsequent design of the controllers.

3.4. A Simplified Preliminary Fractional-Order Controller Design for Analgesia

A Proportional Integral–Fractional-Order Proportional Derivative (PI-FOPD) controller is designed for analgesia models as indicated in (17). Fractional-order controllers are generalizations of their integer-order counterparts [32,33,34]. As such, they come with increased flexibility and better closed-loop performance, especially in terms of robustness [32,33,34,38,39,40,41]. The challenges of controlling the DoA, mathematically represented by the fractional-order model in (17), reside in the parameter uncertainty and great variability, as shown in Table 4. These parameter variations are due to patient variability [53,54]. A controller designed for a nominal patient model needs to be robust to tackle these modeling uncertainties. A PI controller is first designed to meet the requirements of zero steady-state errors. The FOPD controller is designed in series with the PI controller to enhance the stability of the control system, and reduce the overshoot and the settling time. The transfer function of the dedicated PI-FOPD controller is given as follows:

$$H_C\left(s\right)=k_p\left(1+{\displaystyle \frac{k_i}{s}}\right)\left(1+k_d{s}^{\mu }\right)$$

(22)

where k_p, k_i and k_d are the proportional, integral and derivative gains, and $\mu \in \left(\mathrm{0,1}\right)$ is the fractional order of differentiation. To tune the parameters of the controller, a set of three performance specifications is addressed. The first one refers to the settling time requirement, addressed via a desired gain crossover frequency ${\omega }_c$:

$$\left|H_{OL}\left(j{\omega }_c\right)\right|=\left|H_C\left(j{\omega }_c\right)\right|\left|H_{DoA}\left(j{\omega }_c\right)\right|$$

(23)

where $H_{DoA}\left(j{\omega }_c\right)$ is the frequency response of the patient model in (17) and $H_{OL}\left(j{\omega }_c\right)$ stands for the loop frequency response. The larger the ${\omega }_c$ is, the faster the settling time would be. The second performance criteria refer to the overshoot requirement. This is mathematically addressed in the frequency domain using the phase margin (PM) equation:

$$\angle H_{OL}\left(j{\omega }_c\right)=\angle H_C\left(j{\omega }_c\right)+\angle H_{DoA}\left(j{\omega }_c\right)=-\pi +PM$$

(24)

Theoretically, the larger the PM is, the smaller the overshoot would be. A large PM typically means a better stability of the closed-loop system as well. Finally, to address patient variability in the model (in terms of gain variations), the iso-damping condition is specified as a third design criterion:

$${\left.{\displaystyle \frac{d\angle H_{OL}\left(j\omega \right)}{d\omega }}\right|}_{\omega ={\omega }_c}={\left.{\displaystyle \frac{d\angle H_C\left(j\omega \right)}{d\omega }}\right|}_{\omega ={\omega }_c}+{\left.{\displaystyle \frac{d\angle H_{DoA}\left(j\omega \right)}{d\omega }}\right|}_{\omega ={\omega }_c}=0$$

(25)

To solve the system of three Equations (23)–(25), the Matlab^® “fsolve” function is used, which implements the Levenberg–Marquardt algorithm [55,56]. The nominal patient model considered in this section for the design of the controller is patient 24 with the fractional-order model given in (20). The imposed gain crossover frequency is ${\omega }_c=0.08$ and PM = 70°. For the nominal patient model in (20), the parameters of the PI-FOPD controller were determined as follows: $k_p=34.17$, $k_i=0.034$, $k_d=51.8$ and $\mu =0.63$. The Bode diagram of the loop transfer function for the designed controller and the nominal model in (20) is given in Figure 8, clearly showing that the design specifications were met.

Figure 9a,b shows the closed-loop simulation results for the pain index q for all patients in Table 1, as well as the corresponding Remifentanil progression rate. The amplitude of the pain index q was normalized for a more accurate comparison. The results in terms of overshoot and settling time are indicated in

Table 5. Performance indexes for the closed-loop simulation results obtained with the proposed fractional-order controller.

Patient No.	Overshoot	Settling Time (s)
24	5%	61
3	10%	43
10	7%	68
14	8%	63
21	0%	68
25	4%	33
33	5%	36
35	0%	30
38	6%	28
39	0%	>300
41	8%	45
48	0%	>300
54	0%	105
56	0%	168
59	10%	28
63	0%	97
64	3%	94
66	8%	>300
70	0%	>300

. The results correspond to the pharmacological properties of Remifentanil, with a rapid onset of action of approximately 1 min [57]. In some cases, the settling time is significantly less than 1 min and well above 5 min (more than 300 s), as indicated in Table 5. The overshoot for the nominal patient P24 is 5% and remains small (maximum 10%) for all other patients included in this case study. The Remifentanil maximum drug rate indicated in Figure 9b corresponds to experimental data, as seen in Figure 4, Figure 5 and Figure 6, for example. However, the drug rate is maintained at a higher level (similarly to [54] for a longer period of time, which leads to a faster settling time for some patients (less than 60 s, as indicated in Table 5). The steady-state value of Remifentanil drug rate is approximately 0.4 $\mathsf{\mu }\mathrm{g}/\mathrm{k}\mathrm{g}/\mathrm{m}\mathrm{i}\mathrm{n}$, corresponding to clinical data [58,59]

4. Conclusions

Extensive research is now being conducted on automatic control in general anesthesia, with an emphasis not just on hypnosis modeling and control but on a comprehensive strategy that also addresses hemodynamic parameters, analgesia, etc. Due to a lack of specialized monitoring devices, little research has been carried out on nociception and how opioid drugs affect it. As a result of some novel devices measuring skin impedance or conductance, the idea of modeling the Depth of Analgesia is starting to gain momentum. Various impedance models have been proposed, but most of them have been developed in the absence of opioid drug. A recent paper has proposed the use of a pain index to evaluate nociception in patients relative to the administered opioid drug rates. In this paper, a new pain index is defined and validated on clinical available data from patients undergoing surgery. Then, this pain index is considered a reference measure for the Depth of Analgesia and its dynamics are correlated to those of the administered Remifentanil in the absence of nociception. The reference model can be considered a calibration model to level the intensity of nociception. Fractional-order models are used to mathematically describe the correlation between the pain index and Remifentanil infusion rates. Standard optimization routines are selected to estimate the parameters of these fractional-order models. Clinical data from 19 patients are used to determine individualized fractional-order models. One of these models is randomly selected as the nominal model and used in the design of a fractional-order controller. The simplified controller manages to achieve good closed-loop performance corresponding to clinical requirements. Improved closed-loop performance despite patient variability could be obtained with a controller tackling more performance specifications. Further research to enhance the closed-loop robustness of the designed controller is required, as well as a decision on the micropumps and the practical implementation of the control system.