Katie Hergenrother, Victoria Madison, and Elaine Smith

Final Project for CBE 30338, Spring 2017

Anesthesia allows a patient to enter into a state of unconsciousness to ensure that he or she does not experience extreme pain from a surgical procedure. An appropriate concentration of anesthesia must be administered for the patient to remain in this unconscious state for the entire procedure, but not longer. It is crucial that brain activity is not permanently affected by prolonged exposure to the drug. The focus of this project is how to properly dose patients with anesthesia for brain-related injuries.

In practice, anesthesiologists manually monitor brain activity and control the drug administration to the patient. However, a Closed-Loop Anesthetic Delivery (CLAD) system could be one solution to transforming a manual drug delivery process into a more systematic and controlled process. The CLAD system is modeled through monitoring brain activity patterns established by electroencephalogram (EEG) signals that are converted into burst suppression probabilities (BSP). Burst suppressions are characterized by alternating periods of suppressed brain activity with periods of electrical activity. By identifying a target burst suppression probability, a CLAD system computes how much anesthetic should be administered to a patient via an IV connected to an infusion pump. A proportional-integral controller alters the flow rate of the drug released from the infusion pump into a peripheral compartment in the body and into the brain. The level of anesthesia detected in the brain affects the brain activity, which produces different EEG signals. The frequency of bursts in brain activity detected by the EEG is then converted to a BSP and sent through a BSP filter algorithm to compute new probabilities that are used to adjust the controller until the desired amount of anesthetic has been properly administered to the patient to keep him or her unconscious for a medical procedure. Every patient has a different epidemiologic profile and therefore cannot be treated in the same way, so a control feedback loop is necessary to achieve an individualized treatment plan.

Burst suppression is a measure of electrical activity and inactivity in the brain, monitored by electroencephalogram (EEG) markers. The EEG markers are connected to a patient’s head externally near the brain. Electrical impulses from the brain interact with the electrodes in the EEG and create wave signals that can be translated into brain activity. A PI controller and a closed-loop system take an electrical impulse reading, analyze the wave pattern, and convert this signal into a control for the entrance of the anesthesia from the IV, through a compartmentalized location within the blood of the body and across the blood-brain barrier as necessary to either stimulate or suppress activity, depending upon the reading. This continues until the error in the brain activity level of the patient is at a minimal level compared to a standard set for the patient based on their needs.

A technological challenge arises when converting the electrical signal into a drug concentration flow rate. A randomized Gaussian probability distribution is incorporated into the control loop design to simulate electrical signals to convert between electrical and concentration readings. Another challenge that arises is that the drug cannot be safely administered directly to the brain to only affect this region of the body. Anesthesia must enter and travel through a centralized vein and to the brain as needed. The amount of anesthesia that enters the brain compared to the amount that remains compartmentalized in the blood must be continually monitored and maintained through feedback control. Overshoot should be limited in this system because if the drug concentration exceeds the target level of burst suppression by too significant of an amount, the patient could have negative effects such as a prolonged coma or possibly death.

The goal of the project will be to develop and implement a PI control model for the CLAD system for a patient with a brain injury. The model ensures that the appropriate amount of anesthesia is administered and manual adjustments are minimized. A diagram of the CLAD system is illustrated in Figure 1. The burst suppression ratio will be obtained through analysis of electrical brain signals from the EEG, flow rate of the anesthesia to the IV, and entrance of the drug from the compartmentalized location across the blood-brain barrier and to the desired area of the brain. Programming of the CLAD system is implemented to simulate the levels of anesthesia in a patient’s blood system.

*Figure 1. CLAD System Block Diagram*

In order to implement the CLAD system, an experiment is conducted over a time interval defined as T. Subintervals of time, k, refer to small divisions of the experiment time. The width of these subintervals is then defined as $\Delta=TK$ where K is a large integer used to divide the experimental time and allows a given time to be expressed as t=k. It is assumed that K is very large in order to ensure that each width is small and that only one event can occur during this width of time. The EEG produces patterns that describe whether or not an event occurs. An event is defined as either a burst or a suppression. Only one event can occur during a sampling of time so that a binary model can be used for which $n_k$ represents the binary event. The convention $n_k=1$ defines a burst and $n_k=0$ defines a suppression.

