Robust Control of a Biophysical Model of Burst Suppression

Burst suppression is a brain electrical phenomenon measured on the electroencephalogram (EEG) that is typically associated with a state of general anesthesia [1], as well as cooling [2] and certain types of coma [3]. It is characterized by bi-phasic dynamics when electrical activity alternates between high voltage activity (bursts) and low voltage activity (suppression). When burst suppression occurs, the suppression intervals are largely devoid of frequency content, while bursts often manifest physiologically meaningful oscillatory power in the 8–30 Hz range [1]. Figures 1(B) and 1(C) from Ref. [4] show a comparison of EEG data while a patient is under general anesthesia but has not displayed burst suppression (B) and then a deeper general anesthesia wherein burst suppression has been induced (C).

A key aspect of burst suppression is that the relative proportion of the EEG signal spent in suppression varies with respect to the depth of anesthesia. Thus, when pharmacologically inducing burst suppression, one can increase the “amount” of suppression by increasing the dosage applied to the given patient. This can be visualized in Fig. 2.

Such phenomenology is of interest clinically because burst suppression is often targeted as part of a medically induced coma for patients that have suffered brain injuries [5]. In this situation, clinicians often attempt to titrate to a specific burst suppression ratio (BSR) “level.” However, such practice is often quite imprecise since clinicians rely on visual assessments of the EEG and relatively infrequent, bolus-type dosing of anesthetic drugs. As a result, interest has grown in the use of closed-loop, continuous systems for controlling BSR [6–9].

There have been many attempts at formulating a solution to this problem [6–9]. These frameworks follow a general framework wherein a diffusion model transforms anesthetic infusion to a drug brain concentration variable (i.e., the pharmacokinetics (PK)). The concentration is in turn mapped either through a static transformation or via a probabilistic decoder, into a BSR estimate. In the current paper, we strive to further the progress made in the robust control of burst suppression by introducing a biophysical, nonlinear dynamical systems model of burst suppression in place of the latter transformation. Furthermore, we introduce actuator dynamics into the closed-loop anesthetic delivery (CLAD) system. In so doing, we are able to introduce uncertainty not only in the patient's physical characteristics (such as patient height, weight, age, and sex) but also in the neural circuit dynamics of the patient. Once coupled with the presence of an actuator (which would simulate an anesthetic pump), we are able to study the robust control design of CLAD for burst suppression.

Our basic modeling framework is built on the canonical Wilson–Cowan neural mass model [10,11], which we have recently augmented to model burst suppression as a result of fast–slow dynamical interactions attributed to neural and metabolic processes [12]. This model provides us with a mathematical basis by which we can study the physiological processes underlying burst suppression.

The specific contributions of this paper are as follows:

1. Introducing a closed-loop simulation architecture with actuator dynamics, pharmacokinetic dynamics, and neurological dynamics.

2. Design of a robust control law for controlling set-point BSR levels despite parameter uncertainty in the pharmacokinetic model and the Wilson–Cowan model.

We begin by providing the model of burst suppression that we will use in order to pursue the two specific contributions noted above. To this end, we describe an existing dynamical systems model from the literature, which will enable our analysis and design.

The Wilson–Cowan equations constitute a mean-field neural model that describes aggregate activity in a large population of neurons [10,11]. The primary mechanism of the model involves the dynamics of excitatory and inhibitory neural populations. In Ref. [12], the basic Wilson–Cowan formulation was extended to modify these equations to give them a modulating input, ϕ_j(t), in addition to their baseline fast dynamics. Within these modified Wilson–Cowan equations [12] there are parameters that describe the recovery dynamics of these neurons, which govern how quickly the neurons transfer from a suppressed state back to a bursting state. By varying these parameters, an internal neurological uncertainty can be created and used to increase the robustness of burst suppression feedback control algorithms. Also within these equations, there are various parameters that can be used to quantify the concentration of a pharmacological substance, such as propofol, present in the brain effect site. We can treat the concentration of propofol as the control input to increase the BSR to the desired state. It is crucial to note that clinically, the level of concentration in the effect site cannot be decreased except by a metabolic process. The control input, a pharmacological drug, can only be used to increase the ratio. Therefore, it is vital to both correctly estimate the BSR and use a control method that does not produce high overshoot or steady-state error. These equations are thus presented and described:

The recovery and consumption processes are primarily a function of the excitatory dynamics. Our main focus will be on the evolution of the variable β, and, more specifically, how small changes in the constants that describe its behavior are crucial in how the system as a whole behaves.

