Real-time Closed-loop Control in a Rodent Model of Medically Induced Coma Using Burst Suppression

MEDICALLY induced coma is an anesthetic state of profound unconsciousness and brain inactivation created to treat status epilepticus and to facilitate recovery after traumatic brain injuries.^1–3  When treating status epilepticus, a hypnotic, such as propofol or a barbiturate, is used to directly inhibit seizure activity.^2,3  After a brain injury these drugs are administered to provide brain protection by reducing cerebral blood flow and metabolism.¹  In both cases the anesthetic is titrated to achieve a specific clinical target that indicates a state of large-scale brain inactivation. A standard approach is to monitor the patient’s brain activity continuously with an electroencephalogram and use a specified level of burst suppression as the target. Burst suppression is an electroencephalogram pattern indicating a state of highly reduced electrical and metabolic activity during which periods of electrical bursts alternate with isoelectric periods termed suppressions.^4–6 

No established guidelines exist for specifying the level of burst suppression required for a medically induced coma. A target level is chosen, and control of that level is managed by continually monitoring the electroencephalogram and manually adjusting the drug infusion rate. A common goal of medically induced coma is maintaining a reduction in brain activity for 24 h or more, periods significantly longer than any human operator can maintain tight control. Defining a precise, quantitative target level of burst suppression and designing a closed-loop anesthetic delivery (CLAD) system for maintaining that target would be a more efficient approach.

Closed-loop anesthetic delivery systems for control of unconsciousness and sedation have been extensively studied.^7–26  Although no CLAD system has been designed to manage medical coma in humans, Vijn and Sneyd²⁷  implemented a CLAD system to test new anesthetics in rodents using as the control signal the burst-suppression ratio; the fraction of time per 15 s that the electroencephalogram is suppressed. For several anesthetics they established nonmodel-based control of burst-suppression ratio levels measured in terms of group averages rather than individual control trajectories. Cotten et al.²⁸  studied methoxycarbonyl etomidate with this paradigm in rodents and also reported only group average control results.

We hypothesize that a CLAD system could precisely control burst suppression as a way to efficiently maintain a medically induced coma. We test this hypothesis by constructing a CLAD system to control burst suppression in real time in a rodent model, using electroencephalogram recordings and a computer-controlled infusion of propofol. The CLAD system uses a two-compartment pharmacokinetics model to characterize the effect of propofol on the electroencephalogram. We introduce as the control signal the burst-suppression probability (BSP), the instantaneous probability of the brain being suppressed computed from the electroencephalogram in real time. We estimate the pharmacokinetic model parameters online for individual rodents and use them to construct proportional-integral (PI) controllers. To evaluate performance of our CLAD system we establish new statistical criteria to assess reliability and accuracy at individual target levels of burst suppression, and a new Bayesian statistical approach to assess overall reliability and accuracy of the control experiments. We use our CLAD system to maintain precise control of burst suppression in individual rats.

These animal studies were approved by the Subcommittee on Research Animal Care, Massachusetts General Hospital, Boston, Massachusetts. Six male Sprague–Dawley rats (Charles River Laboratories, Wilmington, MA) weighing 377–460 g were used for these studies. Animals were kept on a standard day–night cycle (lights on at 7:00 am, and off at 7:00 pm), and all experiments were performed during the day. We use rats as our experimental system because they are an established model for study of burst suppression.^27,28 

Extradural electroencephalogram electrodes were surgically implanted at least 7 days before experimentation as previously described.^29,30  Briefly, general anesthesia was induced and maintained with isoflurane. A microdrill (Patterson Dental Supply Inc., Wilmington, MA) was used to make four holes at the following stereotactic coordinates: A0L0, A6L3, A6L-3, and A10L2 relative to the lambda.^29,30  An electrode with mounting screw and socket (Plastics One, Roanoke, VA) was screwed into each hole, and the sockets were inserted in a pedestal (Plastics One). The screws, sockets, and pedestal were all permanently fixed with dental acrylic cement, and the animal underwent a minimum recovery period of 7 days. The potential difference between electrodes A0L0 and A6L-3 (left somatosensory cortex) was recorded. The signal was referenced to A10L2 and recorded using a QP511 Quad AC Amplifier System (Grass Instruments, West Warwick, RI) and a USB-6009 14-bit data-acquisition board (National Instruments, Austin, TX). The sampling rate was 512 Hz. A line filter with cutoff frequencies of 0.3 and 50 Hz was used.