A target BSP value, $p_{target}$, is determined based on the size and gender of the patient to which anesthesia is being administered. This target value is converted into an amount of anesthetic that is targeted to administer to the patient’s brain, $x_2^{target}$. The target amount of anesthesia to be administered is compared to the amount of anesthesia computed from the CLAD system controller, $x_2$, and the error between the two values is minimized. The CLAD system controller takes continuous readings of the EEG signals and computes a probability of bursts in brain activity based on the frequency of these events. A BSP during a given subinterval is defined as

$$p_k=\frac{1-\exp(-x_2)}{1+\exp(-x_2)} [1]$$

which is a sigmoidal function for which $p_k$ is the probability of a suppression in a given subinterval, k, and $x_2$ is the amount of anesthesia present in the patient’s brain in the subinterval. The BSP and drug amount are inverse functions of eachother, so the amount of propofol present in the brain can be represented by

$$x_2=\ln(\frac{1+p_k}{1-p_k})[2]$$based on the probability of a suppression in a subinterval. Instead of comparing BSP values of the target and control system, amounts of anesthesia are compared. An amount of the drug can be more easily translated into release of anesthesia from the infusion pump into the patient’s body.

A two compartment model, illustrated in Figure 2, is used to simplify the system. In other applications, three or four compartmental models are used. The model is simplified based on the assumption that burst suppression has a very limited dynamic range. This means that once burst suppression is achieved, the only transitions that occur are varying degrees of burst suppression. If the flow rate of the drug is increased the burst suppression will increase until eventually an isoelectric EEG occurs.

*Figure 2. Two Compartment Model Diagram*

If the flow rate of the drug is decreased, burst suppression decreases. The change in amount of anesthesia in each compartment can be expressed by

$$\frac{dx_1}{dt}=u(t)+ x_2a_{21}-x_1(a_{10}+a_{12})[3]$$$$\frac{dx_2}{dt}=x_1a_{12}-x_2a_{21}[4]$$

For which the rate parameters, $a_{12}$, and $a_{21}$ indicate the flow from the central compartment into the brain compartment and the flow from the brain to the rest of the body, respectively, while $a_{10}$ determines the rate that anesthesia leaves the body. The only way for anesthesia to leave the body is through the central compartment.

The CLAD system controls the amount of anesthesia entering the brain to achieve the desired BSP. The controller uses proportional integral control where the infusion rate, $u(t)$, is

$$u(t)=a_pE(\tau)+a_i\int_0^\tau\mathrm E(\tau)\,\mathrm{d}\tau[5]$$

using $a_p$ as the coefficient for proportional control and $a_i$ as the coefficient for integral control. The error, $E(t)$, is given by

$$E(t)=x_2^{target}-x_2(t)[6]$$where $x_2(t)$ is the amount of anesthetic in the brain at a given time, and $x_2^{target}$ is the target value for the concentration. This target value relates to a target BSP, such that

$$x_2^{target}=f^{-1}(p^{target})[7]$$

from the relationship between Equations 1 and 2.

The values for the rate parameters and coefficients are individualized for the patient to ensure that the system is critically damped because overshoot is harmful in anesthetic delivery. A transfer function to convert an error reading in the amount of anesthesia administered to the brain into an infusion pump flow rate can be determined from a second order differential equation. For modeling purposes, this second order equation is written as a system of two first order differential equations given by

$$dx=x_2$$$$dx_2=-(a_{10}a_{21})x-(a_{10}a_{21})x_2+(ba_{12})u$$where $b$ is the gain coefficient for the system and the rate parameters are chosen rate parameters that are adjusted for the patient's needs, determined by their age, gender, weight, and height.

Spikes present on an EEG represent the electrical impulse bursts of brain activity. As the BSP increases, there is an observed decrease in the number of spikes meaning that the patient’s brain is increasing in inactivity over time. This trend is as expected because for a greater dosage of propofol, the patient is desired to be in a deeper state of sedation, and therefore the brain activity should be less than if the person was not under an anesthetic.

