Animal studies were approved by the Subcommittee on Research Animal Care, Massachusetts General Hospital, Boston, Massachusetts, which serves as our Institutional Animal Care and Use Committee. 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.

BMI Design

Overview. We use a stochastic optimal control paradigm to design a real-time BMI to control medical coma using burst suppression (Figure 1a). As our measure of the burst suppression level, we use the burst suppression probability (BSP), a number between 0 and 1, which defines the instantaneous probability of the EEG being suppressed. The BSP is computed in one-second intervals in real time by filtering and thresholding the EEG to convert it into binary observations (Figure 1b). To estimate the BSP from the binary observations we first formulate a two-dimensional compartment model that relates the BSP to the concentrations of the anesthetic in the central compartment and the effect site compartment (Figure 1c). We next estimate the parameters of the compartment model based on the EEG observations recorded in a systems identification experiment conducted prior to initiating real-time control (Figure 2). We carry out our stochastic control framework by developing from the two-dimensional compartment model a recursive Bayesian estimator of the concentration states and consequently of the BSP from the binary observations in real time (14)–(17). We develop a LQR controller that takes the concentration estimates as feedback and determines the drug infusion rate in real time (25). In addition to the LQR control strategy, we also implement a model predictive controller (MPC) that allows us to explicitly impose constraints on the anesthetic infusion rates (27). We present the mathematical details of the system identification, formulation of the Bayesian estimator and the two controllers for the interested readers below. These mathematical details in this subsection are not necessary to follow the remainder of the paper beginning with the Results. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Figure 1. The BMI system. (a) The BMI records the EEG, segments the EEG into a binary time-series by filtering and thresholding, estimates the BSP or equivalently the effect-site concentration level based on the binary-time series, and then uses this estimate as feedback to control the drug infusion rate. (b) A sample burst suppression EEG trace. Top panel shows the EEG signal, middle panel shows the corresponding filtered EEG magnitude signal (orange) and the threshold (blue) used to detect the burst suppression events, and bottom panel shows the corresponding binary time-series with black indicating the suppression and white indicating the burst events. (c) The two-compartmental model used by the BMI to characterize the effect of propofol on the EEG. https://doi.org/10.1371/journal.pcbi.1003284.g001 PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Figure 2. System identification. (a) and (b) show two sample fitted system responses. The measured BSP trace in response to a preliminary bolus of propofol is shown in grey and the response of the second-order system model in (2) fitted using nonlinear least-squares is shown in red. https://doi.org/10.1371/journal.pcbi.1003284.g002

Problem formulation. Our goal is to control the anesthetic state of the brain in burst suppression, which depends on the effect-site (i.e., brain) drug concentration. The burst suppression state or the effect-site concentration, however, are not directly observable. What we observe is the EEG signal, a stochastic process that depends on the burst suppression state. To design the closed-loop BMI, we present a certainty-equivalent optimal feedback control approach [8] by deriving an estimator for the burst suppression state based on the EEG observations and designing an optimal feedback controller that takes this estimate as a feedback signal to control the drug infusion rate in real time (Figure 1a). As our measure of the burst suppression state, we use the burst suppression probability (BSP) by filtering and thresholding the EEG signal in small intervals to identify the activity in each interval as a burst or a suppression event (see Experimental Procedure; Figure 1b). BSP is then defined as the brain's instantaneous probability of being in the suppressed state at a given time interval. We denote the BSP at time by . The BSP is in turn related to the effect-site drug concentration. Since higher levels of effect-site concentration should result in higher levels of BSP and since BSP should be a number between , in this work we relate the BSP to a measure of the effect-site concentration, , using a hyperbolic transform (1)Hence our goal is to control the BSP, or equivalently to control this measure of effect-site concentration, . To develop the estimator and the controller, we construct a state model for the drug concentration state that describes its dynamics in response to propofol infusion. Pharmacokinetic models characterize the dynamics of a drug's absorption, distribution, and elimination in the body (e.g. [9], [10]). We adapt a simplified two-compartment linear pharmacokinetic model [11] to describe the anesthetic drug's dynamics in burst suppression. In this model, one compartment represents the central plasma and the other compartment represents the effect-site or brain (Figure 1c). The anesthetic drug enters the body and is eliminated from the body through the central compartment, and can flow in both directions between the two compartments. In the Results section we show that this model is sufficient to achieve reliable and accurate control of burst suppression. Given the two compartments in the model, the concentration state is two-dimensional and is denoted by , where as before, is the brain's anesthetic concentration and is a measure of the central plasma concentration at time . Denoting the sequence of drug infusion rates by , the sequence of anesthetic concentration states in the two-compartment model, , are generated according to the linear dynamical system (2)where (3) (4) (5)Here is the discretization time step, and (or equivalently , , and ) are parameters of the two-compartment model that we need to estimate (Figure 1c). We estimate this model for each animal from the EEG data prior to initiating real-time control as discussed in detail in the System Identification section. We first derive a recursive Bayesian estimator of the burst suppression level from the EEG thresholded and segmented into a binary time-series. We then derive an optimal feedback-controller that uses this estimate as a feedback signal to decide on the drug infusion rate in real time.