Rats were anesthetized in an induction chamber with 2–3% isoflurane in oxygen before the placement of a 24-gauge intravenous catheter in the tail. Isoflurane was discontinued and propofol anesthesia was initiated only after the animal regained the righting reflex and returned to a normal level of activity. A Physio 22 syringe pump (Harvard Apparatus, Holliston, MA) was used to deliver propofol. After the rat had loss of righting under propofol anesthesia, a rectal temperature probe was inserted and the animal was placed in the supine position on a heating pad to maintain core temperature between 36.5° and 37.4°C. Oxygen was provided by face mask to prevent hypoxemia.

Before describing the experimental protocol, we give a conceptual overview of the components of our CLAD system and how the system works. We constructed a CLAD system whose primary components were (fig. 1A): an electroencephalogram recording system; a computer-controlled infusion pump; a real-time segmentation algorithm to convert the continuous electroencephalogram into a binary time series; a real-time signal-processing algorithm to estimate the state of burst suppression from the binary time-series; and a PI control algorithm that issues commands to the infusion pump based on the burst-suppression state estimate.

To understand how our CLAD system works, we explain it schematically (fig. 1A). We introduce the concept of the BSP; a number between 0 and 1, which describes the instantaneous probability of the electroencephalogram being in a state of suppression.³¹  A BSP value of 0 corresponds to an active electroencephalogram with no suppression, whereas a value of 1 corresponds to a completely isoelectric or suppressed electroencephalogram. We assume that a target level of burst suppression, BSP_target, has been set as a number between 0 and 1 (step 0). We further assume that propofol is being administered by an infusion pump and that this infusion is producing a state of burst suppression, which we wish to control at BSP_target (steps 1 and 2). The time course of the BSP computed from the electroencephalogram is the quantity that our CLAD system tracks. To do so, we segment the electroencephalogram into a binary time series (step 3) with a 20 ms resolution in which a burst is a 0 and a suppression is a 1 (Eqs. 19–21). The binary time series is input to the BSP filter algorithm (Eqs. 7–13), which computes a real-time BSP estimate (step 4).

The estimated BSP is the feedback signal that our controller compares with BSP_target (step 5). The difference between the estimated BSP and BSP_target is the error signal. The error signal, transformed to concentration (Eqs. 14 and 15), is passed to a PI controller (step 6). The objective of the controller is to keep the error as close to 0 as possible, which means that the CLAD is maintaining the target BSP level. Therefore, the PI controller issues commands to the infusion pump to change the infusion rate based on the magnitude and sign of the error signal (step 7). The entire cycle from steps 1–7 takes 1 s, the update interval of our CLAD system.

A central component of our CLAD system is a second-order pharmacokinetics state model, which characterizes how the infusion of propofol changes the concentration in the brain compartment or effect site (Eqs. 1–2). We link this model to the BSP through a binomial probability model (Eqs. 3–4). By equation 4, the probability of being suppressed increases monotonically with the effect-site concentration. To use the state model in the PI controller we must estimate its parameters for each rat so that the model’s response properties are tailored to each animal. We estimated the model parameters for each animal by conducting a preliminary experiment before starting the control experiment. We accomplished this process, termed system identification, by administering one or more bolus doses of propofol, estimating the time course of the effect-site concentrations from the time course of the BSP and then fitting the pharmacokinetic state model to the estimated effect-site concentration data by nonlinear least squares (fig. 1B).

After the system identification, which required 10–15 min, we initialized the controller by requiring it to track a BSP_target of 0.2 for 5–15 min before starting the control experiment (fig. 1B). This allowed us to ensure correct communication between the software and the infusion pump. The initialization also ensured that each control experiment started with each animal in the same, low-level state of burst suppression. After the initialization, we set BSP_target to the target BSP trajectory selected for that animal and let the CLAD system continue BSP control.

