runPropofolSystem.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% runPropofolSystem.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Declares system parameters and runs PropofolSystem.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all
clc

% choose parameters and bounds
% from Eleveld et al. 

    % deposition volumes [mL]
    V1 = 6.28 * 1000;               % plasma 
    V2 = 25.5 * 1000;               % fast periph.
    V3 = 273 * 1000;                % slow periph.

    % clearances [mL⁻¹ min⁻¹]
    CL = 1.79 * 1000;               % plasma
    Q2 = 1.75 * 1000;               % fast periph.
    Q3 = 1.11 * 1000;               % slow periph.

    % transport coefficients [min⁻¹]
    k10 = CL/V1;                    % plasma -> elimination
    k12 = Q2/V1;                    % fast periph. -> plasma
    k13 = Q3/V1;                    % slow periph. -> plasma
    k21 = Q2/V2;                    % plasma -> fast periph.
    k31 = Q3/V3;                    % plasma -> slow periph. 
    ke0 = 0.146;                    % plasma -> effect site

    % pharmacodynamic parameters
    ce50 = 3.08;                    % ce50 at 50% BIS [μg mL⁻¹]
    BIS0 = 93;                      % undrugged BIS [%]
    gamma = 1.47;                   % Hill cooperativity constant [-]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set up controllers

    uCap = 5e4;                     % actuation limit, which is
                                    % around maintenance bolus for 
                                    % this patient; lower than induction
                                    % bolus [μg mL⁻¹ min⁻¹]
    BISmax = 60;                    % maximum safe BIS [%]
    BISmin = 40;                    % minimum safe BIS [%]
    BISdes = 50;                    % desired BIS [%]

% exponential controller

    a1 = sqrt(10); a2 = sqrt(10);   % eigenvalues         
    wExp = 5;                       % slack var. pareto weight for min

