5. PID control of atracurium

This example reproduces the example presented in Figure 5 of the paper Merigo et al., 2018. All the code detail to perform the simulation is available on the GitHub project here.

# Third party imports
import numpy as np
import pandas as pd
from bokeh.models import HoverTool
from bokeh.layouts import row, column
from bokeh.plotting import figure, show
from bokeh.io import output_notebook
from bokeh.models import Span

# Local imports
from pid_sim_functions import Patient_table_atracurium, PID_NMB
import python_anesthesia_simulator as pas

output_notebook()

Loading BokehJS ...

5.1. Perform table population simulation

Test the PID on the 12 patients of the tuning dataset.

p1 = figure(width=600, height=300)
p2 = figure(width=600, height=300)

# PID tuning parameters
Kp = 64.9721/60  # [ug/s]
Ti = 17.5456*60  # [s]
Td = 6.3897*60   # [s]
PID_param = [Kp, Ti, Td]

# Simulation parameters
ts = 10  # Sampling time (s)
TOF_cible = 10
pump_maximum_speed = 1200  # [ml/h]
atracurium_concentration = 10  # [mg/ml]
u_max = ((pump_maximum_speed*atracurium_concentration)/3600)*1000  # [ug/s]
u_min = 0       # [ug/s]

# Total duration of the simulation (in seconds)
T_simu = 250 * 60
N_simu = int(T_simu/ts)
# Control from 20 min to 150 min
T_control = 20 * 60
N_control = int(T_control/ts)
T_end = 150 * 60
N_end = int(T_end/ts)


for j in range(1, 13):
        # Load patient info
        Patient_info = Patient_table_atracurium[j-1]
        # Create patient
        Age = 30
        Height = 160
        Weight = Patient_info['Weight']
        sex = 1
        George = pas.Patient(patient_characteristic= [Age, Height, Weight, sex],atracurium_model_params= Patient_info, ts=ts)
        # Initial atracurium bolus of 500 ug/kg
        x0 = np.array([500/Patient_info['V1'], 0, 0, 0])
        George.initialized_at_given_state(x0)
        simulator = pas.Simulator(George)

        # Run simulation
        init = True
        for i in range(N_simu):

            if i >= N_control and i < N_end:
                if init:
                    PID_controller = PID_NMB(Kp, Ti, Td, Ts=ts, umax=u_max, umin=u_min, last_BIS=George.tof)
                    init = False
                uA = PID_controller.one_step(George.tof, TOF_cible)
                uA = min(u_max, max(0, uA))
            else:
                uA = 0

            simulator.one_step(input_atra= uA)


        p1.line(simulator.dataframe['Time']/60, simulator.dataframe['TOF'],
                        line_color="#006d43", legend_label='tof (%)')
        p2.line(simulator.dataframe['Time']/60, (simulator.dataframe['u_atra']*60)/Patient_info['Weight'],
                        line_color="#006d43", legend_label='Atracurium (ug/kg/min)')
        init = True

hline = Span(location=TOF_cible, dimension='width', line_color='black',
             line_dash='dashed', line_width=2)
p1.add_layout(hline)
p1.title.text = 'TOF response'
p2.title.text = 'Infusion rate'
p2.xaxis.axis_label = 'Time (min)'
grid = row(column(p1, p2))

show(grid)