Improved Individualized Patient-Oriented Depth-of-Hypnosis Measurement Based on Bispectral Index

Total intravenous anesthesia is an anesthesiologic technique where all substances are injected intravenously. The main task of the anesthesiologist is to assess the depth of anesthesia, or, more specifically, the depth of hypnosis (DoH), and accordingly adjust the dose of intravenous anesthetic agents. However, it is not possible to directly measure the anesthetic agent concentrations or the DoH, so the anesthesiologist must rely on various vital signs and EEG-based measurements, such as the bispectral (BIS) index. The ability to better measure DoH is directly applicable in clinical practice—it improves the anesthesiologist’s assessment of the patient state regarding anesthetic agent concentrations and, consequently, the effects, as well as provides the basis for closed-loop control algorithms. This article introduces a novel structure for modeling DoH, which employs a residual dynamic model. The improved model can take into account the patient’s individual sensitivity to the anesthetic agent, which is not the case when using the available population-data-based models. The improved model was tested using real clinical data. The results show that the predictions of the BIS-index trajectory were improved considerably. The proposed model thus seems to provide a good basis for a more patient-oriented individualized assessment of DoH, which should lead to better administration methods that will relieve the anesthesiologist’s workload and will benefit the patient by providing improved safety, individualized treatment, and, thus, alleviation of possible adverse effects during and after surgery.