% graceful controller

    omega = 4; zeta = 1;        % natural freq., damping ratio
    wGrc1 = 1e8; wGrc2 = 1e2;       % slack var. pareto weights
                                    % for min, max
    ceMinGrc = -0.001;              % lowermost barrier  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set up interval estimator

    d = 0.05;                       % disturbance magnitude bound
                                    % [μg mL⁻¹]
    V = 0.5;                        % sensor noise magnitude bound 
                                    % [μg mL⁻¹]
                                    % note that this corresponds to
                                    % about 10-20% BIS error points
    initError = 0.3;                % initial estimation error:
                                    % for the given graceful pareto 
                                    % weights, one state blows up at 0.5

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initialize patients
    
    x0 = [0; 0; 0; 0]; xLoEst0 = x0; xHiEst0 = x0; 
    x0Color = '#000000'; 
    patient0 = PropofolSystem(x0Color, x0, ...
                k10, k12, k13, k21, k31, ke0, V1, ...
                BIS0, ce50, gamma, ...
                BISdes, BISmin, BISmax, ...
                uCap, a1, a2, wExp, ...
                ceMinGrc, zeta, omega, wGrc1, wGrc2, ...
                xLoEst0, xHiEst0, d, V);
    patient0.plotCharacteristicCurve(); 

    x0a = [patient0.ceMin-2; 0; 0; patient0.ceMin-1]; 
    xLoEst0a = [x0a(1)-initError; 0; 0; x0a(4)-initError];
    xHiEst0a = [x0a(1)+initError; 0; 0; x0a(4)+initError];
    x0aColor = '#666666'; 
    patient0a = PropofolSystem(x0aColor, x0a, ...
                k10, k12, k13, k21, k31, ke0, V1, ...
                BIS0, ce50, gamma, ...
                BISdes, BISmin, BISmax, ...
                uCap, a1, a2, wExp, ...
                ceMinGrc, zeta, omega, wGrc1, wGrc2, ...
                xLoEst0a, xHiEst0a, d, V);

    x0b = [patient0.ceMin-2; 0; 0; patient0.ceMin-initError]; 
    xLoEst0b = [x0b(1)-initError; 0; 0; x0b(4)-initError];
    xHiEst0b = [x0b(1)+initError; 0; 0; x0b(4)+initError];
    x0bColor = '#b7b7b7';
    patient0b = PropofolSystem(x0bColor, x0b, ...
                k10, k12, k13, k21, k31, ke0, V1, ...
                BIS0, ce50, gamma, ...
                BISdes, BISmin, BISmax, ...
                uCap, a1, a2, wExp, ...
                ceMinGrc, zeta, omega, wGrc1, wGrc2, ...
                xLoEst0b, xHiEst0b, d, V);

    x01 = [patient0.ceMax+initError; 0; 0; patient0.ceMax-initError]; 
    xLoEst01 = [x01(1)-initError; 0; 0; x01(4)-initError];
    xHiEst01 = [x01(1)+initError; 0; 0; x01(4)+initError];
    x1Color = '#66da54';
    patient1 = PropofolSystem(x1Color, x01, ...
                k10, k12, k13, k21, k31, ke0, V1, ...
                BIS0, ce50, gamma, ...
                BISdes, BISmin, BISmax, ...
                uCap, a1, a2, wExp, ...
                ceMinGrc, zeta, omega, wGrc1, wGrc2, ...
                xLoEst01, xHiEst01, d, V);

    x02 = [patient0.ceMax-1; 0; 0; patient0.ceMax-initError]; 
    xLoEst02 = [x02(1)-initError; 0; 0; x02(4)-initError];
    xHiEst02 = [x02(1)+initError; 0; 0; x02(4)+initError];
    x2Color = '#cada54';
    patient2 = PropofolSystem(x2Color, x02, ...
               k10, k12, k13, k21, k31, ke0, V1, ...
                BIS0, ce50, gamma, ...
                BISdes, BISmin, BISmax, ...
                uCap, a1, a2, wExp, ...
                ceMinGrc, zeta, omega, wGrc1, wGrc2, ...
                xLoEst02, xHiEst02, d, V);

    x03 = [patient0.ceMax; 0; 0; patient0.ceMax-0.7]; 
    xLoEst03 = [x03(1)-initError; 0; 0; x03(4)-initError];
    xHiEst03 = [x03(1)+initError; 0; 0; x03(4)+initError];
    x3Color = '#dcd62b';
    patient3 = PropofolSystem(x3Color, x03, ...
                k10, k12, k13, k21, k31, ke0, V1, ...
                BIS0, ce50, gamma, ...
                BISdes, BISmin, BISmax, ...
                uCap, a1, a2, wExp, ...
                ceMinGrc, zeta, omega, wGrc1, wGrc2, ...
                xLoEst03, xHiEst03, d, V);

    x04 = [patient0.ceDes-initError; 0; 0; patient0.ceDes-initError]; 
    xLoEst04 = [x04(1)-initError; 0; 0; x04(4)-initError];
    xHiEst04 = [x04(1)+initError; 0; 0; x04(4)+initError];
    x4Color = '#dda010';
    patient4 = PropofolSystem(x4Color, x04, ...
                k10, k12, k13, k21, k31, ke0, V1, ...
                BIS0, ce50, gamma, ...
                BISdes, BISmin, BISmax, ...
                uCap, a1, a2, wExp, ...
                ceMinGrc, zeta, omega, wGrc1, wGrc2, ...
                xLoEst04, xHiEst04, d, V);

    x05 = [patient0.ceMin-1.85; 0; 0; patient0.ceMin+0.12]; 
    xLoEst05 = [x05(1)-initError; 0; 0; x05(4)-initError];
    xHiEst05 = [x05(1)+initError; 0; 0; x05(4)+initError];
    x5Color = '#dd5510';
    patient5 = PropofolSystem(x5Color, x05, ...
                k10, k12, k13, k21, k31, ke0, V1, ...
                BIS0, ce50, gamma, ...
                BISdes, BISmin, BISmax, ...
                uCap, a1, a2, wExp, ...
                ceMinGrc, zeta, omega, wGrc1, wGrc2, ...
                xLoEst05, xHiEst05, d, V);

    x06 = [patient0.ceMin-1; 0; 0; patient0.ceMin+0.04]; 
    xLoEst06 = [x06(1)-initError; 0; 0; x06(4)-initError];
    xHiEst06 = [x06(1)+initError; 0; 0; x06(4)+initError];
    x6Color = '#980000';
    patient6 = PropofolSystem(x6Color, x06, ...
                k10, k12, k13, k21, k31, ke0, V1, ...
                BIS0, ce50, gamma, ...
                BISdes, BISmin, BISmax, ...
                uCap, a1, a2, wExp, ...
                ceMinGrc, zeta, omega, wGrc1, wGrc2, ...
                xLoEst06, xHiEst06, d, V);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% simulation settings

    tStart = 0; tEnd = 70;          % [min]    
    nTimes = 8000;
    tHist = zeros(1, nTimes); 
    runObserver = true;             % perform estimation in simulation
    useObserver = true;             % feed estimates into graceful 
                                    % controller

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% simulate patients with exponential controller

    % graceful = false;
    % disp('RUNNING EXPONENTIAL CONTROLLER...'); 
    %     [patient0, patient0a, patient0b, ...
    %     patient1, patient2, patient3, patient4, patient5, patient6] = ...
    % simulatePatients(patient0, patient0a, patient0b, ...
    %     patient1, patient2, patient3, patient4, patient5, patient6, ...
    %     tStart, tEnd, nTimes, graceful, runObserver, useObserver);
    % plotInputHistories(patient0, patient0a, patient0b, patient1, ...
    %     patient2, patient3, patient4, patient5, patient6, graceful);
    % plotBISHistories(patient0, patient0a, patient0b, patient1, ...
    %     patient2, patient3, patient4, patient5, patient6, graceful) 
    % plotPhasePortraits(patient0, patient0a, patient0b, patient1, ...
    %     patient2, patient3, patient4, patient5, patient6, graceful);

