| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
BACKGROUND: Implementing Bayesian methods in a model-based closed-loop system requires the integration of a standard response model with a patient-specific response model. This process makes use of specific modeling weights, called Bayesian variances, which determine how the specific model can deviate from the standard model. In this study we applied simulations to select the Bayesian variances yielding the optimal controller for a Bayesian-based closed-loop system for propofol administration using the Bispectral Index (BIS) as a controlled variable. METHODS: The relevant Bayesian variances determining the modeling process were identified. Each set of such Bayesian variances represents a potential controller. The set, which will result in optimal control, was estimated using calculations on a simulated population. We selected 625 candidate sets. Similar to our previous closed-loop performance study, we applied a simulation protocol to evaluate controller performance. Our population consisted of 416 virtual patients, generated using population characteristics from previous work. A BIS offset trajectory similar to a surgical case was used. RESULTS: We were able to develop, describe, and optimize the parameter setting for a patient-individualized model-based closed-loop controller using Bayesian optimization. Selection of the optimal set yields a controller performing with the following median absolute prediction errors at BIS targets 30, 50, and 70: 12.9 ± 2.87, 7.59 ± 0.74, and 5.76 ± 1.03 respectively. CONCLUSIONS: We believe this system can be introduced safely into clinical testing for both induction and maintenance of anesthesia under direct observation of an anesthesiologist.
Automated control of anesthesia using continuous feedback of a patient's clinical state might improve the quality of anesthetic drug administration, resulting in a better control of anesthesia and faster recovery (1,2). Various closed-loop systems have been proposed to guide IV drug administration using spontaneous electroencephalographic (2–8) or auditory evoked-potential indices (9) as controlled variables. However, limitations of these control systems have precluded their adoption into clinical practice. Most controllers are designed to maintain the anesthetic only during stable periods of anesthesia, and thus exclude induction. Model-based controllers reported in the literature may be used during induction; however, these controllers typically assume an initial drug-free patient, and thus cannot be used in already anesthetized patients (2,6,7). Finally, the concentration versus response profile is likely to change over the course of anesthesia, mostly as a result of changing levels of noxious stimulation, but also because of possible pharmacodynamics (i.e., acute tolerance to opioids) or because of changes to the patient (e.g., liver transplantation). In an attempt to resolve these shortcomings, the Bayesian approach was included in our existing control algorithm (2,7,10). Bayesian optimization, as proposed by Sheiner et al. (11) is classically used in pharmacokinetic-dynamic modeling to individualize a drug dosing regimen by combining individual information with the knowledge of an a priori probability density function containing the statistical properties of the parameter to be estimated (12). The Bayesian method starts from a standard, population-based response model providing the previous distribution of parameter values. These values are adjusted to reflect the patient's own parameters over time, based on the observed response of the individual patient under varying circumstances. This process makes use of specific modeling weights, called Bayesian variances, which determine how the patient-specific model can deviate from the population model. These Bayesian variances need to be optimized for control performance in a target population. The aim of this study was to apply simulations to select the Bayesian variances yielding the optimal controller for a Bayesian-based closed-loop system for propofol administration using the Bispectral Index (BIS) as a controlled variable.
We have previously described a closed-loop controller applying a patient-individualized sigmoid Emax model as controller input (2). This controller estimated the patient-specific pharmacodynamic relationship (sigmoid Emax model) using the measured BIS from the induction phase. The modified controller developed and evaluated in this study uses a patient-specific sigmoid Emax model continuously adjusted to the observed responses over the entire course of the anesthetic using a Bayesian technique. The steps in this process are described below.
Closed-Loop System Structure
Overall Controller Action
The controller always uses the latest available sigmoid Emax model calculated by the Bayesian model estimator, as described below. The output of the controller is a predicted propofol effect-site concentration, which will be used as the input parameter of an effect-compartment controlled infusion system. This computer-assisted continuous infusion device targets an effect-site concentration using a three-compartment model completed with an effect-site compartment (RUGLOOP1). For propofol, the pharmacokinetic-dynamic model published by Schnider et al. was used (13,14). The propofol effect-site concentration was computed to yield a time to peak effect of 1.6 min after bolus injection, as published by Schnider et al. (15,16) and clinically confirmed by Struys et al. (17)
Sigmoid EMAX Model Estimator Integrating Bayesian Approach
over a set of measured points. In Eq. 2, BISsample is the observed value and BISestimated is the estimated value based on the model to be fitted. The sum is called the objective function. For applying the Bayesian technique in the estimator, a reference sigmoid Emax model is used as a priori information, whereas the patient-specific propofol effect site-concentration versus measured BIS data pairs are used to tune this reference or population model towards the sigmoid Emax model in the individual patient. To incorporate the a priori knowledge while using the least-squares technique, we introduced additional terms in the objective function, representing a penalty for the difference between the reference model (the a priori knowledge) and the resulting sigmoid Emax model. This yields the following objective function
The second to fifth terms account for the covariates that determine the sigmoid Emax model to be estimated. "Population" is the original population reference model parameter (= a priori information), and "Estimated" is the estimated value of parameter for the individual. The minimization of Eq. 3 without any measured samples will yield the reference model only. As observations are made over the course of a procedure, the model in the individual will gradually diverge from the original population model. With many samples, the optimal solution becomes primarily one that minimizes the difference between the measured and predicted observations, and the reliance on the original population model diminishes. This is the nature of Bayesian inference. Both the observations and the model parameters may vary in magnitude over several orders of magnitude. When this is the case, the model primarily works to minimize the error in the large observations, and minimizes the deviation in the large parameters from the initial population estimate, as the larger numbers typically have much larger squared errors. This behavior is undesirable, and is most appropriately handled by weighting the squared error by the expected variance. In the case of the model parameters, the expected variance represents the Bayesian uncertainty about the parameter estimate. As a consequence, the variance used to weight the squared difference between parameter estimate and population estimate will be defined as the Bayesian variances. The patient-specific sigmoid Emax model's characteristics as calculated continuously by the estimator will depend on the population model, the patient's drug response as well as the applied Bayesian variances (Fig. 1). The optimal values for the Bayesian variances depend on the specific terms to which they are applied. For the observations, they typically reflect the variance in the measurements themselves. For the parameters, they reflect the interindividual variability in the reference model parameters. Or they could be selected empirically using simulations to obtain the best controller behavior according to a predefined criterion in a closed-loop application. This is the intent of the present investigation. As the patient response during induction might not be representative for the response during surgery, and because the concentration versus response relationship may change over time, it may be undesirable to calculate the current patient-specific drug-response model using the full history of patient-specific samples. Therefore, the samples considered for the modeling are time-limited using a forgetting factor: as samples get older, their contribution in the modeling is reduced, until they vanish completely. The resulting objective function incorporating these effects is shown in Eq. 4.
In this equation, the VARs are the Bayesian variances for each covariate, sampleTO represents the sample time out, which can also be defined as the forgetting factor, t is the current time and tsample is the time when each sample was captured. MAX[] means that the contribution of a specific sample in the objective function is only appreciated while it is larger than zero. The actual minimization of the objective function uses the Levenberg–Marquardt method for fitting using nonlinear models (18).
Applied Estimator Parameters Shortlist As learned from pilot trials (data not shown), the overall delay or timing mismatch of the controlled system was important in our estimator. This total delay of the system is an overall estimate that includes misspecification of plasma-effect site equilibration, IV line flow delay, and monitoring delay. The initial trials showed a vast improvement in the modeling (not yet control) when an attempt was made to estimate the overall delay as well.
Whereby D is the delay to be estimated. The patient-specific sigmoid Emax model resulting from the objective function shown in Eq. 5 is used in the control algorithm as shown in Figure 1.
Population Reference Model Selection
The reference value for the total delay in the system, Dpopulation, was set at 10 s, which represents the cumulative delay generated from the delay from the A 2000 BIS-XP monitor (Aspect Medical, Newton, IL), infusion delay, plus a possible mismatch in the pharmacokinetic-dynamic link model (represented by a population ke0).
Bayesian Variances Optimization for Closed-Loop
As indicated earlier, we used the Bayesian technique to tune two parameters of the sigmoid Emax model to be used in closed-loop: c50 and
Each set of Bayesian variances and SampleTO as defined above represents a potential controller. The set that will result in optimal control for a Bayesian-based, closed-loop system for propofol administration using the BIS as the controlled variable, was estimated using calculations on a simulated population (defined below). To limit the number of evaluations to be done, we selected 625 candidate sets, assuming continuous performance in-between. The candidate sets were chosen after preliminary simulations (data not shown), which indicated good starting points, while still broad enough to cover all useful values. The evaluated sets are shown in Table 2.
Simulation Population
Case Simulation
Processing of Simulation Results Based on the method of Varvel et al. (21), previously applied by Kansanaho et al. for the performance of a closed-loop system for muscles relaxants (22) and hypnotics (2,5,23), the overall performance of all controllers was characterized on the basis of the following parameters. First, using all observations within the period, the performance error (PE) was calculated according to the formula
Subsequently, the inaccuracy (median absolute performance error [MDAPE]) was calculated as follows:
where Ni is the number of values PE obtained for the i-th subject. Standard errors (SE) were calculated using the two-stage approach as described by Varvel et al (21). Additionally, as shown in Figure 3, the induction phase parameters average first time to target (AVGFTTT) and worst time to target (MAXFTTT), required to reach lowest BIS (from start of induction), average overshoot (AVGOVS), and worst overshoot (MINOVS) respectively, were used to compare both the speed and stability of induction between the controllers. The overshoot was used as a negative value.
With larger values for the Bayesian variances, the freedom of the model to be calculated for the patient is increased, but the chances to get an "off " model is increased as well. Therefore, we introduced a cost function (COST), which penalized a higher degree of freedom (larger Bayesian variance values). To ease the search for the best controller, we grouped all criteria in an overall quality number. This one metric combined all the performance parameters in a balanced way by normalizing the parameters for each controller versus the minimum (Min) and maximum (Max) of those in the group of controllers, and by multiplying with a power factor afterwards. The weight of each term in the sum is chosen to obtain an equal importance of the overall MDAPE versus all other parameters (induction characteristics and overall cost).
This formula results in a Q value between 0 and –1000, with 0 the best performing controller. Still, to avoid that controllers performing out of range for a criterion would show up in an acceptable range, they were discriminated with a large value (–1000) for the corresponding criterion, widening the scale from 0 to –6000.
The controller with the best Q factor was found using Excel (Microsoft, Redmond, WA) by listing all controllers and their Q factors. As it might be undesirable to use a different controller at every different target, an overall best controller will be suggested. To asses a systematical evolution of the Q factor versus the separate variables "VARc50," "VAR Finally, the overall MDAPE for best and the worst controller and the final selected controller (called "overall best controller") at each target was studied.
All 780,000 simulation cases (i.e., 416 patients x 625 controllers x 3 target BIS values) were included in the analysis. We investigated the controller performance during both the induction phase and the maintenance phase. Table 3 shows the Q factors and the Bayesian variance values for the best and worst behaving controllers for the three targets. It can be observed that only VARDelay and VARc50 differ between the best controllers for the three targets.
To further explore these differences, Figures 4, 5, and 6 plot a pooled ranking of the controllers at each target. The four charts within a figure each represent a ranking of all the Q factors for the BIS target versus one of the Bayesian variances (SampleTO, VARc50, VAR
As a secondary measurement, one could observe the sequence of the rightmost peaks indicating the worst controller for each group, but we are little interested in the evolution of the worst controller. For a target BIS of 30, the evolution of the best controllers is most marked as a function of VAR , even though the other variables still show a clear trend. For a target BIS of 50, the most express dependency is in both VAR and VARc50. Minimal dependency is observed for SAMPLETO. Lastly, a target 70 displays a clear downgrade in performance of the best controller with increasing value for all variables. To select the overall best performing controller, we considered the three BIS targets simultaneously, and weighed the performance of each BIS target's best controller at the other BIS targets. It showed that the best controller for BIS target 50 was ranked 2nd best performing at BIS target 70 and 3rd best performing at BIS target 30. However, the best controller at BIS target 30 was ranked 7th and 23rd at BIS targets 50 and 70 respectively. Finally, the optimal controller for BIS target 70 performed as 32nd and 31st best at the BIS targets 30 and 50. As such, the overall preferred controller clearly was the one which performed best at BIS target 50. The resulting controller Bayesian Variances are shown in Table 4.
Finally, the behavior of the best, the worst, and the final selected controller (called "overall best controller") at each target is detailed. Table 5 shows MDAPE for each specific controller.
We applied simulations to estimate the optimal Bayesian variances for a Bayesian-based closed-loop system for propofol administration using the BIS as a controlled variable. Our final goal was to develop a practical, patient-individualized model-based closed-loop controller using Bayesian optimization that can be introduced safely into clinical testing. The introduction of Bayesian methods in our previously published (2) model-based adaptive control strategy will allow the use of the controller in a more universal way. The use of a priori information will permit the use of the controller during anesthesia induction, a distinct advantage compared to previously reported control algorithms (24). It can be used on an already anesthetized patient, where the controller will keep the total anesthetic state at a steady level, whereas the patient's drug history (even unknown) is gradually washing out. The mismatch between the controller-perceived required drug concentration and the total patient drug concentration is actually irrelevant, and compensated for by the Bayesian modeling. Similar to other control algorithms (24), a Bayesian controller's performance is influenced by multiple parameters, called Bayesian variances, which are constant numbers requiring fine tuning to guarantee optimal behavior. This fine-tuning implies testing under extreme clinical circumstances to establish fully the safety, efficacy, and reliability of closed-loop anesthesia before introduction into clinical testing. As it might be considered morally irresponsible to stress human subjects under target effect settings and surgical stimuli beyond those accepted as good clinical practice, we felt that computer simulation of patients and intraoperative events would enable a more thorough characterization of controller responses to variation inpatient types and interventions. The applied simulation methods were previously used to compare different closed-loop control algorithms (7). As we use the sigmoid Emax model obtained using Bayesian approach for our closed-loop algorithm, we need to optimize the Bayesian variances towards this application. The optimal selection should result in a Bayesian sigmoid Emax model estimator providing the best model for accurate closed-loop control for all BIS targets. After optimization to three clinically significant BIS targets (30,50,70), we selected a single overall optimal controller for future testing. Even though part of the selection criteria was based on overshoot, undershoot and avoiding oscillations, we do not intend to imply that this study tested robustness of a particular controller, which must be done in further research.
The optimal values for the Bayesian variances resulting from our simulations were 5 for VARc50, 0 for VAR
Verification of the dependency of the Q factor on VARc50 in Figures 5 and 6 for targets 50 and 70 shows that the controller performance decreases when VARc50 is zero. An underlying explanation is that the modeling algorithm cannot properly model incoming data centered around a different c50 by varying As for the VARDelay, one might wonder why a delay to be estimated was introduced. As the total delay is an overall estimate that groups pharmacokinetic-dynamic modeling mismatch or "ke0 mismatch," IV line flow delay, and artifact-dependent monitoring delay, several uncertainties contribute to an overall mismatch between the expected time behavior versus the actual time behavior. One could try to separately model the ke0 and the delay, but this only leads to more difficult identification of the various model parameters from the measured data in the modeling procedure. Moreover, using a varying ke0 would lead to difficult interpretation of effect-site concentration by the anesthesiologist in future clinical applications. The value of 120 for SampleTO means that only the last 2 min of patient data are considered to build the model. One can wonder whether the patient history of 15 min ago reflects the current patient state better than the population model. The simulation results seem to reveal that this is not the case. Most probably, the population model has the advantage of being a considered starting value providing good basic control. Selection of accurate a priori information is the foundation of good Bayesian control. All Bayesian methods start from relevant measured data (11). The population model, which was used by the controller when no patient-specific measurements were acquired, determines the control behavior during induction. We selected our a priori information based on previous work as shown in Table 1. Although this selection resulted in satisfying results in the current simulations, a (relatively) high c50 could mean a possible threat to people who cannot tolerate high propofol concentrations due to hemodynamic instability. One might make the initial population models' c50 dependent on patient covariates including age, ASA classification, and concurrent disease. The possible benefit of this might be the subject of further research. To introduce the final selected controller into clinical testing, its safety performance needs to be stressed. As shown in Table 5, MDAPE of the final selected controller are close to the best at each target and are better than the worst controller at each target. Even more, the MDAPE results from this simulation study show better results than the MDAPEs from our own previous work using similar simulations but a model-based controller without the Bayesian approach (7). The MDAPEs are also better or comparable with various clinically applied closed-loop systems (2,8,25) but this Bayesian-based controller offers much more freedom during clinical applications as stated above. In conclusion, we were able to develop, describe, and optimize the parameter setting for a patient-individualized model-based closed-loop controller using Bayesian optimization. We believe this system can be introduced safely into clinical testing under direct observation of an anesthesiologist.
1RUGLOOP, written by Tom De Smet, MSc (Electronic Engineering) and Michel MRF Struys, MD, PhD (Professor in Anesthesia). More information at http://www.anesthesia-uzgent.be. Accepted for publication July 24, 2007. Dr. Steven Shafer, Editor-in-Chief, was recused from all editorial decisions related to this manuscript. Disclosure: Michel M.R.F Struys, MD, PhD; Tom De Smet, MSc; and Steven L Shafer, MD are patent holders of the "model-based" closed-loop system. Scott Greenwald, PhD is an employee of Aspect Medical, the developers of BIS®. The new controller framework has been published in detail in the patent application WO 2005/072792 A1 (see http://www.wipo.int/pctdb/en/search-adv.jsp). The authors Tom De Smet and Michel M.R.F Struys have equally contributed to this manuscript.
This article has been cited by other articles:
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|