Introduction
Adequate general anesthesia (GA) is a prerequisite in surgeries as well as in various other medical procedures. The anesthesiologist must take care of three main aspects of the patient state during the procedure: besides considering the vital signs, they must administer substances that keep the patient deeply unconscious, prevent the patient from feeling pain, and keep the patient's muscles adequately relaxed. Furthermore, the patient must not be aware of, or even remember, what was happening during GA. Therefore, it is essential to properly administer the needed substances during the medical procedure.
GA and the dynamic response of the patient's body to anesthetic substances can be regarded as quite complex dynamic systems. Various pharmacokinetic (PK) and pharmacodynamic (PD) mechanisms take place inside the patient's body; however, it is unfortunately generally not yet possible to claim that these PK and PD systems have been completely studied or adequately modeled. Furthermore, the anesthetic depth is a controversial concept, which involves the three main aspects of hypnosis, analgesia, and muscle relaxation. In addition, amnesia must also be ensured. However, anesthetic depth cannot be directly measured.
The main job of the anesthesiologist is to monitor the patient's vital functions and properly maintain the functions of the patient's organs. Anesthetic agents are administered in various ways into the patient's body so as to achieve adequate GA: the two most commonly used methods in clinical practice are the intravenous administration induction of the anesthetic agent, i.e., injection of the anesthetic agent into a patient's vein, or the inhalatory administration of anesthetic agent, in which the substance is induced by the patient inhaling a prepared breathing mixture. If the anesthetic agents are injected intravenously, the anesthesiologic technique is known as total intravenous anesthesia (TIVA) [1,2]. In this article, we focus on TIVA administration exclusively.
During the medical procedure, the anesthesiologist aims to maintain the appropriate depth of hypnosis (DoH) by adjusting the dosage of the anesthetic agent. Clearly, the PK and PD of the anesthetic agent and the type of procedure must be taken into account, e.g., long-term sedation in an intensive care unit requires deeper DoH than GA during surgery. Too-deep anesthesia is manifested as a drop in blood pressure level and heart rate frequency as well as slow post-operative awakening of the patient from GA [3]. On the other hand, the opposite is also true. Moreover, inadequate depth of anesthesia results in the activation of sympathetic nerves, or, in the most unlikely event, with the patient awakening, which must be avoided at all costs [4]. In modern clinical practice, DoH is determined by assessing the relevant clinical signs (iris, sweating, movements), by interpreting hemodynamic measurements [5] and by estimating the DoH from EEG signals [6]. The latter is made possible by several established measurement systems, such as BIS index, NeuroSense, Narcotrend, Entropy (Scale and Response), WAV CNS , and Patient State Index.
BIS index is a noninvasive measurement method. A BIS monitor is connected to electrodes on the patient's head and the bispectral index is calculated from the measured EEG signals. The BIS monitor provides a single dimensionless number, which ranges from 0 (equivalent to EEG silence) to 100. A BIS value between 40 and 60 indicates an appropriate level for GA, whereas a value below 40 is appropriate for long-term sedation due to head injuries. The reference can thus be set to the applicable value; the manner and speed of approaching the reference value depend on the specific characteristics of the procedure and the pharmacokinetics and pharmacodynamics of the substance in the patient's body. The BIS value can, therefore, be considered as a representation of the DoH, although some papers question the relevance of BIS measurement for representing hypnosis as a graded state during surgery [7,8]. The authors in [7] pointed out the poor correlation between BIS and serum concentrations of propofol, which calls into question the relevance of the BIS measurement in representing hypnosis as a graded state during surgery. When they included all values during anesthesia, there was a significant correlation between BIS and serum propofol, but the correlation was poor and disappeared when the outliers were removed. Their results do not seem to support the idea that BIS represents a continuous measure of DoH, assuming that elevated serum concentrations of propofol correspond to deeper hypnotic levels. However, because they have too few data points to test the correlation between serum concentration and BIS in individual patients, the reason could be that there are interindividual differences in sensitivity to propofol. They noted that a correlation could exist in a single patient but is not evident in data from many subjects.
In the literature, there are several approaches to modeling the effect of propofol, which is a hypnotic anesthetic agent. For this purpose, a number of PK and PD models have been developed, e.g., [9][10][11][12]. The models of propofol response define the general structure of the dynamic system, whereas the particular parameters depend on the individual patient's characteristics, such as gender, age, height, weight, etc., as well as the particular patient's individual sensitivity to the effects of propofol and their ability to have propofol eliminated from the body. Certain infusion pumps employ the PK models to enable target controlled infusion (TCI), where the pump sets the proper flow of the medication with regard to the model [13].
In the last years, an emerging paradigm in medicine seems to have grasped the idea of personalized medicine [14]. In the field of drug delivery systems, this includes modeling, control, analysis, and pharmacological studies, as well as development of novel medical devices and conducting of targeted clinical trials. In this regard, the systematic employment of dynamic-system analysis along with control theory offers a wide range of application opportunities in the medical domain [15][16][17][18][19][20]. Despite the fact that for TCI various PK models can be implemented, all having their own advantages and drawbacks, the Marsh [21] and Schnieder [22,23] models are mainly used in clinical practice. However, these models often do not reflect the actual dynamic properties that depend on individual sensitivity to the substance of the particular patient, which is generally not considered. For instance, the authors in [24] found that titration of propofol based on BIS monitoring allows a reduction in drug consumption that is associated with a similar decrease in propofol plasma concentrations compared with TCI. An inverse relationship between cardiac output and plasma propofol concentration was reported in [25]. In addition, the authors in [26] showed that remifentanil plasma concentrations during remifentanil and propofol anesthesia are influenced by cardiac output in a similar manner to propofol, although the metabolic sites are different. Because cardiac output is known not to be constant during anesthesia, this also significantly affects plasma concentrations and, thus, the effects of propofol and remifentanil. The authors concluded that actual blood concentrations of remifentanil and propofol may differ significantly from expected concentrations, especially when cardiac output is low. Cardiac output is therefore an important factor whose influence should be investigated in PK and PD models. This again supports the idea that individualized PK and PD modeling could improve anesthetic delivery methods. Therefore, approaches using population-based-data models often cannot ensure optimal performance, especially when treating a patient with a particularly considerable discrepancy from the population-databased models. A population-data-based approach cannot yield universally usable models, even if a very large number of patients would be taken into account and despite the fact that only relatively healthy patients are considered.
The dynamic mathematical models can be directly used in clinical practice for assessing DoH-related variables by implementing them in soft sensors [27] or state observers [28], which would improve the anesthesiologist's assessment of the patient state regarding anesthetic agent concentrations and, consequently, the effects, hence benefiting the patients in the long run. In addition, the models can enable in silico testing of various potentially clinically implementable closed-loop control algorithms, e.g., PID-controller-based approaches [29][30][31]. Furthermore, they represent the basis for a number of advanced control approaches, such asrobust control [32][33][34], model-predictive control [35,36], fuzzy-rulebased decision system [37], event-based control [38], etc. Despite difficulties in objective pain measurement [18,39], a number of articles also consider the inherent MIMO (or MISO) nature of the controlled system, which is due to drug interactions, especially when analgesics (such as remeifentanil) are considered [40][41][42].
The article introduces a novel structure for modeling DoH. The modeling framework results in an improved individualized assessment of DoH, which is reflected in better predictions of the related anesthetic agent concentrations as well as the measured BIS signal. The article is structured as follows. First, we introduce the basic three-compartmental model and the upgraded population-data-based model, namely, the PK and PD parts. Next, we introduce the improved model structure and present the residual dynamic model add-on. In Section 4, the identification procedure is presented, including the identification signals, filtering, and parameter estimation. This is further extended to online recursive parameter estimation. Next, we validate the improved model by comparing the simulated data for the particular patient to real clinically acquired data, which were logged during a surgery in a real clinical setting. We end the article with some concluding remarks.

