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.

METHODS
We have previously described a closed-loop controller applying a patient-individualized sigmoid E _{max} model as controller input (^{2} ). This controller estimated the patient-specific pharmacodynamic relationship (sigmoid E _{max} model) using the measured BIS from the induction phase. The modified controller developed and evaluated in this study uses a patient-specific sigmoid E _{max} 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
Figure 1 shows the schematic representation of the system.

Figure 1.:
Flow chart of the closed-loop system. The straight lines represent the closed-loop control system. At each time the required effect-site concentration is calculated by the controller. This value is sent to an additional algorithm considering the safety limits. The result of these calculations is the required effect-site concentration sent to the target-controlled infusion (TCI) algorithm which steers a pump injecting propofol to the patient. The measured Bispectral Index (BIS) is used as the input of the closed-loop controller. The dotted lines represent the Bayesian sigmoid E _{max} model estimator. The estimator receives a priori information from the population sigmoid E _{max} model, the optimal Bayesian variances for control and the patient measured BIS values.

Overall Controller Action
The closed-loop controller applies a sigmoid E _{max} model describing the relation between BIS and propofol effect-site concentration, and is characterized by four parameters: E _{0} equals the BIS value at no drug effect; E _{max} equals the change in BIS between no drug effect and maximum drug effect); c _{50} represents the propofol effect-site concentration at 50% of effect and γ represents the steepness of drug effect around 50%. The model can be presented as:

The controller always uses the latest available sigmoid E _{max} 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 (RUGLOOP^{1} ). 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 E_{MAX} Model Estimator Integrating Bayesian Approach
The model estimator will continuously calculate a sigmoid E _{max} model based on the patient's drug response. A standard approach to fit a model to a set of measured samples is using a least-squares technique by minimizing

over a set of measured points. In Eq. 2 , BIS_{sample} is the observed value and BIS_{estimated} 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 E _{max} 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 E _{max} 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 E _{max} model. This yields the following objective function

The second to fifth terms account for the covariates that determine the sigmoid E _{max} 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 E _{max} 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, sample_{TO} represents the sample time out, which can also be defined as the forgetting factor, t is the current time and t _{sample} 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
The BIS values, which we used as the controlled variable, can have any value between 100 and 0. Therefore, both E _{0} and E _{max} are fixed at 100 to simplify the identification of the patient's own sigmoid E _{max} model. This leaves only two parameters to estimate using the Bayesian technique: c _{50} and γ .

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 E _{max} 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 controller uses the reference sigmoid E _{max} model whenever the model estimator does not obtain a valid model. (e.g., at start of case when no patient-specific information is available). A good selection of the reference model should result in safe, effective control in the absence of patient-specific data samples. The reference model was calculated from previous work regarding propofol and BIS (^{19} ) as shown in Table 1 . As our closed-loop controller uses a fixed E _{0} and E _{max} (Eq. 4 ) both at 100, new c _{50} and γ typical values were modeled on the raw data using NONMEM V pharmacokinetic-dynamic modeling software (Globomax, Ellicott City, MD) whereby the values for E _{0} and E _{max} were fixed. This yields the modified reference model as applied in our controller.

Table 1: Original and Modified Reference Sigmoid E _{max} Model

The reference value for the total delay in the system, D _{population} , 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
Identification of the Bayesian Variances to Be Optimized
The aim of this study was to apply simulations to select the Bayesian variances, represented by the various VAR (^{2} ) terms in Eq. 5 , yielding the optimal controller for a Bayesian-based closed-loop system for propofol administration using the BIS as a controlled variable.

As indicated earlier, we used the Bayesian technique to tune two parameters of the sigmoid E _{max} model to be used in closed-loop: c _{50} and γ . The following terms from Eq. 5 are retained in the optimization:

Bayesian variance for squared deviation in c _{50} , called “VARc _{50} ”
Bayesian variance for squared deviation in γ , called “VAR_{γ} ”
Bayesian variance for the squared deviation in delay, called “VAR_{Delay} ”
The optimal number of samples to consider for model estimation, represented by “Sample_{TO} ”
Bayesian variance for the squared deviation in samples, called “VAR_{Sample} .“ The sum of the variances is only a scaling factor for the objective function, so we can arbitrarily choose the Bayesian variance for the samples equal to 1. We only need to optimize for the three others and Sample_{TO} .
Each set of Bayesian variances and Sample_{TO} 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 .

Table 2: Values for Bayesian Variances Used in Simulations

Simulation Population
Similar to our previous closed-loop performance study (^{7} ), we applied a simulation protocol to evaluate controller performance by adjusting parameters within the three components that comprised a study: 1) the virtual patient, 2) a stimulus profile, and 3) the BIS target level. Our population consisted of 416 virtual patients, generated using population characteristics from previous work not connected to the closed-loop controller algorithms (independent dataset) (^{20} ). This study provided the pharmacological characteristics for 13 patients. Of these 13 patients, we created the simulation population of 416 virtual patients by varying the c _{50} , γ , ke0 and clearance parameters considering the coefficient of variation numbers.

Case Simulation
As shown in Figure 2 , a BIS offset was composed to emulate a typical stimulus trajectory similar to a surgical case and was already found to be accurate in previous work (^{7} ). The total case time is exactly 3 h. The simulation trial for each virtual patient was run at three different control targets: target BIS values of 30, 50, and 70.

Figure 2.:
Bispectral Index (BIS) offset over time emulating a typical stimulus trajectory of surgical case. The total case time is exactly 3 h, including a virtual induction, maintenance and skin closure. Stimulus A simulates the arousal due to laryngoscopy/intubation; B represents surgical incision followed by a period of no surgical stimulation; C represents an abrupt stimulus after a period of low level stimulation; D shows on- and offset of normal surgical stimulation and; E simulates the withdrawal of stimulation during the closing period.

Processing of Simulation Results
A total of 625 different controllers were simulated, each against all patients. Each controller was tested at three different target values for BIS: 70, 50, and 30, representing the target levels for sedation, anesthesia and deep anesthesia, respectively.

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 N_{i} 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.

Figure 3.:
Schematic representation of the induction characteristics “first time to target ” and “overshoot.”

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 “VARc _{50} ,” “VAR_{γ} ,” “VAR_{Delay} ,” and “Sample_{TO} ,” the Q factors were clustered in bins for each value of a single “VAR” parameter, whereas Q is ranked within those bins. In the figures (see below), the leftmost peak in every bin represents the best controller within each bin.

Finally, the overall MDAPE for best and the worst controller and the final selected controller (called “overall best controller”) at each target was studied.