The graph of the amount of propofol in the brain, $x_2$, versus time without the EEG connected to the control loop exhibits the expected shape. The concentration of propofol is initially set to zero and rapidly increases and oscillates around the target value of 1.7 mcg. There is overshoot which could not be reduced because the coefficient for integral control, $a_i$, was severely sensitive. There is minimal ringing and the graph settles and stays at the target amount of propofol in the brain. Adjusting the parameters $a_i$ and $a_p$ from the literature values gave rise to a better model that reaches the target value in less time with less ringing.

The graph of $x_2$ versus time with the EEG connected to the loop has a significant amount of ringing present because this section introduced the randomized Gaussian distribution based on the fact that the brain’s response to anesthesia cannot be accurately predicted. Again the graph has the expected behavior of increasing quickly from zero to the target amount of propofol in the brain. The adjusted parameter values of $a_i$ and $a_p$ improve the dosage of propofol from the suggested literature values.

A maximum flow rate of 200 mcg/sec is imposed to ensure the safety of the patient. The flow rate swiftly increases from zero to 200 mcg/sec and the flow rate stays saturated at 200 mcg/sec for about 100 seconds then quickly drops down to no flow rate and then rings with slower changes in flow rate around the value of 60 mcg/sec. Again the improved parameters of $a_i$ and $a_p$ give a better model because the flow rate reaches the slower ringing in less time. The graph for flow rate has less ringing than the graph for amount of propofol with the EEG in the loop because the randomized contribution is weighted less than the initial condition value so that quick,dramatic over-corrections are not present.

The burst suppression probability versus time graph has the desired behavior of increasing from zero to the target burst suppression probability of 0.7 with minimized ringing. As can be seen on the graph, the proposed parameters for $a_i$ and $a_p$ improve the response of BSP compared to the literature values.

The literature values for proportional and integral control resulted in higher overshoot than desired and longer time to dose to the target range for propofol in the brain. Different control values were proposed in the hopes of reducing the overshoot and the time to reach the target range. The amount of overshoot was affected by the amount of integral control $a_i$, making the overshoot difficult to reduce because $a_i$ was very sensitive. When $a_i$ was greater than 0.98 1/s, the target value of $x_2$ was reached with overshoot. When $a_i$ was less than 0.98 1/s, $x_2$ never reached the target value but was more damped. A proposed value of $a_i$ was set to 1 1/s, and increasing $a_i$ further had no discernable effect on control.

It is desired to reach the respective target values quickly so that the patient reaches the appropriate level of brain inactivation and unconsciousness for the procedure to begin. The amount of time it took to reach the target value was decreased by increasing proportional control, $a_p$. However, a high value of $a_p$ also produced more ringing. A proposed value of $a_p$ was set to an optimal value 183.3 1/s. Increasing $a_p$ further did not significantly reduce the time to reach the target value, and a value of 183.1 1/s had low enough ringing. These trends were most clearly seen in the Infusion Pump Compartment Model without EEG, but also applied to the Full CLAD System model.

The brain responds to propofol in a systematic, yet random manner. As the amount of propofol in the brain increases, the activity of the brain decreases and vice versa. Although this trend is always observed, there is not a consistent electrical impulse burst and suppression pattern based on a given amount of propofol. The CLAD model introduces Gaussian distributed randomized bursts in the EEG portion of the control loop to account for this unpredictability. The ringing found in the plots of amount of propofol in the brain, pump infusion rate of propofol, and burst suppression probability versus time plots illustrates the randomization. With each run of the model, the amount of the ringing differs slightly. Since the variation in plots due to randomization and ringing made interpreting and tuning the controller difficult, the model that is only affected by the infusion pump and two-compartment body and brain system was used to determine the appropriate proportional and integral control constants. It was assumed that introducing the EEG would not alter the overall shape or behavior of the system, so the parameters could be kept consistent between the simplified and more complex simulations.