Propofol Pharmacokinetic and Pharmacodynamic Modeling
The pharmacokinetics of the basic population-data-based model is described by the Schnider model [22,23]. A well-established three-compartmental model structure is used as the basis for dynamic relations. For details, see [43,44]. The internal dynamics of the model can be formulated using Equations (1)-(3).
In Equations (1)-(3), the variables x 1 , x 2 , and x 3 represent the amount of the drug in compartments V 1 , V 2 , and V 3 , respectively. The infusion flow rate is denoted as φ. As noted above, the parameters k 12 , k 21 , k 13 , and k 31 represent the partition coefficients that determine the speed at which the drug moves from one particular compartment to another. Finally, k 10 is the rate of elimination of the drug from the patient's body. Note that the concentration in the central compartment is often referred to as plasmatic concentration.
The effect site for the drug propofol is basically the central nervous system. The effect site is thus part of the central compartment, but the effect of the drug is subject to some dynamics with regard to the (theoretical) concentration in the central compartment [44]. Therefore, a first-order model was used to describe the effect-site concentration dynamics, as given in Equation (4).
The effect-site concentration of propofol is considered as the main influence on the DoH. Despite acknowledging that DoH is a multivariable and a not very easy to grasp concept, involving deep unconsciousness, analgesia, amnesia, and muscle relaxation, we want to keep the model as simple as possible, yet not too simple for our requirements. Therefore, we first assume that BIS index is an adequate measure for DoH, despite the fact that the assumption might be debatable [8,45,46]. Furthermore, we consider an SISO model, where the input represents propofol inflow, and the output represents the value of the BIS index.
Despite the fact that the PD effect mechanism has not been fully studied yet, in the literature, a sigmoid E max model based on the general Hill equation [47] is usually considered, as given in Equation (5).
In Equation (5), BIS 0 denotes the characteristic value for a fully awake patient, BIS min stands for the minimum value of BIS index, and γ is the parameter that defines the nonlinear shape of the response curve.
Therefore, the combined model can be regarded structurally as a Wiener nonlinear model.
The parameters of the phamacokinetic model are taken from [22,23,48]. The values are given in Table 1. Note that the values of model parameters depend on the patient's age, weight, height, and gender. Parameter LBM is calculated as given in Equation (6).
The parameters of the BIS-index-effect output submodel (see Equation (5)) are set as suggested in [49] and are given in Table 2.