Finally, it should be noted that the main variables e and i should not be misunderstood as EEG activity, although there are fundamental similarities between the two. This distinction was made in Ref. [12] while also stating that it is a suitable variable for burst–suppression type studies. Therefore, in this paper, we will be analyzing the excitatory firing rate as a surrogate for EEG data.

To succinctly summarize the equations, variables, and parameters above, Tables 2 and 3 are given in the  Appendix. These tables provide descriptions of the state variables and the parameters as well as the initial conditions and nominal values according to Ref. [12]. Additionally, in Ref. [12], a time scale of milliseconds was used. Since we will be studying the model in the span of hours, we have proportionally increased the time constant and metabolic recovery and consumption rates in order to allow the dynamics to scale properly to a timescale of seconds.

The purpose of using this model is to take advantage of its inherent burst suppression qualities as the parameter c₂ is varied. We can use this parameter to model the effects of propofol or other inhibitory agonists on the dynamical model. In a clinical setting, a patient will become more sedated as the propofol infusion rate is increased. Therefore, it is vital to describe bounds with which the model will behave properly as c₂ is varied. In Ref. [12], the bounds were set by the system's stability about the c₂ parameter. The system has a bifurcation point around the value of 7.7 after which the system develops a stable limit cycle envelope where burst suppression behavior occurs until a c₂ value of 69. At this value, a stable steady-state solution can be found, and the model becomes completely suppressed. Therefore, we wish to analyze how burst suppression emerges and changes as the value of c₂ increases from 7.7 to 69. Figures 3 and 4 give an example for how the model behaves with multiple c₂ values.

The values of c₂ were chosen specifically to show how the bursting amplitude changes with time. It is also notable within the spectrogram how when the model reaches near-complete suppression, frequencies in the range above 8 Hz are no longer present. These differences in the model expression will be considered when the burst suppression estimation algorithm is presented and used in later sections.

It is important to emphasize that our goal is not to validate the above model since this has already been carried out in the prior literature. Rather this model forms the basis of the control analysis and design that is our primary contribution.

In addition to simply running a propofol (c₂) increase on this model, we wish to generate some uncertainty for different patient types. We focused primarily on the variation of the recovery parameters, $kβ,ζ$⁠, and ν to simulate how different patients may respond to propofol as it enters their system. Figure 5 shows the same c₂ value of 30, but with a change in these recovery parameters. The changes in these recovery parameters cause a notable difference in the suppression length for a given value of c₂. Each simulation was performed with each parameter changed independently, while the others remained at their nominal value.

This shows that the system is much more sensitive to a change in either $kβ$ or ζ, while a significant change must be made to ν to see a similar percent change in suppression time. This uncertainty will be used to analyze the effectiveness of our control architecture during the simulation results of this paper.

As described in the Introduction, the burst suppression ratio quantifies the level of suppression manifest in the EEG of a patient. It is measured as the ratio between the time the signal is suppressed over the total time measured. We constructed a method to estimate the BSR from our model output. The algorithm uses the defining frequency characteristics of burst suppression for BSR estimation. In the simulation studies that follow, we set the model output sample rate (512 Hz) that would emulate sampling commonly used for EEG processing.

A positive value for the exponent was chosen because with lower ranges of c₂, higher values of m_s and m_b are typically found, thus giving a value of S(n) closer to zero which indicates a burst event. In order to make S(n) suitable for a wide range of c₂, k_sb was made equal to 1, x_0sb was made equal to −30, and k_b was given a value of 500.

With this parameterization, a nearly sigmoidal relationship was found between the BSR and c₂. Figure 6 shows the BSR calculated as a function of c₂. This algorithm calculated a higher value for the lower ranges of c₂, specifically between 8 and 10. After a c₂ value of 10, the BSR increases in a near monotonic fashion until it reaches a near 1 BSR at a c₂ value of 69.

Despite these limitations on the lower ranges of c₂, the remainder of the values with this parameterization produced results consistent with what would be expected for burst suppression as the propofol concentration is increased. Since this is the type of physiological behavior the model is attempting to reproduce, we can lower bound the effective value of c₂ where the BSR stops decreasing and begins its monotonic increase to 1.