For the control experiment, we defined six BSP target trajectories by permuting the order of three BSP levels, 0.4, 0.65, and 0.9, to test the ability of our CLAD system to achieve and maintain BSP control for each individual animal. We randomly assigned each of six rats to one of the six control trajectories. For each permutation, the duration of each level was 15 min. We set an approximate 5–10 min linear ramp between each level to give a total target control duration of approximately 60 min. We set the controller update interval at 1 s.

The electroencephalogram acquisition and segmentation, the BSP estimation, and the PI control were carried out in real-time using custom MATLAB (Version R13; Natick, MA) software run on an HP Probook 5430s (Hewlett Packard, Palo Alto, CA) laptop computer (fig. 1A). The controller software issued commands through an RS-232 serial connection to control the syringe pump (Harvard-22; Harvard Apparatus) infusion rate. We give the mathematical details of each component of the CLAD design and implementation in the appendix.

We measured performance of the CLAD system in terms of the error defined by

For each target level and for each transition we used (Eq. a) to obtain percent error defined as

where the SD is taken across all sampled data points. For transitions, we took BSP_target as the mean BSP traversed during the transition.

In addition to (Eq. b) we determined, for each level and transition, the median absolute deviation (MAD)

along with two common performance metrics used in the evaluation of clinical pharmacokinetics models.³² 

They are the median performance error (MDPE):

and the median absolute performance error (MDAPE):

where, in equations c–e, the median is taken across all sampled data points. The MDPE provides a measure of bias (in our case, steady state tracking offset), whereas the MDAPE is an alternative to (Eq. b), which is less sensitive to outlying errors. The MDAPE is simply the ratio of the MAD to the target level. To summarize performance we computed the median of each metric across animals at each level and the transitions. All statistics are computed in MATLAB.

We used the performance measures to define specific statistical criteria to assess the reliability and accuracy of our CLAD system.³³  A preliminary study of our system indicated that the absolute errors were less than 0.2. Therefore, for a given target level, we defined reliability of the CLAD system as the absolute error being less than 0.15 with high probability. We set that probability at 0.95. This criterion can be easily evaluated as it is equivalent to the 95th percentile of the absolute error distribution being 0.15 or less. We computed the absolute error distribution at a given level as the absolute values of e-tracking (Eq. a) at that level. We used 900 data points (60 points per min × 15 min per level) to compute the absolute error distribution at a level. There were six animals and three levels per animal or 18 levels in total. We assessed reliability on each of the 18 levels separately and overall by considering all levels across animals.

We defined accuracy of the CLAD system for a given level by the error distribution at that level being indistinguishable from zero. We computed the error distribution at a given level as the values of e-tracking (Eq. a) at that level. We also used 900 data points to compute the error distribution at a level. We considered the CLAD system control to be accurate at a given level if zero was inside the 2.5th and 97.5th percentile of the error distribution (95% CI). If zero was outside this 95% CI for a given level then we rejected the hypothesis that the system is accurate with a P value less than 0.05. If zero was below (above) the 2.5th (97.5th) percentile the systems had a positive (negative) bias. If zero was inside the 50% CI, i.e., between the 25th and 75th percentiles, we considered the system to be highly accurate. We assessed accuracy on each of the 18 levels separately and overall by considering all levels across animals.

We used Bayesian analysis to assess overall reliability of the CLAD system by combining data across levels.^29,34  We performed the overall reliability analysis across the 18 levels assuming levels within animals were independent. Independence is a reasonable assumption because if we assume a high first-order serial correlation of 0.98 between adjacent data points separated by 1 s and if we allow between-level transitions of 5–10 min then, the maximum correlation between the closest two points in immediately adjacent levels is between ((0.98)⁶⁰⁰ = 5.4 × 10⁻⁶; (0.98)³⁰⁰ = 2.3 × 10⁻³), where 300 (600) = 5 (10) min × 60 data points per min.³⁵  That is, control activity separated by 5 min or more is unrelated.