Residual Model Introduction-Improvement of the Population-Data-Based Model
Population-data-based mathematical models, such as the one described in the previous section, are regularly used in clinical practice, namely, for TCI propofol infusion, which has become a standard approach in administration of the anesthetic agent. In spite of this, we need to consider the fact that the available models are derived from population-based measurements, both regarding the pharmacokinetic and the pharmacodynamic part.
Prior to a medical procedure, the implemented mathematical model is always tuned to the particular patient's properties (e.g., age, weight, height, and gender). Despite this, the model is still based on a broader population sample that was involved in the measureddata gathering. Hence, it is impossible for such a model to take into account specific interpatient variabilities. Therefore, the basic population-data-based model can broadly predict DoH, but it can exhibit severe discrepancies from the actual BIS-indicated DoH of a particular patient, especially when the particular patient's sensitivity to the anesthetic agent is considerably different from the dynamics assumed in the mathematical model.
In order to improve the population-data-based model accuracy by taking into account the specific patient's individual dynamic properties, we propose an extension to the basic model structure. The whole model is then structured as follows: besides the population-data-based model, we introduce an additional residual model, which is intended to mathematically describe the dynamic discrepancy of the patient's particular sensitivity to anesthetic-agent infusion. The discrepancy signal represents the ideal residual model output. The improved model structure, with the additional residual dynamic addon, is shown in Figure 1. The output of the improved model BIS sim is, therefore, calculated as given in Equation (7).

PDB model
Residual model Here, BIS PDB denotes the basic population-data-based model output, and y res is the residual model output.
The goal of introducing the residual model is thus to improve the combined model output. This means that the combined model includes the population-data-based model as the basis for modeling the patient dynamic response to propofol. In addition, the residual model considers the particular patient's individual dynamic response to propofol, namely, the individual patient discrepancies from the population-data-based model, thus enabling a more accurate assessment of DoH expressed by BIS index. Note that the following conditions prevent a patient from being included in the study: patients with poor general condition (ASA > 3), patients with BMI > 35, drug addicts, patients taking psychotropic medicines or opioid analgesics (including tramadol), patients with a severe psychiatric disease or central nervous system disease (except the reason for surgery), patients with arrhythmia affecting or preventing the measurements (e.g., chronic atrial fibrillation), and patients that received benzodiazepines.

Identification Signals
The residual dynamic model is identified based on the data measured during a medical procedure treating an individual patient. The input of the residual dynamic model is the actual propofol inflow φ propo f ol , whereas its output is obtained by subtracting the actual measured BIS signal BIS data from the population-data-based model output BIS PDB , as shown in Equation (8). The latter is calculated according to the predefined population-databased PK and PD model, taking into account the individual patient's properties.
Both signals are sampled using T s = 1 s sampling rate. The whole experiment is therefore represented by two discrete-time signals in Equations (8) and (9). y res,id (z) = BIS PDB (z) − BIS data (z) (8) u res,id (z) = φ propo f ol (z)

Residual Dynamic Model Structure
The structure of the residual dynamic model should be as simple as possible, but at the same time complex enough to be able to adequately model the patient's discrepancy from the population-data-based model. Despite that higher complexity of the models and possibly evolving identification could yield favorable approximations [50][51][52], especially if the identification data are noise-free, we established that an appropriate structure for our case is a second-order model, which is structured as an affine autoregressive model with exogenous inputs. Its discrete-time formulation is given in Equation (10).
Here, a 1 , a 2 , b 1 , and c represent the parameters to be estimated. Note that t represents a particular time-instant of the model. A single step of the model T s is not necessarily equal to the data-sampling rate T s,data . In our case, it is T s = 5 · T s,data = 5 s.
The model parameters can be gathered in the parameter vector Θ as given in Equation (11).
The regressor is defined in Equation (12).
In this way, Equation (10) can be rewritten in Equation (13).

Parameter Estimation
The parameters of the residual dynamic model a 1 , a 2 , b 1 , and c are estimated from the measured data concerning the particular patient.
As the measured BIS data signal is prone to significant noise, the first step is to apply a suitable filter, which should ensure better identification results. We established that a simple first-order filter (with a filtering time-constant τ f ) is adequate. The filter can be represented by the transfer function in Equation (14). Its discrete-time equivalent is given in Equation (15).
Finally, the filtered identification signals are given in Equations (16) and (17).
y res,id, f (z) = H f (z)y res,id (z) (16) u res,id, f (z) = H f (z)u res,id (z) (17) The filtered signals y res,id, f and u res,id, f are used for estimating the parameters of the residual dynamic model Θ.
The output data vector Y contains the output variable as given in Equation (18).
The regression matrix Ψ is obtained by using the whole set of measured data, as given in Equations (19) and (20).
Here, t i represents a time instant concerning a particular identification data pair (i = 1, . . . , P).
According to Equations (18)-(20) and (13), the parameters of the residual dynamic model Θ can be obtained using the least-squares identification method, as given in Equation (21).