In order to characterize the suppression of the system over time, a running mean of S(n) is calculated. This results in a delay before a certain confidence level in the value can be achieved. This is useful information from a control design perspective because these algorithms are very nonlinear in nature and we desire to implement a linear control method. Thus, if we can find the expected delay values for this algorithm to converge, our control design can account for it by using a linear approximation for the delay associated with time to converge to a confidence level. Since we desire tight control on our BSR to eliminate fluctuating the patient between multiple BSR levels, a convergence value for the mean signal variance was chosen to be 0.0001.

In order to design the gains for the system, a time for this algorithm to converge had to be selected. Thus, the time for this algorithm to converge with nominal and varied recovery parameters of a limited subset to that of the range run in Sec. 2.1 were calculated. The results for these convergence times are shown in Fig. 7.

This variation on convergence times shows some very sharp increases in the lower ranges of c₂ but very quick convergences for most values of c₂. This is likely due to the complete absence of higher frequency bursting data at high c₂ values, which make it much easier for the algorithm to converge given the lower frequency data are more consistently defined and are more present at these points. Nonetheless, a delay for this algorithm to be implemented in linear gain design was chosen to be the average of the values above which came out to be approximately 0.7 s.

We now proceed to enact our actual closed-loop strategy, using the BSR estimation procedure described above as a feedback signal to a controller. To facilitate design, we used a classical control system architecture designed around the Wilson–Cowan dynamics in order to more accurately simulate a CLAD system. All prototyping and subsequent simulations were performed in Simulink and matlab 2019a. This model was used in Rapid Accelerator mode to be able to run 18 h of simulation in a realistic time. For reference, a block diagram schematic is shown in Fig. 8.

The parameters shown in Eqs. (12) and (13) above are functions of the patient's height, weight, age, and sex. Typically, these parameters are set in the time scale of 1/min, however, since we are modeling delay and performing our control design in the timescale of seconds, the parameters were scaled to be in the correct time scale. The equations for nonstatic parameters are shown below with the variables m, h, and a representing mass in kilograms, height in meters, and age in years, respectively. The static parameters are as follows: k₁₃ is set to 3.3e−3, k₃₁ is set to 5.83e−5, and k_eo is set to 7.6e−3:

Before we send the effect-site concentration straight to the Wilson–Cowan equations, we must develop a relationship to convert this parameter to the nondimensional representation that the Wilson–Cowan equations use. Based on the relationship established in the BSR algorithm in Fig. 8, we can use a similar relationship described in Ref. [9] between effect-site concentration and burst suppression probability (BSP) to give an approximation for this transformation. A sigmoid model was made to approximate the BSP versus effect-site concentration shown in Ref. [9] using Eq. (24):

It is shown in Ref. [9] that this relationship holds for a certain patient parameterization and is also used for all later experiments run with various patients. The concentration to c₂ values derived from this method is shown in Fig. 9.

To choose the threshold value, we can take a linear fit of the sigmoid shown in Sec. 2.1 where we computed the nominal relationship between c₂ and BSR. Figure 10 shows the original BSR curve with the linear approximated plotted over it.

Prior to the plant model, we have a second-order actuator model and our controller. For this problem, a proportional-integral-derivative (PID) controller architecture was chosen. A specific goal of this study was to minimize overshoot. Implementing derivative control in addition to a classic proportional-integral (PI) controller is an excellent way to minimize overshoot in a system as well as minimize steady-state error. The equations which describe the controller and actuator are below.

Finally, two switches were added to disallow unrealistic results to come from the controller or the BSR output. First, a switch was put before the PK model that limited the input to the system to be positive (i.e., always providing a propofol infusion rate to the system). This was done because the only inputs we can put into the system are a positive one or zero. However, a useful quality of the PK model is that it is asymptotically stable about a zero-concentration equilibrium point. In this way, it models the system's ability to effectively decay its concentration if no additional propofol is fed into the system. Second, a max block was put at after the noise is added to the output of the BSR estimation algorithm. This was done so that if the BSR estimation algorithm is estimating a near-zero BSR, the noise will not allow a negative BSR to be estimated as that is not physiologically possible.

In this section, the PID control design and actuator specs will be discussed as well as final simulation results with the gains and actuator included.