Let p denote the probability that the CLAD system is reliable at a level. The analysis of reliability across levels yields a binomial probability model with n = 18 assessments of which on k levels the system was reliable and on n − k levels the system was not reliable. Being reliable at a given level is defined as the system satisfying the reliability criterion that the 95th percentile of the absolute error distribution was 0.15 or less. If we assume a uniform probability distribution on the interval (0, 1) for the previous distribution of p then it is well known that the posterior distribution of p is a β distribution.^29,34  We estimate the probability that the CLAD system is reliable across all levels as the mode of the posterior distribution. We consider the experiment to have established reliability of our CLAD system across levels if 0 is less than the left endpoint of the 95% Bayesian credibility (confidence) interval for p.

We performed similar Bayesian analyses to assess the accuracy of the CLAD system across levels. We consider the experiment to have established accuracy of our CLAD system across all levels if 0 is less than the left endpoint of the 95% Bayesian credibility (confidence) interval for P.

We chose 18 levels because if there were a minimum of 14 levels of reliable (accurate) control giving a posterior probability of estimate of 0.78, the lower endpoint on the 95% Bayesian credibility interval would be 0.55, meaning that overall it was more likely that the CLAD system was reliable (accurate) than not.

Using our CLAD system, we successfully segmented the electroencephalogram into bursts and suppressions in real time for each animal. The system was robust to electroencephalogram signal quality and the overall morphology of burst activity (fig. 2). In some animals, we encountered good electroencephalogram signal quality and sharp, easily discernible bursts (rat 6; fig. 2A) whereas in others, we found noisier electroencephalogram signals and broader bursts (rat 1; fig. 2B). In each case, we were able to use the design parameters of our electroencephalogram filtering (Eqs. 19–20 and panel 1 in fig. 2, A and B) and thresholding (Eq. 21 and panel 2 in fig. 2, A and B) to obtain good segmentation of the electroencephalogram into the binary time series (panel 3 in fig. 2, A and B). In these representative examples, we illustrate the ability of the filtering, governed by the forgetting factor α (Eqs. 20 and 21), to distinguish the amplitude envelope of bursts from the background suppression level. Similarly, we show that the amplitude threshold v_threshold partitioned bursts from suppressions assigning a 0 (1) when the filtered electroencephalogram was less than (exceeded) v_threshold. The forgetting factors were 0.995 for animals 1 and 2, 0.595 for animals 3 and 4, and 0.695 for animals 5 and 6 whereas, the thresholds were 3 × 10⁻⁵ μV for animals 1 and 2, 1.5 × 10⁻⁶ μV for animal 3, and 4.3 × 10⁻⁶ μV for animals 4, 5, and 6. A larger forgetting factor (more forgetting) corresponds to more filtering, meaning that the animal had a noisier electroencephalogram (e.g., fig. 2B), whereas a smaller forgetting factor (less forgetting) is more suitable for easily discernible bursts, meaning the animal had sharper electroencephalogram signals (fig. 2A).

In each animal, we were able to effectively record the dynamic response of the electroencephalogram to infusions of propofol by applying the BSP filter to the binary time series. The BSP responded reliably in all six animals to bolus infusions, which we demonstrate with a representative example (fig. 3). There were clear changes in the unprocessed electroencephalogram (fig. 3A), the segmented electroencephalogram (fig. 3B), and the corresponding BSP time course (fig. 3C) estimated by the BSP filter algorithm (Eqs.7–13). This bolus rapidly brought the animal to a nearly isoelectric electroencephalogram, as evidenced by a long period of suppression (fig. 3A). Consequently, the BSP rose rapidly (fig. 3D). As the effect of the bolus subsided, bursts reappeared, were detected by our segmentation method, and the BSP began a slow decay. The success of the BSP algorithm allowed us to conduct the system identification (fig. 4) and to estimate in real time the animal’s instantaneous state of burst suppression for the CLAD system (fig. 5).

