## Abstract

Burst suppression is a phenomenon in which the electroencephalogram (EEG) of a pharmacologically sedated patient alternates between higher frequency and amplitude bursting to lower frequency and amplitude suppression. The level of sedation can be quantified by the burst suppression ratio (BSR), which is defined as the amount of time that an EEG is suppressed over the amount of time measured. Maintaining a stable BSR in patients is an important clinical problem, which has led to an interest in the development of BSR-based closed-loop pharmacological control systems. Methods to address the problem typically involve pharmacokinetic (PK) modeling that describes the dynamics of drug infusion in the body as well as signal processing methods for estimating burst suppression from EEG measurements. In this regard, simulations, physiological modeling, and control design can play a key role in producing a solution. This paper seeks to add to prior work by incorporating a Schnider PK model with the Wilson–Cowan neural mass model to use as a basis for robust control design of biophysical burst suppression dynamics. We add to this framework actuator modeling, real-time burst suppression estimation, and uncertainty modeling in both the patient's physical characteristics and neurological phenomena to form a closed-loop system. Using the Ziegler–Nichols tuning method for proportional-integral-derivative (PID) control, we were able to control this system at multiple BSR set points over a simulation time of 18 h in both nominal, patient varied with noise added and with reduced performance due to BSR bounding when including patient variation and noise as well as neurological uncertainty.

## 1 Introduction

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:

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

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

## 2 Burst Suppression Ratio Modeling and Simulation

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.

### 2.1 Modeling Burst Suppression Due to Neuronal and Metabolic Dynamics.

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:

*j*represents the column of neuron tissue that is being described. Additionally, the subscript

*k*describes a column that would be coupled with column

*j*. This specific study seeks only to analyze a single column of tissue and thus the equations will be hereafter described without the additional subscript or summation for the interconnectivity between columns. Thus, the de-coupled equations for describing neuron activity are as follows:

*a*and

*θ*represent tuning variables that are typically used in such functions to alter the slope and midpoint of the sigmoid. The modulating input mentioned above is described by the following equation:

*μ*

_{1}is the time constant for the autonomous part of the dynamics and

*μ*

_{2}is the time constants of the sigmoidal part of the dynamics. $k\varphi $ is a tuning variable for the sigmoid slope shape,

*η*represents the midpoint of the sigmoid when the input

*M*is zero. The input

*M*acts as a gating variable for

*ϕ*. The input

*M*is essential in describing the key features of this model as it attempts to accurately model the underlying physiological events occurring. It is a function of two variables, a consumption variable,

*g*

_{c}, and a recovery variable,

*g*

_{r}

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*_{2} 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*_{2} is varied. In Ref. [12], the bounds were set by the system's stability about the *c*_{2} 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*_{2} 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*_{2} increases from 7.7 to 69. Figures 3 and 4 give an example for how the model behaves with multiple *c*_{2} values.

The values of *c*_{2} 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.

### 2.2 Modeling Physiological Variability.

In addition to simply running a propofol (*c*_{2}) 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\beta ,\zeta $, and *ν* to simulate how different patients may respond to propofol as it enters their system. Figure 5 shows the same *c*_{2} 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*_{2}. 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\beta $ 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.

### 2.3 Characterizing Burst Suppression Ratio Versus Modeled Inhibition.

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.

*m*

_{s}for the suppressed mean magnitude and

*m*

_{b}for the bursting mean magnitude). Finally, the gain

*k*

_{b}was put on

*m*

_{b}to place the bursting and suppressed magnitudes in a similar value range to make both effective inputs. The sigmoid equation below then gives an output from 0 to 1 to describe the suppression at that time-step

A positive value for the exponent was chosen because with lower ranges of *c*_{2}, 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*_{2}, *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*_{2}. Figure 6 shows the BSR calculated as a function of *c*_{2}. This algorithm calculated a higher value for the lower ranges of *c*_{2}, specifically between 8 and 10. After a *c*_{2} value of 10, the BSR increases in a near monotonic fashion until it reaches a near 1 BSR at a *c*_{2} value of 69.