For our PID controller and actuator, there were several goals that needed to be met in order to claim the gains and actuator parameters had been sufficiently designed to meet the system requirements. From a time-domain perspective, we desired minimal overshoot, a modest rise time, and sufficient robustness in terms of the noise added at the input and output of the plant, as well as parameter uncertainty. In terms of the time domain requirements, the rise time was a much lower priority for these systems than the steady-state error or the overshoot. This is due largely in part to the long periods of time that the system will need to be held at a certain level. Of course, it will have limits on how slow of a system response it can be, but it will not be a driving requirement. Overshoot and steady-state error on the other hand are much more important requirements for our application. Overshoot is important because of our control system's lack of ability to remove any propofol from the system. It cannot command a negative infusion rate. Therefore, if a target BSR is overshot, our closed-loop system must rely on the natural release of the propofol from the effect-site compartment to get down to the desired level. A large steady-state error is undesirable for obvious reasons that apply to any control system, but for our application, a large steady-state error could command a much higher propofol dose than the patient requires and thus result in patient overdose.

This equation tells us that the closed-loop response, T, must be balanced with the sensitivity of the system, S. It is desirable for our system to respond to low-frequency inputs (i.e., set-point values or step responses in BSR) and be robust to any high-frequency oscillations that could occur due to noise. In this way, the magnitude of the closed-loop system, T, should be 1 for low-frequency inputs and 0 at high-frequency inputs. The opposite should be true for S in order to be robust to noise. An easy tuning method to ensure this will be the case is to look at the loop gain throughout the frequency range. The loop gain should “roll off” at high frequencies which ensure our closed-loop system is not responding to high-frequency input. In order to maximize gain and phase margin, it is desirable to look at how quickly the loop gain of the system “rolls off” at the loop gain crossover frequency. Making this value near to −20 dB/decade ensures that we are getting sufficient gain and phase margin for our closed-loop system.

The thicker line is the step response while the thinner line is the tangent line at the steepest point in the step response. L is the time from 0 to the x-axis intersection with the tangent line and T is the time from the tangent line intersection with the x-axis and the tangent line intersection with the final value of the step response. Using this method a value of 98.66 s was determined for L and a value of 501.95 s was determined for T. This yielded a K_p of 6.1, a K_I of 0.31, and a K_d of 301.2. It should be noted that the step response for this system reaches steady-state at nearly 10,000 s which is why we do not see the step response reach its final value in this figure.

While using these coefficients for design, we also swept the actuator parameters to determine the overall closed-loop performance from a linear perspective. A range of damping coefficients from 0 to 1 was run as well as a logarithmic range of time constants.

At each iteration, a linear closed-loop system was formed to estimate the nonlinear closed-loop system that will ultimately be simulated. The PK model, being linear initially, was coupled with a Pade approximation for the delay to form the plant model. A Pade approximation effectively desires to fit a certain equation order to a delay in order to approximate it in a linear sense. In this way, we can linearly approximate the effects of algorithmic delay in the system. A third-order approximation was found to be sufficient for our application.

This was then used to form the sensitivity and co-sensitivity systems as shown in Eqs. (31)–(33). The co-sensitivity, T, also being the closed-loop response, was then used to get the time domain performance out of the system as well as gain and phase margin. The loop gain was also taken independently, and the loop gain crossover frequency was calculated. The slope of the loop gain (dB/decade) was then determined at this point. Of course, if the chosen actuator and controller destabilized the system then the iteration would be considered a failure, and the design would iterate to the next range of values. The stability of the system was simply determined by taking the eigenvalues of the closed-loop system and ensuring they were all negative, which indicates all the poles of the system are on the left half-plane and the system is stable. The various combinations of actuator time constant and damping coefficient yielded a variety of results. We desired to maintain a gain margin above 6 dB and a phase margin greater than 45 deg in addition to minimal overshoot and a roll-off steepness near −20 dB. To meet these design requirements and shape a co-sensitivity and sensitivity maximum value that was relatively small as well as minimize overshoot, a damping ratio of 1 and a time constant of 91.03 s was chosen.

These actuator parameters yielded excellent gain and phase margin results. With the Ziegler–Nichols tuning method shown in Eq. (32), the system yielded a gain margin of 26.4 dB and a phase margin of 160.45 deg as well as a roll-off steepness of −21 dB at the loop gain crossover frequency. These excellent stability margins should not be altogether surprising since our system was already quite stable.

In addition, Fig. 12 shows the step response and sensitivity/co-sensitivity functions. The figures show a steady-state overshoot with a 2% steady-state error. Both functions show magnitudes over 1; however, they are relatively small and since our system has an excessive margin, these magnitudes should not pose any stability issues.