We completed successful system identification in each animal, estimating the parameters of its pharmacokinetics model in response to one or more propofol boluses. To illustrate, we present a representative example (rat 1) with two boluses of propofol (fig. 4B). We administered the second bolus once the electroencephalogram returned to a continuously active state, that is, all suppression had subsided. The second propofol bolus produced a similar BSP response as the first bolus (fig. 4A). After this second bolus, we fit the pharmacokinetics model by nonlinear least squares to the time course of effect-site concentrations computed from the BSP estimates using Eq. 12 (see the appendix). The fit accurately described the time course of the BSP in response to the bolus sequence, demonstrating that our pharmacokinetics model captured well the dynamics of the state of burst suppression induced by propofol (fig. 4A). In particular, we found that the second-order state models are sufficient to account for the rodents’ BSP responses, obviating the need for more detailed four-compartment pharmacokinetics models³⁶  commonly used to represent propofol.

We performed all system identification for each animal online, whereby the pharmacokinetic model fitting did not depend on the specific shape, size, or design of the bolus sequence. In each animal we used two or three boluses. After estimating the model parameters, we computed each animal’s parameters for its PI controller (Eqs. 17–18), which we then used for closed-loop tracking.

For each of the six animals the CLAD system tracked the target levels closely by making control corrections every second (fig. 5). The PI controller rapidly changed the infusion rates to maintain control at the target levels. The infusion rates changed most to maintain the BSP level of 0.9, and least to maintain the level of 0.4. The larger changes at 0.9 were expected because a BSP of 0.9 means administering larger amounts of propofol to keep the electroencephalogram isoelectric 90% of the time yet, allowing for burst activity 10% of the time. The controllers made transitions from a lower to a higher target level by increasing the infusion rate, whereas they made the transitions to lower rates by decreasing or frequently, setting the infusion rate to zero. The controller responded rapidly to correct divergence between the estimated and target BSP levels. In two cases (fig. 5, D and E) in which the BSP estimate showed transient excursions from the target trajectories, the controller compensated immediately and restored control. In all cases, the computer regulated the infusion rates with second-to-second dynamics that could not be equaled by a human operator.

In agreement with the plots in figure 5, the magnitude of the tracking errors as evaluated by the SD, and the MADs were small (table 1). For the target levels 0.4, 0.65, 0.9 and the between-level transitions, the SDs (MADs) were respectively 0.039 (0.34), 0.062 (0.4), 0.042 (0.026), and 0.057 (0.56). The smaller values of the MAD relative to the SD and the percent error relative to the median absolute percent error are expected, given the insensitivity of the MAD to large errors. The errors were approximately equal on levels 0.4 and 0.65, lowest on level 0.9, and highest during the transitions. The percentage errors relative to the SD and the MAD, as well as the bias, were consistent with this observation.

The better control at the 0.9 target level is most likely because the BSP state of 0.9, i.e., a probability of 0.9 of being suppressed, is easier to estimate from the electroencephalogram recordings than the BSP states of 0.4 and 0.65. The variance of a binomial random variable (Eq. 4) is lowest for p_t closer to either 0 and 1, and highest for p_t close to 0.5. The larger errors during the transitions likely reflect anticipated limitations in the controller performance (see Discussion).

Our statistical analysis (see Materials and Methods) shows that on 17 of the 18 levels the 95th percentile of the absolute error distribution was less than 0.15 (fig. 6A). The exception was animal 2, at level 0.65 for which the 95th percentile was 0.155. The posterior probability that the CLAD system was reliable across all levels was 0.94 = 17/18 (95% CI, 0.77–1.00; n = 18). Because zero is well below the left endpoint of these Bayesian credibility (confidence) intervals we conclude that our CLAD system is reliable.

On 18 of the 18 levels, the 95% CIs for the errors included 0, whereas on 17 of the 18 levels, the 50% CIs included 0 (fig. 6B). For the latter, the only exception was animal 5 at level 0.4. The posterior probability that the CLAD system was accurate across all levels was 1.00 =18/18 (95% CI, 0.84–1.00; n = 18). Similarly, the posterior probability that the CLAD system was highly accurate across all levels was 0.94 = 17/18 (95% CI, 0.77–1.00; n = 18). Therefore, because 0 is well within these CIs we conclude that our CLAD system is highly accurate.