% simulate patients with graceful controller

    graceful = true; 
    disp('RUNNING GRACEFUL CONTROLLER...'); 
    if runObserver && ~useObserver
        disp('Estimation ON, NOT USED by controller.');
    elseif runObserver && useObserver
        disp('Estimation ON, USED by controller.');
    elseif ~runObserver
        disp('Estimation OFF.');
    end
        [patient0, patient0a, patient0b, ...
        patient1, patient2, patient3, patient4, patient5, patient6] = ...
    simulatePatients(patient0, patient0a, patient0b, ...
        patient1, patient2, patient3, patient4, patient5, patient6, ...
        tStart, tEnd, nTimes, graceful, runObserver, useObserver);
    printPatientMetrics(patient0, patient0a, patient0b, patient1, ...
        patient2, patient3, patient4, patient5, patient6);
    plotInputHistories(patient0, patient0a, patient0b, patient1, ...
        patient2, patient3, patient4, patient5, patient6, graceful);
    plotBISHistories(patient0, patient0a, patient0b, patient1, ...
        patient2, patient3, patient4, patient5, patient6, graceful) 
    plotPhasePortraits(patient0, patient0a, patient0b, patient1, ...
        patient2, patient3, patient4, patient5, patient6, graceful);
    if runObserver
        plotEstimation(patient0, patient0a, patient0b, patient1, ...
            patient2, patient3, patient4, patient5, patient6); 
        disp('OBSERVER DYNAMICS:')
        fprintf('A - LC =\n')
        disp(patient0.A - patient0.L*patient0.C); 
        fprintf('eig(A - LC) =\n')
        disp(eig(patient0.A - patient0.L*patient0.C));
    end

   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [patient0, patient0a, patient0b, ...
        patient1, patient2, patient3, patient4, patient5, patient6] = ...
    simulatePatients(patient0, patient0a, patient0b, ...
        patient1, patient2, patient3, patient4, patient5, patient6, ...
        tStart, tEnd, nTimes, graceful, runObserver, useObserver)
    disp('Running Initial Condition 0...'); 
    patient0 = patient0.simulate(graceful, runObserver, useObserver, ...
        tStart, tEnd, nTimes);
    disp('Running Initial Condition 0a...'); 
    patient0a = patient0a.simulate(graceful, runObserver, useObserver, ...
        tStart, tEnd, nTimes);
    disp('Running Initial Condition 0b...'); 
    patient0b = patient0b.simulate(graceful, runObserver, useObserver, ...
        tStart, tEnd, nTimes);
    disp('Running Initial Condition 1...'); 
    patient1 = patient1.simulate(graceful, runObserver, useObserver, ...
        tStart, tEnd, nTimes);
    disp('Running Initial Condition 2...'); 
    patient2 = patient2.simulate(graceful, runObserver, useObserver, ...
        tStart, tEnd, nTimes);
    disp('Running Initial Condition 3...'); 
    patient3 = patient3.simulate(graceful, runObserver, useObserver, ...
        tStart, tEnd, nTimes);
    disp('Running Initial Condition 4...'); 
    patient4 = patient4.simulate(graceful, runObserver, useObserver, ...
        tStart, tEnd, nTimes);
    disp('Running Initial Condition 5...'); 
    patient5 = patient5.simulate(graceful, runObserver, useObserver, ...
        tStart, tEnd, nTimes);
    disp('Running Initial Condition 6...'); 
    patient6 = patient6.simulate(graceful, runObserver, useObserver, ...
        tStart, tEnd, nTimes);
    fprintf('\n'); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% print per-IC summary metrics: induction time (first time BIS reaches
% the desired level), min/max BIS over the run, and total propofol dose
% administered (integral of the infusion rate over time)