In order to test the model's ability to control the BSR of a certain patient, a BSR trajectory was created that correlates with the values in Table 1.

This trajectory was chosen to stress a large range of BSR values while also allowing us to view the steady-state values at each of these conditions by commanding the system to hold at that BSR for an extended period of time. The first simulation done is for the nominal recovery parameters shown in Table 2 with no noise added and with the patient data that the controller was designed to.

The results of this simulation, in Fig. 13, show the BSR closely tracking the step responses as they are commanded separately throughout the trajectory. An obvious first takeaway is that the initial BSR comes out to a value of roughly 0.15. This is because the minimum burst suppression ratio our algorithm could estimate while maintaining a relatively monotonic increase throughout the range of c₂ is also 0.15. Due to this, our conversion from concentration in the effect-site compartment was limited to the minimum c₂ value that could be attained with near monotonic performance. These results show that there is a steady-state error after the system is commanded to a set point of 0.8 as well as a small initial overshoot. This is not entirely unexpected due to the step response results from our gain selection and actuator design in Sec. 2.2.

The second simulation, shown in Fig. 14, was performed with a different patient than what the gains were designed to. In order to model what would result in less propofol infusion than what the original patient needed in an attempt to cause the system to overshoot, a patient type was chosen with a female gender, age of 65, a weight of 50 kg, and a height of 125 cm. Due to the decreased weight and age and increased height in comparison to the original patient, this patient requires less propofol to reach a certain effect-site concentration.

The thinner line in Fig. 15 represents the nominal patient concentration while the thicker line shows the concentration of the patient with added uncertainty. The main difference shown between the two is that the system with patient uncertainty does pose a faster rise time and higher overshoot than the original patient. However, this overshoot is quickly dissipated, and a steady-state result like that of the nominal case is reached quickly. Aside from this, the closed-loop system performance shows the evident noise that was added to the system but maintains a mean value around the commanded BSR.

Finally, a simulation with a drastically different patient and the maximal increase in recovery parameters was run (shown in Fig. 16) to see how the system deals with maximal uncertainty in the system it was designed for. The maximal increase in recovery parameters is as shown in Fig. 5. The immediate impact of this increase is the relationship between BSR and c₂ which will in turn result in a change in the effect-site concentration required to achieve a certain BSR. This change also increased the effective lower bound of control due to the increased suppression that we witnessed in Fig. 5 when the recovery parameters were all simultaneously increased. This increased suppression, coupled with the static lower bound set on the concentration to c₂ conversion, caused the system to be unable to both attempt a lower c₂ value than 10 for the purpose of searching for a smaller BSR to be controlled and, therefore, unable to reach a BSR at the target level that we desired. The added noise to the system resulted in similar behavior to the model which we added patient uncertainty and noise. In addition, larger spikes than would be expected by noise generation occurred at high BSRs. It is likely that the algorithm had seen a c₂ value or range of c₂ values that caused this jump which could be unique to the recovery parameter uncertainty that was imparted on the model.

This paper presented a control systems design for detecting and controlling burst suppression in a low-dimensional biophysical model [12] in conjunction with a four-dimensional pharmacokinetics model. Our design employed robust control design strategies to account for the presence of patient and neurological uncertainty. Design goals of this system were sufficiently met in the nominal case and the case with added noise and patient uncertainty where the maximum observed overshoot was roughly 6% and the system's steady-state error was roughly 2%.

One potential caveat with regards to our paper is whether a detailed model is necessary to enable our control design, which is ultimately based on a classical PID backbone. While PID controllers do not overtly embed a model, they nonetheless typically require tuning to achieve satisfactory performance when deployed for a given model. In particular, for a model with several nontrivial nonlinearities such as ours, it is not obvious a priori that a generic tuning strategy will work. In this regard, a systematic analysis and design evaluation is both necessary and productive for most PID designs. Because testing these designs on actual patients is difficult, due to obvious safety concerns, having a plausible model upon which to assess performance issues is a key design step. This is, ultimately, the goal and contribution of our paper. Certainly, other robust control strategies may also be possible and could be worth exploring in the future, though our results suggest that a relatively interpretable PID strategy can be effective for this application.

There are no conflicts of interest.

Table 2 provides a description for the state variables used in the Wilson–Cowan model as described by Liu and Ching [12] and gives the initial conditions for each state variable for the simulations performed in this paper.

Table 3 gives the nominal parameterization of the Wilson–Cowan model as described by Liu and Ching [12].