To study the feasibility of automating control of medically induced coma, we developed a CLAD system to control burst suppression in a rodent model. We demonstrated that our CLAD system can reliably and accurately control burst suppression in individual animals across dynamic target trajectories.

Closed-loop anesthetic delivery system development started more than 60 yr ago¹⁰  and later reappeared in the 1980s.²⁴  There have been several clinical studies of CLAD systems, and versions are now commercially available.³⁷  The most frequently used control signal has been the Bispectral Index (BIS).^{7–9,11,13,14,17,18,21–23,38–40} Other control signals have included a wavelet-based index,¹²  entropy measures,¹⁶  an auditory evoked potential index,¹⁵  and the spectrogram median frequency.^24–26  These systems have been constructed with standard and nonstandard control paradigms^{15,16,23,24,40}  and used principally to control unconsciousness.^{7,9,12,14,16,18,23,24,38,39}  A recent report investigated control of both antinociception and unconsciousness.¹⁶  The criteria for successful control differed across these studies. Schwilden et al.²⁴  demonstrated control of median frequency in individual human subjects. In contrast, several of the studies that used BIS as the control signal defined successful control as a BIS value between 40 and 60, and reported BIS time courses averaged across subjects.^{7–9,11,13,14,17,18,21–23,38–40} Control using BIS as a control signal is achieved with a 20- to 30-s delay required to compute the BIS updates.⁴¹  In contrast, our system updates the control input every second with no delay. Those CLAD systems that have been developed to study burst suppression have also only shown results for time courses averaged across subjects.^27,28  None of these studies considered control of dynamic trajectories nor conducted a formal statistical assessment of reliability and accuracy.

Our work makes important improvements on current CLAD systems. We chose burst suppression as a control state because, unlike the brain states defined by the BIS score,⁴²  the state of the brain in burst suppression is well-defined neurophysiologically.^4,5  Furthermore, burst suppression has a well-defined electroencephalogram signature that can be quantified in real time, and therefore, controlled. Our two-dimensional pharmacokinetics model (Eqs. 1–2) provides a simple and sufficient representation for capturing the essential properties of burst suppression. This model was the starting point for designing our CLAD system. We formulated this model based on observations made while monitoring the electroencephalogram of patients under general anesthesia in the operating room. We noticed that once the state of burst suppression is achieved increasing or decreasing the rate of a propofol infusion directly increases or decreases the rate of suppression events. However, because intravenous injection of propofol does not deliver the drug directly to the brain, a burst-suppression model must have a minimum of two compartments. Therefore, for control of burst suppression, our simpler second-order model can replace more detailed four-compartment models³⁶  because the objective is to control a single brain state and not the wide range of brain states that could be represented by the higher-order model.

The segmentation and BSP filter algorithms and our system identification procedure were critical for implementing our CLAD system in real time. Tuning the segmentation algorithm allowed us to robustly detect bursts over a broad range of signal qualities. The BSP filter algorithm, derived from our state-space estimation paradigm for point processes and binary observations,⁴³  uses the one-to-one relationship between the BSP (Eq. 3) and the effect-site anesthetic level (Eq. 12) to estimate in real time the effect-site anesthetic level. We have previously shown that the BSP filter algorithm gives more credible burst-suppression estimates than the burst-suppression ratio.³¹  The burst-suppression ratio can require up to 5 min to estimate the brain’s burst-suppression state,⁴⁴  a feature that would substantially limit its use in real-time control.

The system identification procedure (fig. 4), allowed us to estimate model (Eqs. 1–2) and control parameters (Eqs. 14–18) for each animal and thereby implement individually tailored PI control strategies. We chose a PI controller because it is widely used and known to be robust noise and parameter uncertainty.⁴⁵ 

We used a statistically efficient design in our experiments by testing control at all three levels within animal. Furthermore, by allowing 5–10 min for the CLAD system to transition between levels, data from different levels within animal are effectively independent. Hence, the 18 levels served as the units of analysis instead of the six animals in our assessments of overall reliability and accuracy.