RESULTS
All 780,000 simulation cases (i.e., 416 patients × 625 controllers × 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 VAR_{Delay} and VARc _{50} differ between the best controllers for the three targets.

Table 3: Q Factor as Influenced by the Most Optimal and Worst Controller Behavior Resulting from the Bayesian variances' Selection

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 (Sample_{TO} , VARc _{50} , VAR_{γ} , and VAR_{Delay} ). Within a chart, the Q factor values are grouped for each single value of the Bayesian variance and sorted within this group from best to worst. The influence of this Bayesian variance on the Q factor shows best from the sequence of the peaks as a function of this variable. The best controller (highest Q factor) is marked with an arrow.

Figure 4.:
Pooled ranking of the controllers at Bispectral Index (BIS) target 30. The four charts each represent a ranking of all the Q factors for BIS target 30 versus one of the Bayesian variances (“Sample_{TO} ,” “VARc _{50} ,” “VAR_{γ} ,” and “VAR_{Delay} ”). Within a chart, the Q factor values are grouped for each single value of the Bayesian variance and sorted within this group from best to worst. The influence of this Bayesian variance on the Q factor shows best from the sequence of the peaks as a function of this variable. The best controller (highest Q factor) is marked with an arrow.

Figure 5.:
Pooled ranking of the controllers at Bispectral Index (BIS) target 50. The four charts each represent a ranking of all the Q factors for BIS target 50 versus one of the Bayesian variances (“Sample_{TO} ,” “VARc _{50} ,” “VAR_{γ} ,” and “VAR_{Delay} ”). Within a chart, the Q factor values are grouped for each single value of the Bayesian variance and sorted within this group from best to worst. The influence of this Bayesian variance on the Q factor shows best from the sequence of the peaks as a function of this variable. The best controller (highest Q factor) is marked with an arrow.

Figure 6.:
Pooled ranking of the controllers at Bispectral Index (BIS) target 70. The four charts each represent a ranking of all the Q factors for BIS target 70 versus one of the Bayesian variances (“Sample_{TO} ,” “VARc _{50} ,” “VAR_{γ} ,” and “VAR_{Delay} ”). Within a chart, the Q factor values are grouped for each single value of the Bayesian variance and sorted within this group from best to worst. The influence of this Bayesian variance on the Q factor shows best from the sequence of the peaks as a function of this variable. The best controller (highest Q factor) is marked with an arrow.

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 VARc _{50} _{.} Minimal dependency is observed for SAMPLE_{TO} . 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 .

Table 4: Bayesian Variances Yielding the Best Results in the Current Trials

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.

Table 5: Overall Median Absolute Prediction Error (MDAPE) for the Best, the Worst, and the Final Selected Controller (Called “Overall Best Controller”) at Each Target

DISCUSSION
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 E _{max} 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 E _{max} 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 VARc _{50} , 0 for VAR_{γ} , 30 for VAR_{Delay} , and 120 for Sample_{TO} . Compared to the proposed universe of Bayesian variances sets, as shown in Table 2 , this set, providing optimal control, allows relatively little intersubject variability. This means that the improvement in control performance resulting from the larger degree of freedom did not compensate for the drawback introduced by the cost function.

Verification of the dependency of the Q factor on VARc _{50} in Figures 5 and 6 for targets 50 and 70 shows that the controller performance decreases when VARc _{50} is zero. An underlying explanation is that the modeling algorithm cannot properly model incoming data centered around a different c _{50} by varying γ only. This leads to absurd models and, accordingly, absurd control behavior results from it. It can be concluded that the model does best when it can differentiate between sensitive patients and resistant patients.

As for the VAR_{Delay} , 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 Sample_{TO} 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 c _{50} could mean a possible threat to people who cannot tolerate high propofol concentrations due to hemodynamic instability. One might make the initial population models' c _{50} 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.

REFERENCES
1. Linkens DA, Hacisalihzade SS. Computer control systems and pharmacological drug administration: a survey. J Med Engl Technol 1990;14:41–54

2. Struys MM, De Smet T, Versichelen LF, Van De Velde S, Van den Broecke R, Mortier EP. Comparison of closed-loop controlled administration of propofol using Bispectral Index as the controlled variable versus “standard practice ” controlled administration. Anesthesiology 2001;95:6–17

3. Schwilden H, Stoeckel H, Schuttler J. Closed-loop feedback control of propofol anaesthesia by quantitative EEG analysis in humans. Br J Anaesth 1989;62:290–6

4. Morley A, Derrick J, Mainland P, Lee BB, Short TG. Closed loop control of anaesthesia: an assessment of the bispectral index as the target of control. Anaesthesia 2000;55:953–9

5. Absalom AR, Sutcliffe N, Kenny GN. Closed-loop control of anesthesia using Bispectral index: performance assessment in patients undergoing major orthopedic surgery under combined general and regional anesthesia. Anesthesiology 2002;96:67–73

6. Mortier E, Struys M, De Smet T, Versichelen L, Rolly G. Closed-loop controlled administration of propofol using bispectral analysis. Anaesthesia 1998;53:749–54

7. Struys MM, De Smet T, Greenwald S, Absalom AR, Binge S, Mortier EP. Performance evaluation of two published closed-loop control systems using bispectral index monitoring: a simulation study. Anesthesiology 2004;100:640–7

8. Locher S, Stadler KS, Boehlen T, Bouillon T, Leibundgut D, Schumacher PM, Wymann R, Zbinden AM. A new closed-loop control system for isoflurane using bispectral index outperforms manual control. Anesthesiology 2004;101:591–602

9. Kenny GN, Mantzaridis H. Closed-loop control of propofol anaesthesia. Br J Anaesth 1999;83:223–8

10. Mortier E, Struys M, De Smet T, Versichelen L, Rolly G. Closed-loop controlled administration of propofol using bispectral analysis. Anaesthesia 1998;53:749–54

11. Sheiner LB, Beal S, Rosenberg B, Marathe VV. Forecasting individual pharmacokinetics. Clin Pharmacol Ther 1979;26:294–305

12. Garraffo R, Iliadis A, Cano JP, Dellamonica P, Lapalus P. Application of Bayesian estimation for the prediction of an appropriate dosage regimen of amikacin. J Pharm Sci 1989;78: 753–7

13. Schnider TW, Minto CF, Bruckert H, Mandema JW. Population pharmacodynamic modeling and covariate detection for central neural blockade. Anesthesiology 1996;85:502–12

14. Schnider TW, Minto CF, Gambus PL, Andresen C, Goodale DB, Shafer SL, Youngs EJ. The influence of method of administration and covariates on the pharmacokinetics of propofol in adult volunteers. Anesthesiology 1998;88:1170–82

15. Minto CF, Schnider TW, Gregg KM, Henthorn TK, Shafer SL. Using the time of maximum effect site concentration to combine pharmacokinetics and pharmacodynamics. Anesthesiology 2003;99: 324–33

16. Schnider TW, Minto CF, Shafer SL, Gambus PL, Andresen C, Goodale DB, Youngs EJ. The influence of age on propofol pharmacodynamics. Anesthesiology 1999;90:1502–16

17. Struys MM, De Smet T, Depoorter B, Versichelen LF, Mortier EP, Dumortier FJ, Shafer SL, Rolly G. Comparison of plasma compartment versus two methods for effect compartment–controlled target-controlled infusion for propofol. Anesthesiology 2000;92: 399–406

18. Press W, Teukolsky S, Vetterling W, Flannery B: Numerical recipes in C: the art of scientific computing. 2nd ed. Cambridge, UK: Cambridge University Press, 1992

19. Vanluchene AL, Vereecke H, Thas O, Mortier EP, Shafer SL, Struys MM. Spectral entropy as an electroencephalographic measure of anesthetic drug effect: a comparison with bispectral index and processed midlatency auditory evoked response. Anesthesiology 2004;101:34–42

20. Vereecke HE, Vasquez PM, Jensen EW, Thas O, Vandenbroecke R, Mortier EP, Struys MM. New composite index based on midlatency auditory evoked potential and electroencephalographic parameters to optimize correlation with propofol effect site concentration: comparison with bispectral index and solitary used fast extracting auditory evoked potential index. Anesthesiology 2005;103:500–7

21. Varvel JR, Donoho DL, Shafer SL. Measuring the predictive performance of computer-controlled infusion pumps. J Pharmacokinet Biopharm 1992;20:63–94

22. Kansanaho M, Hynynen M, Olkkola KT. Model-driven closed-loop feedback infusion of atracurium and vecuronium during hypothermic cardiopulmonary bypass. J Cardiothorac Vasc Anesth 1997;11:58–61

23. Absalom AR, Kenny GN. Closed-loop control of propofol anaesthesia using bispectral index: performance assessment in patients receiving computer-controlled propofol and manually controlled remifentanil infusions for minor surgery. Br J Anaesth 2003;90:737–41

24. O'Hara DA, Bogen DK, Noordergraaf A. The use of computers for controlling the delivery of anesthesia. Anesthesiology 1992;77:563–81

25. Liu N, Chazot T, Genty A, Landais A, Restoux A, McGee K, Laloe PA, Trillat B, Barvais L, Fischler M. Titration of propofol for anesthetic induction and maintenance guided by the bispectral index: closed-loop versus manual control: a prospective, randomized, multicenter study. Anesthesiology 2006;104: 686–95

^{1} RUGLOOP, 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 . Cited Here