Despite these limitations on the lower ranges of *c*_{2}, 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*_{2} where the BSR stops decreasing and begins its monotonic increase to 1.

### 2.4 Burst Suppression Ratio Moving Average.

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*_{2} but very quick convergences for most values of *c*_{2}. This is likely due to the complete absence of higher frequency bursting data at high *c*_{2} 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.

## 3 Control System Architecture and Design

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.

### 3.1 Pharmacokinetic Model.

*x*

_{1}, describes the concentration of the drug in the central compartment, or the blood, of the patient. Compartments

*x*

_{2}and

*x*

_{3}describe the fast and slow compartments, respectively. This can be easily seen by the magnitude of the coefficients that describe the second-order relationship between these compartments and the central compartment. Finally, there is

*x*

_{4}, or the effect-site compartment concentration. This represents the concentration of propofol in the brain, which we can then use as an estimate for the variable,

*c*

_{2}, which is an input to the Wilson–Cowan equations

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*_{13} is set to 3.3e−3, *k*_{31} is set to 5.83e−5, and *k*_{eo} is set to 7.6e−3:

### 3.2 Effect-Site Concentration to Modeled Inhibition.

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):

### 3.3 Control Resolution.

*c*

_{2}values were then sent to the Wilson–Cowan equations for the dynamics to be calculated and fed through to the BSR estimation algorithm. The BSR estimation algorithm here is implemented in the same way that is shown in Eqs. (11)–(15), with the exception that a reset line has been incorporated for when a new

*c*

_{2}value is being fed through to the Wilson–Cowan equations. This is necessary since we need to achieve a certain confidence level for each measured BSR. Hence, when a new value of

*c*

_{2}is fed through the system, the mean and signal variance blocks must be reset. In order to determine the

*c*

_{2}increments at which the delay is reset, we must look at the performance of the control system. If the system naturally increases

*c*

_{2}so quickly that the delay caused by the algorithm cannot complete before a new

*c*

_{2}value is fed through, the system will be caught in a state of constant delay. In this state of constant delay, the BSR will never claim convergence at the current propofol concentration before the concentration changes so drastically that the estimated BSR no longer accurately reflects the state of the system. Hence, a control resolution can be found for BSR as a function of the

*c*

_{2}delay reset value. To coincide with the confidence level we set of 0.0001 BSR variance, we set criteria for each new

*c*

_{2}value that is sent into the system as described in Eq. (26):

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*_{2} and BSR. Figure 10 shows the original BSR curve with the linear approximated plotted over it.

*c*

_{2}, there is an approximate increase of 0.0143 BSR. Thus, we can use a similar relationship that the BSR uses for convergence to calculate what the two-norm of the difference between the next

*c*

_{2}value and the previous

*c*

_{2}should be to reset the delay

### 3.4 Controller and Actuator Model.

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.

### 3.5 Constraints.

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.

## 4 Control Design, Results, and Discussion

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.

### 4.1 Control Law and Actuator Specifications.

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.

*T*and

*L*may be described using the step response of the plant. As an example, the plant model with both the PK model and the selected delay is shown in Fig. 11.

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.

### 4.2 Simulation Test Case Results and Discussion.

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.

BSR | 0.8 | 0.8 | 0.5 | 0.5 | 0.2 | 0.2 |

Time (s) | 0 | 25,000 | 25,001 | 45,000 | 45,001 | 65,000 |

BSR | 0.8 | 0.8 | 0.5 | 0.5 | 0.2 | 0.2 |

Time (s) | 0 | 25,000 | 25,001 | 45,000 | 45,001 | 65,000 |

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*_{2} is also 0.15. Due to this, our conversion from concentration in the effect-site compartment was limited to the minimum *c*_{2} 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*_{2} 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*_{2} conversion, caused the system to be unable to both attempt a lower *c*_{2} 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*_{2} value or range of *c*_{2} values that caused this jump which could be unique to the recovery parameter uncertainty that was imparted on the model.

## 5 Conclusion

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.

## Conflict of Interest