The final novel feature of our approach is the use of specific statistical criteria to assess reliability and accuracy of our CLAD system within levels and across the entire experiment. We adapted accepted concepts in reliability theory to define these criteria and computed the overall assessments of reliability and accuracy in a Bayesian framework.³³  Current performance measures for CLAD systems have been adapted from those used to assess performance of target-controlled infusion systems.³²  Recently, CLAD systems have been evaluated by comparing their performance to performance using manual control.¹⁹  Our Bayesian paradigm should facilitate design and testing of future CLAD systems by making it possible to assess performance in terms of specified properties on a system’s error distribution.

Many technical improvements can be made in our CLAD system. The second-order state model performed well in our PI controller. A more detailed state model could be constructed from our recently developed neurophysiological metabolic model of burst suppression.⁵  As an alternative to our deterministic PI controller, we could model system noise explicitly and apply a stochastic control strategy.⁴⁶  A model predictive control strategy could be adopted to formally impose constraints such as nonnegative infusion rates.⁴⁷  The BSP filter algorithm (Eqs. 7–10) could be improved by using instead of the first-order random walk model (Eq. 6), a stochastic version of our two-dimensional state model (Eq. 1) to estimate simultaneously the BSP and its rate of change. This is equivalent to estimating both the peripheral and effect-site anesthetic levels from the binary time series. Modifying the BSP filter in this way could improve tracking performance during transitions by allowing more rapid increases in transitions from lower to higher target levels (fig. 5A) and preventing undershoot in transitions from higher to lower levels (fig. 5D). By doing so, we would be adopting the common practice of using the same model for state estimation and control. Tight control during transitions was not a primary design consideration in the current work because large and frequent level changes are not usually required to manage medically induced coma.

There are several possible benefits of using a CLAD system to control medically induced coma using burst suppression. These include maintaining tight control of brain states for extended periods, providing adequate brain protection with the least amount of anesthetic, facilitating periodic arousals for neurological assessments, reducing the likelihood of propofol overdose syndrome,⁴⁸  and using intensive care unit staff more efficiently.

To realize these benefits will require development of our CLAD system in humans. Because in the intensive care unit the controller will have to function over several hours or days, it would be prudent to next test our CLAD system in a rodent model for periods longer than an hour. Experiments lasting for several hours will require changes in our current protocol to allow intubation and mechanical ventilation, as well as invasive monitoring of blood pressure and use of vasoactive drugs. Though challenging, these experiments would substantially test system robustness and better approximate actual management of intensive care unit patients using our system. System identification in the intensive care unit would be conducted when the initial doses of propofol are administered to induce burst suppression. This approach would obviate the administration of propofol solely for parameter estimation, a maneuver that could further destabilize a hemodynamically unstable patient. Over several hours or days, the assumption that the model and control parameters remain constant will likely not hold. If the model and control parameters were to drift then, we would estimate them adaptively.⁴⁹  This is a tractable yet, nontrivial estimation problem.

In conclusion, we have demonstrated that a CLAD system using a computer-controlled infusion of propofol can reliably and accurately control burst suppression in a rodent model. Our work opens the possibility of implementing a CLAD system to control burst suppression for maintenance of medically induced coma in intensive care unit patients, and eventually, the anesthetic state of patients in the operating room.

To construct a closed-loop anesthesia delivery system we specify four components: the state model, the state estimation algorithm, the system identification procedure, and the controller. Here, we summarize the essential details.

We assume that two is the minimal dimension of a state model required to define the brain’s state of burst suppression in response to administration of an anesthetic. We construct our state model by adapting a two-compartment pharmacokinetics system for intravenous drug infusion. We let x(t) = (x_1t, x_2t) be the state of the system at time t where x_1t is the amount of anesthetic in the peripheral compartment, and x_2t is the amount of the anesthetic in the brain or the effect-site compartment. We let I(t) denote the infusion rate of the anesthetic at time t. We assume that the anesthetic enters into the peripheral compartment; the anesthetic flows back and forth between the peripheral and the effect-site compartments; the anesthetic is eliminated from the body only through the peripheral compartment; and the amount of anesthetic in the effect-site compartment determines the electroencephalogram level of burst suppression.

The differential equation defining this state model (fig. 1) is

where