Online Recursive Parameter Estimation
In order to implement the proposed modeling framework as an intelligent soft sensor for assessing DoH, we must be able to acquire the improved combined model for DoH during surgery. The population-data-based submodel is based on patient data and can be derived prior to surgery. On the other hand, the particular patient's individual dynamic response to propofol, namely, the individual patient discrepancies from the populationdata-based model, can only be assessed after the measured data becomes available, i.e., either after finishing the medical procedure, as proposed in Section 4.3, which can rarely be used, or during the medical procedure as soon as new data become available. The latter approach enables the use of the improved model during the medical procedure.
The recursive parameter estimation is carried out as described below. The residual model in Equations (10) and (13) is linear in the parameters, therefore it is possible to analytically derive a least-squares estimate of the parameters. Furthermore, if the identified system parametersare expected to be time-varying, the online estimation algorithm can place more emphasis on newly acquired data and gradually discard older data. Therefore, the proposed approach implements a recursive least-squares identification with exponential forgetting [53]. The algorithm can consider a least-squares loss function for exponentially discarding older data as time passes. The model parameters are, therefore, estimated using Equations (22)- (24).
In Equations (22)-(24), P(t) denotes the covariance matrix (P(t) ∈ R n×n ), where n represents the length of the regressor.Θ(t) denotes the vector of the identified or estimated process parameters, λ denotes the forgetting factor, and I is the unity matrix, where I ∈ R n×n . The recursive parameter estimation is carried out online-in every time step t-and returns the calculated parameters of the modelΘ(t).
The forgetting factor λ is defined in Equation (25), where t λ stands for the time constant for the exponential forgetting and T s is the sampling time.
The regressive parameter estimation method can consider time-varying dynamics of the identified process; therefore, exponential forgetting can be employed. The forgetting factor has to be set between 0 < λ ≤ 1. In this manner, the data used for parameter estimation are pondered, so that the last data are pondered by factor 1, whereas the data that are k s · T s time steps old are pondered only by a factor of λ k s .
The initial covariance matrix P(0) has to be positive-definite and properly sized. If the regressive parameter estimation method is interpreted as a Kalman filter, it can be regarded as ensuring that the parameters are distributed with an initial covariance P(0) and initial mean-valuesΘ(0) [53].
In each time step t, a new estimation is calculated. The new parameter estimationsΘ(t) are recursively based on the estimations from the previous time step t − T s and the online newly measured data. Therefore, the resulting model adapts to the new measurements as soon as they are available.

Model Validation Based on Real Clinical Data
In order to validate the proposed modeling approach, we gathered real data during a TIVA medical procedure. In our case, the treated patient was a 46-year-old male who weighed 59 kg and was 175 cm tall.

Recorded Signals
The signal denoting the inflow of propofol φ propo f ol was parsed from the clinically recorded data and is shown in Figure 2. One can see that, firstly, a bolus dose is introduced in order to rapidly increase the concentration of propofol in the body. This phase is called the induction of anesthesia and results in the patient losing consciousness. Later, a suitable dose of propofol is continuously administered in order to keep the proper anesthetic depth.
The recorded BIS data signal for the particular treated patient is shown in Figure 3. The plasmatic concentration of propofol c p and effect-site concentration of propofol c e are obtained by simulating the dynamic model defined in Equations (1)-(4) based on the measured inflow signal of propofol φ propo f ol . Next, the population-data-based BIS trajectory BIS PDB is obtained according to Equation (5). Finally, it is possible to obtain the identification signals for the residual dynamic model y res,id, f and u res,id, f using Equations (8), (9), and (15)- (17). The identification signals are shown in Figure 4. Note that only the relevant interval of the measured data is considered.