function printPatientMetrics(patient0, patient0a, patient0b, ...
    patient1, patient2, patient3, patient4, patient5, patient6)
    labels = {'0', '0a', '0b', '1', '2', '3', '4', '5', '6'};
    patients = {patient0, patient0a, patient0b, ...
        patient1, patient2, patient3, patient4, patient5, patient6};
    fprintf('\n%-4s %16s %10s %10s %16s\n', ...
        'IC', 'Induction [min]', 'BISmin', 'BISmax', 'Total dose [ug]');
    for i = 1:length(patients)
        p = patients{i};
        inductionIdx = find(p.BISHist <= p.BISdes, 1, 'first');
        if isempty(inductionIdx)
            inductionTime = NaN;
        else
            inductionTime = p.tHist(inductionIdx);
        end
        totalDose = trapz(p.tHist, p.uHist);
        fprintf('%-4s %16.2f %10.2f %10.2f %16.1f\n', labels{i}, ...
            inductionTime, min(p.BISHist), max(p.BISHist), totalDose);
    end
    fprintf('\n');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plotInputHistories(patient0, patient0a, patient0b, ...
    patient1, patient2, patient3, patient4, patient5, patient6, graceful)
    figure('Color', [1 1 1]); 
    hold on
    yline(patient0.uCap, 'k--', 'maintenance bolus', ...
        'FontSize', 15)
    plot(patient0.tHist, patient0.uHist, ...
        'LineWidth', 2, 'Color', patient0.color);
    plot(patient0a.tHist, patient0a.uHist, ...
        'LineWidth', 2, 'Color', patient0a.color);
    plot(patient0b.tHist, patient0b.uHist, ...
        'LineWidth', 2, 'Color', patient0b.color);
    plot(patient1.tHist, patient1.uHist, ...
        'LineWidth', 2, 'Color', patient1.color);
    plot(patient2.tHist, patient2.uHist, ...
        'LineWidth', 2, 'Color', patient2.color);
    plot(patient3.tHist, patient3.uHist, ...
        'LineWidth', 2, 'Color', patient3.color);
    plot(patient4.tHist, patient4.uHist, ...
        'LineWidth', 2, 'Color', patient4.color);
    plot(patient5.tHist, patient5.uHist, ...
        'LineWidth', 2, 'Color', patient5.color);
    plot(patient6.tHist, patient6.uHist, ...
        'LineWidth', 2, 'Color', patient6.color);
    if graceful == false
        title('Control input history — Exponential controller'); 
    else 
        title('Control input history — Graceful controller'); 
    end  
    xlabel('Time [min]'); 
    ylabel('Infusion rate [μg mL⁻¹ min⁻¹]');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plotBISHistories(patient0, patient0a, patient0b, patient1, ...
        patient2, patient3, patient4, patient5, patient6, graceful)