The response of the brain to the anesthetic is monitored by the electroencephalogram. We assume that sufficient doses of an anesthetic are administered to a subject to induce burst suppression and that the effects of the anesthetic can be observed continuously in time by recording the electroencephalogram. We further assume that the electroencephalogram can be filtered and thresholded in real time to identify burst and suppression events at a resolution of Δ. To link the electroencephalogram to the state of burst suppression, we let n_t be the binary time series constructed from the filtering and thresholding, where n_t= 1 if there is a suppression at time t, and n_t= 0 if there is a burst at time t (Eqs. 19–21). We define the burst-suppression probability (BSP) as

Because p_t defines the probability of a suppression event at time t given x_2t it follows that n_t obeys the Bernoulli process

Equation (3) is a hyperbolic transformation that maps the amount of the anesthetic in the brain, defined on the interval (0, ∞) into p_t, a well-defined probability on the interval (0, 1). In this way, p_t provides an instantaneous output of the probability of the brain being suppressed.

To implement the CLAD system we require a way to estimate p_t or equivalently, the brain state x_2t, from the binary time series n_t. We develop a version of a binary filter algorithm^31,43  to compute estimates of p_t and x_2t in real time. We assume a simplified, stochastic version of the state model in equation 1 by taking

and assuming that z_t obeys the Gaussian random walk model

where the v_t are independent, zero mean Gaussian random variables with variance . The transformation in equation 5 ensures that x_2t remains nonnegative. Given estimates of z₀ and , the following binary filter algorithm can be applied to the n_t to compute p_t and x_2t.^31,43  It is

where

for t = 1, ..., T and the notation z_ts denotes the estimate of z_t, given the data up through time s. It follows from equations 3, 7, and 9 that at time t the estimates of x_2t and p_t are, respectively,

We term equations 8–13 the BSP filter algorithm.

System identification for our CLAD entails estimating z₀ and , the parameters of the BSP filter algorithm and A, b, the parameters of the state model, in a two-step procedure. First, we assume that a preliminary experiment is conducted in which a bolus dose of the anesthetic sufficient to induce burst suppression is administered to the subject and the electroencephalogram is converted into the binary time series by filtering and thresholding (see Electroencephalogram Segmentation Algorithm below). We estimate z₀ and from the binary times series derived from the bolus experiment by applying an approximate expectation maximization algorithm^31,43  to the state space model defined by equations 4 and 6. The expectation maximization algorithm also provides z_tT the estimate of z_t given all of the binary observations in the bolus experiment, and as a consequence, by equation 12, x_2tT the estimate of x_2t. In the second step, we use the estimated state, x_2tT as data to estimate A and b by nonlinear least squares.

If p_target is the target level of burst suppression, then it follows from equation 3 that the corresponding target effect-site concentration of the anesthetic is

and hence, that the error signal for our controller at time t is

Our objective is to construct a proportional-integral controller of the form⁴⁵ 

where is u(t) is the control signal at time t, t₀ is the start time of the control interval and a_p, a_i are control parameters to be determined. If we take as our design criterion the implementation of a proportional-integral controller that achieves a fast rise time up to a specified level of burst suppression while minimizing the overshoot then it follows from standard control theory arguments that⁴⁵  we take

and

where the coefficients a₁₀, a₁₂, a₂₁, and b are defined in equation 2. Equations 17 and 18 show that once we have estimated the parameters of the pharmacokinetics model, the parameters for the controller are completely defined.

In this control problem we observe the effect-site concentrations x_2t through the binary time series n_t. We convert the electroencephalogram signal y_t into n_t using the following algorithm. At each observation time we compute the following time-varying mean and variance, and evaluate the threshold criterion

where α is a forgetting factor between 0 and 1 and v_threshold is a threshold voltage we set to define a burst. A value of α closer to 0 (1) corresponds to less (more) forgetting. The algorithm in equations 19–21 tracks and s², the time-varying mean and variance, respectively. If the time-varying variance exceeds the threshold, then the electroencephalogram is in a burst and n_t= 0, whereas if the time-varying variance does not exceed the threshold, then the electroencephalogram is in a suppression and n_t= 1.