The ringing in the pump infusion rate of propofol was reduced through the PI controller. By taking the difference between the error in the current time step of the control loop and the error in the previous time through the loop and using this new error to improve proportional control, the error was reduced and the flow rate fluctuated less than it would with a larger, more unpredictable error value. Ringing was also reduced by determining an optimal amount of pump infusion flow that was affected by the error reading. The overall infusion flow rate was split into one portion that was affected by the error and one portion that was affected by the initial condition for the infusion flow rate at a given time. By reducing the contribution of the term affected by the error, which was dependent upon the randomized EEG reading, the ringing was reduced. The ringing should be kept to a minimum because the fluctuation in the flow rates is harmful to the patient and anesthesiologists would prefer more consistent flows that are adjusted slightly, not as drastically as if the flow were not balanced by the initial condition contribution.

The patient is a 40 year old, 165 cm, 60 kg woman. The model rate parameters can be adapted for additional patient ages, genders, heights, and weights.

The following Python libraries are used for computations of the CLAD system model.

In [1]:

```
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
from scipy.integrate import odeint
from scipy.integrate import quad
import random
```

The parameters defined in this CLAD system model are found through Liberman's scholarly paper "A Closed-Loop Anesthetic Delivery System for Real-Time Control of Burst Suppression". The rate parameters determine the rate by which propofol enters and exits the two compartments, the body and brain. The gain coefficient, b, is used to build a PI Controller that monitors the administration of propofol from an infusion pump.

A BSP probability of 0.7 was chosen to be the target value as determined through literature research. This value was maintained most frequently during an anesthetic procedure. Since the BSP probability is inversely related to the concentration of propofol in the brain, a target value of approximately 1.74 mcg propofol was established for the setpoint of the control loop.

Experimental parameters exist for a proof-of-concept 10 minute trial. EEG readings from the brain are taken every second. It is assumed that the time of the trial starts when the patient is wheeled into the operating room and therefore does not have any initial amount of propofol in her brain or peripheral compartment.

In [2]:

```
#Patient 4 rate parameters, scale factor, and control parameters
a12 = 2.26e-4/60 #sec-1
a21 = 0.731/60 #sec-1
a10 = 0.732/60 #sec-1
b = 0.1177
#Target conditions
p_target = .7
x2_target = np.log((1+p_target)/(1-p_target))
#Experimental parameters
T = 600.0 #sec
K = 120000.0
delta = T/K #sec
k = np.linspace(1,K,K)
t = k*delta #sec
dt = 1 #sec
tgrid = np.arange(0, 600, dt)
# Initial Conditions
x2 = 0 # mcg must be in the range (0,5.5)
x1 = 0 # mcg must be in the range (0,5.5)
u = 0 #mcg/sec
```

A second order differential equation, rewritten as a system of two first order equations, models a transfer function from the error to the infusion rate of propofol.

In [3]:

```
E = x2_target - x2
IC = [0, 0, 0] #assume steady state model initially
def deriv(Z,t):
X,X2,In = Z
dX = X2
dX2 = -(a10*a21)*X-(a10*a21)*X2+(b*a12)*u
dIn = E
return[dX,dX2,dIn]
sol = odeint(deriv,IC,t)
X,X2,In = sol.transpose()
```

For initial testing, only the mechanics of the infusion pump were modeled. This model accounted for the infusion of propofol from the pump and the concentration administered to the peripheral body compartment and brain. This simplified model focused only on the target concentration desired in the brain, but not the brain's activity monitored by electrical responses from an EEG. The error in the concentration of propofol in the brain computed after the compartmental model was compared to the target value.

The simplified model was run twice. The first run utilizes literature obtained proportional and integral control parameters. The second run was adapted to reduce the dosing time of propofol to the brain. New values for the proportional and integral control constants were determined through trial and error and are proposed for improving upon the previously existing model. The results from both PI controllers are compared and displayed in the following plots.

In [4]:

```
# LITERATURE VALUES FOR CONTROLLER
#PI Controller
ap = 8028/60 #sec-1 Proportional Control Constant
ai = 90/60 #sec-1 Integral Control Constant
IC2 = [0, 0, 0] # Initial Conditions, assuming that the patient begins the procedure without any anesthetic present in the body
def deriv2(W, t):
x1k, x2k, In = W
E = x2_target - x2k #Compute error between target and measured brain propofol concentrations
u = ap*E +ai*In #Adjust the PI controller for a new infusion rate of propofol from the infusion pump into an IV
dx1 = u - a12*x1k + a21*x2k - a10*x1k #Change in concentration of propofol in the peripheral compartment
dx2 = a12*x1k - a21*x2k #Change in concentration of propofol in the brain
dIn = E #Integrate the error to update the controller in time
return[dx1, dx2, dIn]
sol2 = odeint(deriv2,IC2,t) #Solve the system of differential equations
x1k, x2k, In = sol2.transpose()
o = np.ones(len(t))
x2_target_vec = x2_target*o
plt.plot(t, x2k) #Plot the Literature Values compared to the target value
plt.plot(t, x2_target_vec)
plt.ylim([0,4])
plt.ylabel('x2 [Amount of Propofol]')
plt.xlabel('Time (sec)')
# ADJUSTED VALUES FOR CONTROLLER
#Target parameters
p_target = 0.7
#Patient 4 rate parameters, scale factor, and control parameters
a12 = 2.26e-4/60 #sec-1
a21 = 0.731/60 #sec-1
a10 = 0.732/60 #sec-1
b = 0.1177
#Experimental parameters
T = 600.0 #sec
K = 120000.0
delta = T/K #sec
k = np.linspace(1,K,K)
t = k*delta #sec
# Initial Conditions
x2 = 0 # must be in the range (0,5.5)
x1 = 0 # must be in the range (0,5.5)
u = 0 # mcg/sec
#PI Controller
ap_adjust = 11000/60
ai_adjust = 60/60
IC2 = [0, 0, 0]
def deriv2(W, t):
x1k, x2k, In = W
E = x2_target - x2k
u_adjust = ap_adjust*E +ai_adjust*In
dx1 = u_adjust - a12*x1k + a21*x2k - a10*x1k
dx2 = a12*x1k - a21*x2k
dIn = E
return[dx1, dx2, dIn]
sol2 = odeint(deriv2,IC2,t)
x1k, x2k, In = sol2.transpose()
o = np.ones(len(t))
x2_target_vec = x2_target*o
plt.plot(t, x2k) # Plot the Literature Values
plt.ylim([0,3])
plt.ylabel('x2 [Amount of Propofol]')
plt.xlabel('Time (sec)')
plt.title('Propofol Control')
plt.legend(['Literature Control Values', 'Target Value', 'Proposed Control Values'])
```

Out[4]:

In order to build a more sophisticated model, the simple model above was adapted to include how electrical signals of brain activity are affected by the amount of anesthesia in the brain. The following portion incorporates a PI controller that is adjusted based on the error between a target amount, determined from a burst suppression probability of interest, and an amount of anesthesia affecting brain activity.

A Burst Suppression Probability (BSP), as explained in the problem statement, measures the frequency of electrical bursts of activity in the brain. A BSP value of 0.7 was determined as the target from literature. This model takes the amount of propofol in the brain, computed from the infusion pump mechanics, and observes how the brain randomly responds in bursts of activity based on the amount of propofol detected. A Gaussian distribution is used to randomize bursts of activity and compute probabilities of brain suppression. The probability measured in the brain is inversely proportional to a measured brain amount, which drives the PI controller.

In [5]:

```
# Proposed Value Controller
#PI Controller
ap_adjust = 12000/60 #sec-1 #Experimentally determined
ai_adjust = 60/60 #sec-1 #Experimentally determined
#Initial Conditions
IC2 = [0, 0] #Initial concentrations of propfol in the body and brain, assuming the procedure has not yet begun.
u_adjust_p = 0
u_adjust = 0
I = 0
log = [] #Variable to log all of the values computed through an iterative control loop approach
def deriv2(W, t):
x1k, x2k = W
dx1 = u_adjust - a12*x1k + a21*x2k - a10*x1k #Change in propofol concentration in the peripheral compartment
dx2 = a12*x1k - a21*x2k #Change in propfol concentration in the brain
return[dx1, dx2]
x1, x2 = IC2
Elast = 0
for t in tgrid:
sol2 = odeint(deriv2,[x1,x2],[t, t+dt])[-1] #Solve the set of differential equations over a given interval
x1, x2 = sol2.transpose()
lam = 200*(1-np.exp(-x2))/(1+np.exp(-x2)) #Probability of a burst
n = random.gauss(lam, np.sqrt(lam)) #Randomize the bursts of activity based on a Gaussian Distribution
p_p = n/200 #Determine the probability of bursts from the number of events observed
x2_EEG_p = np.log((1+p_p)/(1-p_p)) #Find the corresponding concentration of propofol in the brain
E = x2_target - x2_EEG_p #Compute the error between the target and measured value
u_adjust_p += ap_adjust*(E-Elast) +ai_adjust*E #Adjust the infusion pump, accounting for the current error in the concentrations of propofol in the brain
u_adjust = 0.07*u_adjust_p + (.93)*u_adjust #Reduce the amount of noise in the infusion pump
u_adjust = max(0,min(200,u_adjust)) #Ensure that the flow rate does not exceed a maximum value of 200
log.append([x2_EEG_p,t,u_adjust,x1,x2,p_p,E]) #Track all values for a given time
Elast = E #Save the previous error
x2_EEG_p,t,u_adjust,x1,x2,p_p,E = np.asarray(log).T
# Literature Value Controller
#PI Controller
ap = 8028/60 #sec-1
ai = 95.90/60 #sec-1
#Initial Conditions
IC2 = [0, 0]
u_p = 0
u = 0
I = 0
log = []
def deriv2(W, t):
x1k, x2k = W
dx1 = u - a12*x1k + a21*x2k - a10*x1k
dx2 = a12*x1k - a21*x2k
return[dx1, dx2]
x1, x2 = IC2
Elast = 0
for t in tgrid:
sol2 = odeint(deriv2,[x1,x2],[t, t+dt])[-1]
x1, x2 = sol2.transpose()
lam = 200*(1-np.exp(-x2))/(1+np.exp(-x2))
n = random.gauss(lam, np.sqrt(lam))
p_lit = n/200
x2_EEG_lit = np.log((1+p_lit)/(1-p_lit))
E = x2_target - x2_EEG_lit
u_p += ap*(E-Elast) + ai*E
u = 0.12*u_p + (.88)*u
u = max(0,min(200,u))
log.append([x2_EEG_lit,t,u,x1,x2,p_lit,E])
Elast = E
x2_EEG_lit,t,u,x1,x2,p_lit,E = np.asarray(log).T
o = np.ones(len(t))
x2_target_vec = x2_target*o
plt.plot(t,x2_EEG_p) #Plot the Proposed Control Values
plt.plot(t, x2_EEG_lit) #Plot the Literature Control Values
plt.plot(t, x2_target_vec) #Compare Values to the Target
plt.ylim([0,5])
plt.ylabel('x2 [Amount of Propofol]')
plt.xlabel('Time (sec)')
plt.title('Propofol Control')
plt.legend(['Proposed Control Values','Literature Control Values','Target Value'])
```

Out[5]:

The rate of propofol administered from the infusion pump is monitored over time. The change in the infusion rate is controlled by the PI Controller and the adjustments made by the controller are shown in the plot below vs time.

In [6]:

```
# Literature Control Values
plt.plot(t,u)
plt.ylabel('u [Infusion Rate of Propofol]')
plt.xlabel('Time (sec)')
plt.title('Infusion of Propofol')
# Proposed Control Values
plt.plot(t,u_adjust)
plt.ylabel('u [mcg/sec]')
plt.xlabel('Time (sec)')
plt.title('Infusion of Propofol')
plt.ylim([0,275])
#Maximum Infusion Rate
Max = 200
Max_vec = np.ones(len(t))*Max
plt.plot(t, Max_vec, '--')
plt.legend(['Literature Control Values', 'Proposed Control Values', 'Maximum Rate'])
```

Out[6]:

By controlling in the rate that propofol is administered to a patient, the concentration of propofol in the brain at a given time is changed. The brain's response to the concentration of the drug can be converted into activity "bursts". The target value of 0.7 was achieved through use of the PI Controller, as shown by the plot.

In [7]:

```
p_target_vec = np.ones(len(t))*p_target #Target value
plt.plot(t,p_p) #Plot proposed BSP values calculated in control loop above
plt.plot(t,p_lit) #Plot literature BSP values calculated in control loop above
plt.plot(t,p_target_vec) #Plot target BSP value
plt.ylabel('BSP')
plt.xlabel('Time (sec)')
plt.ylim([0,1.3])
plt.title('Burst Suppression Probability')
plt.legend(['Proposed BSP', 'Literature BSP','Target BSP'])
```

Out[7]:

An EEG creates a pattern of alternating bursts and suppressions while monitoring the brain's activity level. A patient who is dosed strongly with a drug will experience fewer bursts than a patient who is not under an anesthetic. The following simulations show an example of a randomized brain activity response to propofol. The first EEG is the result of a 0.1 target BSP, the second is a result of a 0.7 target BSP, and the final simulation is the result of a 0.9 target BSP.

In [8]:

```
pk_EEG = .1
n = []
j= 0
for i in k:
if random.random() >= (1-pk_EEG)*delta: #have to have inversefunc to find pk and make this loop run
nk = 0 #suppression
n.append(nk)
j = j + 1
else:
nk = 1 #burst
n.append(nk)
t_EEG =np.linspace(1,len(n),len(n))
plt.plot(t_EEG[0:20000],n[0:20000], color = 'black')
plt.xlabel('Time (sec)')
plt.ylabel('EEG Burst')
plt.title('EEG (p=0.1)')
```

Out[8]:

In [9]:

```
pk_EEG = .7
n = []
j= 0
for i in k:
if random.random() >= (1-pk_EEG)*delta: #have to have inversefunc to find pk and make this loop run
nk = 0 #suppression
n.append(nk)
j = j + 1
else:
nk = 1 #burst
n.append(nk)
t_EEG =np.linspace(1,len(n),len(n))
plt.plot(t_EEG[0:20000],n[0:20000], color = 'black')
plt.xlabel('Time (sec)')
plt.ylabel('EEG Burst')
plt.title('EEG (p=0.7)')
```

Out[9]:

In [10]:

```
pk_EEG = .9
n = []
j= 0
for i in k:
if random.random() >= (1-pk_EEG)*delta: #have to have inversefunc to find pk and make this loop run
nk = 0 #suppression
n.append(nk)
j = j + 1
else:
nk = 1 #burst
n.append(nk)
t_EEG =np.linspace(1,len(n),len(n))
plt.plot(t_EEG[0:20000],n[0:20000], color = 'black')
plt.xlabel('Time (sec)')
plt.ylabel('EEG Burst')
plt.title('EEG (p=0.9)')
```

Out[10]:

The target values for the amount of propofol in the brain and the burst suppression probability were achieved, and the infusion pump flow rate stayed within the safe range of 0 to 200 mcg/sec. Values for proportional and integral control parameters that differ from literature values were proposed to reduce the time to reach the target values for the appropriate BSP in the brain, decrease ringing, and keep overshoot low. As expected, an increase in amount of propofol in the brain corresponds to a higher BSP and less brain activity, as shown by the EEG pattern. Since the electrical signals cannot be replicated to an exact pattern, randomization is introduced to the model to account for the unpredictability. The ringing in the values was reduced by altering the contribution of the infusion flow rate from the error in the setpoint and the measured value at each time interval.

* Literature Cited: *

[1] Brown E, Chemali J,Ching S, Liberman M. A Closed-Loop Anesthetic Delivery System for Real-Time Control of Burst Suppression. J Neural Eng. 2013 August

[2] "Propofol Dosage Guide with Precautions." Drugs.com, N.p., n.d. Web. 30 April 2017.

[3] "DIPRIVAN Infusion Rate Calculator for Sedation." About DIPRIVAN. N.p., n.d. Web. 30 April 2017.