There are no conflicts of interest.

### Appendix: Wilson–Cowan Model State Description and Parameterization

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.

Variable | Description | Units | Initial condition |
---|---|---|---|

e | Excitatory firing rate | Substrate/s | 0.3 |

i | Inhibitory firing rate | Substrate/s | 0.1 |

ϕ | Slow process firing rate | Substrate/s | 0.25 |

M | Consumption/recovery gating variable | Substrate/s | 10 |

β | Recovery evolution | Substrate/s | 0.3 |

Variable | Description | Units | Initial condition |
---|---|---|---|

e | Excitatory firing rate | Substrate/s | 0.3 |

i | Inhibitory firing rate | Substrate/s | 0.1 |

ϕ | Slow process firing rate | Substrate/s | 0.25 |

M | Consumption/recovery gating variable | Substrate/s | 10 |

β | Recovery evolution | Substrate/s | 0.3 |

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

Parameter | Description | Units | Value |
---|---|---|---|

k_{e}, k_{i} | Maximal value of the excitatory and inhibitory response functions | Substrate/s | 1, 1 |

r_{e}, r_{i} | The absolute refractory period of the excitatory, inhibitory subpopulation | Nondimensional | 1, 1 |

ω_{e}, ω_{i} | Wilson–Cowan time constants | 1/s | 500, 500 |

P, Q | Level of background excitation in the excitatory, inhibitory subpopulation | Substrate/s | 900, 0 |

c_{1}, c_{3} | Average number of excitatory synapses per cell | Nondimensional | 16, 15 |

c_{2}, c_{4} | Average number of inhibitory synapses per cell | Nondimensional | 12, 3 |

θ_{e}, θ_{i}, a_{e}, a_{i} | Maximal slope parameters of the logistic curve for the excitatory, inhibitory subpopulation | Nondimensional | 4, 3.7, 1.3, 2 |

μ_{1}, μ_{2} | Modulation time constants | 1/s | 2, 2 |

$k\varphi $ | Sensitivity to the variations of the metabolic substrate | 1/substrate | 75 |

$k\beta $ | Sensitivity to the variations of the neuronal activity | 1/substrate | 4.5 |

η, ζ | Average number of inhibitory synapses per cell | Substrate | 0.25, 0.27 |

k_{r}, k_{c} | Metabolic recovery and consumption rates | Substrate/s | 900, 700 |

ν_{1}, ν_{2} | Homeostatic autoregulation time constants | 1/s | 0.08, 0.08 |

Parameter | Description | Units | Value |
---|---|---|---|

k_{e}, k_{i} | Maximal value of the excitatory and inhibitory response functions | Substrate/s | 1, 1 |

r_{e}, r_{i} | The absolute refractory period of the excitatory, inhibitory subpopulation | Nondimensional | 1, 1 |

ω_{e}, ω_{i} | Wilson–Cowan time constants | 1/s | 500, 500 |

P, Q | Level of background excitation in the excitatory, inhibitory subpopulation | Substrate/s | 900, 0 |

c_{1}, c_{3} | Average number of excitatory synapses per cell | Nondimensional | 16, 15 |

c_{2}, c_{4} | Average number of inhibitory synapses per cell | Nondimensional | 12, 3 |

θ_{e}, θ_{i}, a_{e}, a_{i} | Maximal slope parameters of the logistic curve for the excitatory, inhibitory subpopulation | Nondimensional | 4, 3.7, 1.3, 2 |

μ_{1}, μ_{2} | Modulation time constants | 1/s | 2, 2 |

$k\varphi $ | Sensitivity to the variations of the metabolic substrate | 1/substrate | 75 |

$k\beta $ | Sensitivity to the variations of the neuronal activity | 1/substrate | 4.5 |

η, ζ | Average number of inhibitory synapses per cell | Substrate | 0.25, 0.27 |

k_{r}, k_{c} | Metabolic recovery and consumption rates | Substrate/s | 900, 700 |

ν_{1}, ν_{2} | Homeostatic autoregulation time constants | 1/s | 0.08, 0.08 |