figure('Color', [1 1 1]); 
    hold on
    plot(patient0.tHist, patient0.BISHist, ...
        'LineWidth', 2, 'Color', patient0.color);
    plot(patient0a.tHist, patient0a.BISHist, ...
        'LineWidth', 2, 'Color', patient0a.color);
    plot(patient0b.tHist, patient0b.BISHist, ...
        'LineWidth', 2, 'Color', patient0b.color);
    plot(patient1.tHist, patient1.BISHist, ...
        'LineWidth', 2, 'Color', patient1.color);
    plot(patient2.tHist, patient2.BISHist, ...
        'LineWidth', 2, 'Color', patient2.color);
    plot(patient3.tHist, patient3.BISHist, ...
        'LineWidth', 2, 'Color', patient3.color);
    plot(patient4.tHist, patient4.BISHist, ...
        'LineWidth', 2, 'Color', patient4.color);
    plot(patient5.tHist, patient5.BISHist, ...
        'LineWidth', 2, 'Color', patient5.color);
    plot(patient6.tHist, patient6.BISHist, ...
        'LineWidth', 2, 'Color', patient6.color);
    yline(patient0.BISmax, 'k--', ...
        'upper bound for surgery', 'FontSize', 15);
    yline(patient0.BISmin, 'k--', ...
        'lower bound for surgery', 'FontSize', 15);
    yline(patient0.BISdes, 'k--', ...
        'desired', 'FontSize', 15);
    if graceful == false
        title('Bispectral index history — Exponential controller'); 
    else 
        title('Bispectral index history — Graceful controller');
    end  
    xlabel('Time [min]'); 
    ylabel('Bispectral index [%]');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plotPhasePortraits(patient0, patient0a, patient0b, ...
    patient1, patient2, patient3, patient4, patient5, ...
    patient6, graceful)
    figure('Color', [1 1 1]); 
    hold on
    grid
    xlim([-0.2 4]); ylim([-0.5 1]);
    if graceful == false
        ylim([-0.3 1.12]);
        fill([patient0.ceMin 20 20 patient0.ceMin], ...
            [-10 -10 10 10], [0.9 0.9 0.9], ...
            'EdgeColor', 'none', 'FaceAlpha', 0.6); 
        ceVals = linspace(-20, 20, 500);
        bound = -1 * patient0.a1 * (ceVals - patient0.ceMin);
        fill([ceVals fliplr(ceVals)], ...
            [bound fliplr(ones(size(ceVals)) * 10)], ...
            [0.6 0.8 1], 'EdgeColor', 'none', 'FaceAlpha', 0.2);
        plot(ceVals, bound, 'k--', 'LineWidth', 1); 
        xline(patient0.ceMin, 'k--', 'lower bound for surgery', ...
            'LineWidth', 1, 'FontSize', 15);
        xline(patient0.ceMax, 'k--', 'upper bound for surgery', ...
            'LineWidth', 1, 'FontSize', 15);
        xline(patient0.ceDes, 'k--', 'desired concentration', ...
            'LineWidth', 1, 'FontSize', 15); 
        yline(0, 'k', 'LineWidth', 0.5);
    else
        fill([patient0.ceMin 20 20 patient0.ceMin], ...
            [-10 -10 10 10], [0.9 0.9 0.9], ...
            'EdgeColor', 'none', 'FaceAlpha', 0.6);
        xline(patient0.ceMin, 'k--', 'lower bound for surgery', ...
            'LineWidth', 1, 'FontSize', 15);
        xline(patient0.ceMax, 'k--', 'upper bound for surgery', ...
            'LineWidth', 1, 'FontSize', 15);
        xline(patient0.ceG1, 'k--', ['lower zeroing barrier ' ...
            '(desired)'],'LineWidth', 1, 'FontSize', 15); 
        xline(patient0.ceG2, 'k--', 'upper zeroing barrier', ...
            'LineWidth', 1, 'FontSize', 15);
    end
    plot(patient0.xHist(4,1), patient0.ceDotHist(1), ...
        'Marker', '.', 'MarkerSize', 18, ...
        'Color', patient0.color);
    plot(patient0.xHist(4,:), patient0.ceDotHist, ...
        'LineWidth', 2, 'Color', patient0.color);
    plot(patient0a.xHist(4,1), patient0a.ceDotHist(1), ...
        'Marker', '.', 'MarkerSize', 18, ...
        'Color', patient0a.color);
    plot(patient0a.xHist(4,:), patient0a.ceDotHist, ...
        'LineWidth', 2, 'Color', patient0a.color);
    plot(patient0b.xHist(4,1), patient0b.ceDotHist(1), ...
        'Marker', '.', 'MarkerSize', 18, ...
        'Color', patient0b.color);
    plot(patient0b.xHist(4,:), patient0b.ceDotHist, ...
        'LineWidth', 2, 'Color', patient0b.color);
    plot(patient1.xHist(4,1), patient1.ceDotHist(1), ...
        'Marker', '.', 'MarkerSize', 18, ...
        'Color', patient1.color);
    plot(patient1.xHist(4,:), patient1.ceDotHist, ...
        'LineWidth', 2, 'Color', patient1.color);
    plot(patient2.xHist(4,1), patient2.ceDotHist(1), ...
        'Marker', '.', 'MarkerSize', 18, ...
        'Color', patient2.color);
    plot(patient2.xHist(4,:), patient2.ceDotHist, ...
        'LineWidth', 2, 'Color', patient2.color);
    plot(patient3.xHist(4,1), patient3.ceDotHist(1), ...
        'Marker', '.', 'MarkerSize', 18, ...
        'Color', patient3.color);
    plot(patient3.xHist(4,:), patient3.ceDotHist, ...
        'LineWidth', 2, 'Color', patient3.color);
    plot(patient4.xHist(4,1), patient4.ceDotHist(1), ...
        'Marker', '.', 'MarkerSize', 18, ...
        'Color', patient4.color);
    plot(patient4.xHist(4,:), patient4.ceDotHist, ...
        'LineWidth', 2, 'Color', patient4.color);
    plot(patient5.xHist(4,1), patient5.ceDotHist(1), ...
        'Marker', '.', 'MarkerSize', 18, ...
        'Color', patient5.color);
    plot(patient5.xHist(4,:), patient5.ceDotHist, ...
        'LineWidth', 2, 'Color', patient5.color);
    plot(patient6.xHist(4,1), patient6.ceDotHist(1), ...
        'Marker', '.', 'MarkerSize', 18, ...
        'Color', patient6.color);
    plot(patient6.xHist(4,:), patient6.ceDotHist, ...
        'LineWidth', 2, 'Color', patient6.color);
    if ~graceful
        title(['Effect site concentration phase portrait' ...
            '— Exponential controller']); 
        legend('1st-order safe set', '2nd-order safe set', ...
            'FontSize', 15);
    else 
        title(['Effect site concentration phase portrait' ...
            ' — Graceful controller']); 
        legend('Safe set'); 
    end  
    xlabel('Effect site concentration [μg mL⁻¹]'); 
    ylabel('Effect site concentration time derivative [μg mL⁻¹ min⁻¹]');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plotEstimation(patient0, patient0a, patient0b, ...
    patient1, patient2, patient3, patient4, patient5, patient6)
    patient0.plotEstimation('0'); 
    patient0a.plotEstimation('0a'); 
    patient0b.plotEstimation('0b'); 
    patient1.plotEstimation('1'); 
    patient2.plotEstimation('2'); 
    patient3.plotEstimation('3'); 
    patient4.plotEstimation('4'); 
    patient5.plotEstimation('5'); 
    patient6.plotEstimation('6'); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%