Recursive Bayesian estimator. We now develop a recursive Bayesian estimator for the drug concentrations and consequently for the BSP based on the binary observations of the thresholded EEG signal. Since the drug concentration state, , is a positive variable, we estimate its logarithm, , from the EEG signal instead. A recursive Bayesian estimator consists of two probabilistic models: the prior model on the time sequence of the concentration states, and the observation model relating the EEG signal to these states [12], [13]. Using the two compartment model in (2), we write the prior model for as (6)where is a zero-mean white Gaussian noise with covariance matrix and summarizes the uncertainties in the state model, and is the th component of . At time , is the known drug infusion rate used by the BMI in the previous time step. Note that our prior state model is nonlinear. The observation in the estimator is the binary time-series of the burst suppression events obtained by thresholding the EEG (see Experimental Procedure; Figure 1b). To construct the observation model, we assume that in each time interval there can be at most suppression events and that the number of such suppression events is binomially distributed with burst suppression probability . Denoting the number of suppression events by , the observation model is given by (7)where we have indicated the dependence of the BSP, , on the states, , explicitly. Using the prior and observation models in (6) and (7), we now derive the recursive Bayesian estimator. The estimator's goal is to causally and recursively find the minimum mean-square error (MMSE) estimate of the state at each time step, which is given by the mean of the posterior density at that time step, . To derive the recursions and denoting the sequence of suppression counts by , using the Bayes rule we can write the posterior as (8)which states the posterior density as a function of prediction density, . Note that we have used , since we assume that the observations of the EEG at a given time step only depend on the concentration states at that time step and hence are conditionally independent of the previous EEG observations. Using the Chapman-Kolmogorov equation, we can in turn write the prediction density as (9)Here we have used the conditional independence, , which comes from the prior model in (6). Now the second term inside the integral is just the posterior density from the previous time step. Hence substituting (9) into (8) generates the recursion. The exact expression in (8) is in general complicated and not easy to find analytically. In the special case when both the prior and observation models are linear and Gaussian, these recursions have exact analytical solutions and result in the celebrated Kalman filter. In our case, however, first, the prior model is nonlinear and second, the observation model is not Gaussian but binomial. Hence we make two approximations at every time step to compute the recursions. First, similar to the case of the extended Kalman filter, we make a linear approximation to the prior model at each time step. Second, we make a Gaussian approximation to the posterior at each time step. We denote the mean of the posterior, i.e., , by and its covariance matrix by . Similarly, we denote the mean of the one step prediction density, , by and its covariance matrix by . As the first approximation, we linearize the prior model in (6) around the posterior mean, . Doing so we have (10)where (11) (12)with denoting the evaluation of the inside expression at value and with (13) As the second approximation, we make a Gaussian approximation to the posterior density. Doing so, from (9) the prediction density will be approximately Gaussian since is approximately Gaussian from (10). Using (10) we can find the prediction mean and covariance as (14) (15)This is the prediction step of the estimator. Now making the Gaussian approximation we get the update step (see Supporting Text S1 for details) (16) (17)where again indicates the evaluation of the inside expression at and (18) (19)with (20) Hence (14)–(17) give the estimator recursions. The estimator finds the MMSE estimate of the state or equivalently the posterior mean at time in two steps: first, before data is observed, it uses the prior model in (6) to make a prediction on the state, i.e., find given and —this is the prediction step in (14) and (15). Once data is observed, it combines the observation model in (7) with the prediction density to find the posterior mean —this is the update step in (16) and (17). Consequently since , we find the concentration state estimate as , and since the BSP is related to by a hyperbolic transform in (1), we estimate it as .

Optimal feedback-controller. The recursive Bayesian estimator derived above provides us with a real-time estimate of the concentration states at each time step. We now design a real-time optimal feedback-controller that takes as feedback this state estimate and decides on the sequence of drug infusion rates to control the BSP. To find in the optimal control framework, we need to quantify the goal of the controller in a cost function that will then be minimized by selecting the optimal . For the linear state model in (2), if we pick the cost function as a quadratic function of the state and control signals given by (21)where is the time duration of anesthesia, and are positive semidefinite and is positive, then the optimal control signal at any time, , is simply a linear feedback of the state at that time given by [8] (22)where the feedback matrices, , can be found recursively and offline [8]. This is the linear-quadratic-regulator (LQR) solution. Moreover, when the state model is controllable (as is the case in our problem using the experimental fits; see Results), there exists a steady-state solution, , for the feedback matrix in (22). This steady-state feedback matrix is the solution to the discrete form of the famous algebraic Riccati equation given by [8] (23)where (24) In the general LQR formulation above, however, the goal is to derive the states close to zero—while limiting the total amount of control—as evident from the cost function in (21). In the control of burst suppression, our goal is to achieve a desired non-zero target BSP level, , or equivalently to take the effect-site concentration state close to a non-zero target level , using as little drug as possible. Hence to find the solution, we additionally shift the origin of the state-space to [14]. This way, the control goal is equivalent to deriving the shifted state variable close to zero, as in the classical LQR formulation. We show in the Supporting Text S1 that, in our problem, it is possible to shift the origin and the optimal drug infusion rate is in turn given by (25)where (26)The value of at each time step is provided by the estimator. Note that the LQR formulation does not impose any constraints, such as positivity of the control variables. In practice we can impose these constraints by bounding the LQR control solution in (25) appropriately (for example if the solution is negative, use zero instead). Another way to solve optimal control problems with constraints is to use a model predictive control approach as we develop next.