Identification Results
The final resulting parameters of the residual dynamic model are given in Equation (26).
When conducting online recursive parameter estimation, the parameters converge to the values in Equation (26). Before starting, the initial values for the algorithm are set as given in Equations (27) As the recursive parameter estimation is carried out, the parameters a 1 , a 2 , b 1 , and c change their values along the trajectories shown in Figure 5.

Model Validation
In order to validate the obtained model, we employ a quantitative measure for predictive quality assessment, which is often used in the literature, namely, the prediction mean square error (PMSE). It is calculated according to Equation (30).
Furthermore, there are two more quantitative measures that are often used in the literature to compare the simulated and the measured signal, such as the BIS-index trajectory: the median performance error (MDPE), which is calculated according to Equation (31), and median absolute performance error (MDAPE), which is calculated according to Equation (32).
In Equations (30)-(32), x sim and x data denote the simulated and the measured data, respectively, whereas N stands for the number of data points in the dataset. Figure 6 shows the output of the residual dynamic model y res compared to the measurements-based trajectory y res,id, f . Finally, the output of the improved model is calculated according to Equation (7). The relevant BIS trajectories are shown in Figure 7. For the case of the particular patient and medical procedure, the calculated criteria for y res and BIS are given in Table 3 and Table 4, respectively. The presented residual dynamic model reduces the PMSE of the difference between y res and y res,id, f from 166.0 to 50.8; therefore, the improvement factor is 3.27. The implementation of the improved model reduces the PMSE from 179.8 (BIS BMS and BIS data ) to 63.3 (for BIS sim and BIS data ). In this case, the improvement factor is 2.84. Furthermore, the MDPE criterion improved by 10.24%, and MDAPE by 5.64%.
The main limitation of the proposed model is that it cannot be rigorously validated for a specific patient before its implementation because of its dependence on online measurements collected during the medical procedure. Moreover, no two anesthetic applications are alike, even if the same patient undergoes the same type of surgery several times, for example. Since the dynamic characteristics of the patient may change over time, it is not possible to compare the performance in two different interventions even if the same patient is treated. Because we cannot claim that the individualized dynamic model in question is time-invariant in the long run, it cannot be validated by using a special validation dataset that is strictly separate from the dataset used for identification. Furthermore, as the goal of the additional residual model is to consider the individual patient's discrepancy from the population-data-based model, it is not sensible to cross-validate it using another patient's measurements.
On the other hand, the results suggest that a significant improvement can be achieved, compared with the approach based only on the population-data-based model. In the future, the modeling framework will be further verified with several more datasets on different individuals treated with TIVA.

Conclusions
The article introduces a novel structure for modeling DoH, resulting in an improved individualized assessment of DoH, which is reflected in better predictions of the measured BIS signal. The presented model structure for modeling DoH dynamics employs the residual dynamic model add-on, which is used to model a particular treated patient's individual dynamic discrepancies from the population-data-based model. In such a manner, the model can take into account the patient's individual sensitivity to the anesthetic agent, which is not the case when using the population-data-based model exclusively, e.g., when using TCI, as is often the case in clinical practice. Therefore, the proposed modeling framework provides an improved mechanism for predicting DoH measured by the BIS index.
The improved model was verified using real clinical data logged during a medical treatment of a particular patient that lasted a little more than one hour. The results show that the predictions of the BIS-index trajectory were, indeed, considerably improved. Hence, the improved model seems to provide a solid foundation for better simulations as well as for the implementation in closed-loop model-based predictive control of DoH. The modeling framework will be further verified with several more datasets concerning various individuals that have been treated with TIVA.
To sum up, the presented framework provides a basis for a more patient-oriented individualized model for assessing DoH. The model seems to provide a deeper insight into DoH dynamics, which should lead to better administration methods that will relieve the anesthesiologist's workload and will benefit the patient by providing improved safety, individualized treatment, and, thus, alleviation of possible adverse effects during and after surgery.

Funding:
The authors acknowledge the financial support from the Slovenian Research Agency [research core funding number P2-0219].

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data used are available upon request to the corresponding author.