Mandel – Propofol Masui 2009

MATLAB Code

Implements the propofol model of Masui et al.[1]

Contributed by Jeff E Mandel

%% MATLAB implementation of Early Phase Propofol model 
% This routine implements the model of early phase propofol
% pharmacokinetics described in
%
% <http://journals.lww.com/anesthesiology/Abstract/2009/10000/Early_Phase_Pharmacokinetics_but_Not.24.aspx 
% Masui, K. and Kira, M. and Kazama, T. and 
% Hagihira, S. and Mortier, EP and Struys, MM  "Early Phase
% Pharmacokinetics but Not Pharmacodynamics Are Influenced by Propofol
% Infusion Rate", Anesthesiology 2009; 111:805-17.>
%
% Two figures from the paper have been reproduced; the PDF was used to
% extract PNG files of the figures, and graphs were overlaid from the
% MATLAB routines to demonstrate the equivalence of results.
%
% While the MATLAB program reproduces the results of Masui, it is unlikely
% that the offset of the TRANSIT model is correct. The user is cautioned to
% avoid making clinical decisions based on a precise interpretation of the
% this model.
%
%     Copyright (C) 2009 Jeff E Mandel MD MS
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%% testMasui
% Generates the graphics

function testMasui

timestep = 1/120; % 0.5 second calculation interval

% Fig 5F - Effect of weight variations on model
tt = 210; %seconds
t= 0:timestep:tt/60;
dose = zeros(size(t));
weight = [40;55;70];
nw = length(weight);
Cp = zeros(length(t),nw);
for i = 1:nw
    rate = 6 * weight(i);   % 60 mg/kg/h -> ml/h
    age = 50;
    [sys] = masui_propofol(weight(i), age, rate);
    % Infuse 60 mg/kg/h for 2 minutes
    dose(1:2/timestep) = 1000*weight(i); % µg/min
    Cp(:,i) = lsim(sys,dose,t);
end
pngName='fig5F.png';
superimposePlot(pngName,0);
plot(t*60,Cp,'LineWidth',1);
set(gca,'XLim',[0 tt]);
set(gca,'YLim',[0 20]);
set(gca,'Color', 'none');
set(gca,'Visible', 'off');
set(gcf, 'MenuBar', 'none');
title('Figure 5F');
snapnow;

% Fig 5B 
tt = 50;
t= 0:timestep:tt/60;
rate = [50;250;1000]; %ml/h
age = 50;
weight = 60;
nw = length(rate);
Cp = zeros(length(t),nw);
for i = 1:nw
    [sys] = masui_propofol(weight, age, rate(i));
    dose = ones(size(t))*rate(i)*10000/60;    %ml/h -> µg/min
    Cp(:,i) = lsim(sys,dose,t);
end

pngName='fig5B.png';
superimposePlot(pngName,0);
plot(t*60,Cp,'LineWidth',1);
set(gca,'XLim',[0 tt]);
set(gca,'YLim',[0 5]);
set(gca,'Color', 'none');
set(gca,'Visible', 'off');
set(gcf, 'MenuBar', 'none');
title('Figure 5B');
snapnow;
 

end

%% masui_propofol
% Creates a state-space model using the values in Table 3

function [sys] = masui_propofol(weight, age, rate)
% Weight in kg
% Volume in L
% rate is ml/h of 10 mg/ml propofol

Th1 = 0.0104;
Th2 = 2.20;
Th3 = 3.2;
Th4 = 4.21;
Th5 = 0.399;
Th6 = 29.2;
Th7 = 0.0471;

V1 = weight * Th1;
V2 = Th2;
Cl1 = Th3;  %L/min
Cl2 = Th4;

volume = [V1;V2];
clearance = [Cl1;Cl2];

LAG = age * Th5/60;    %minutes

%% 
% 6th order lag, see
%
% <http://www.springerlink.com/content/wm7763462310778q/ Savic RM, Jonker DM, Kerbusch T, Karlsson MO: Implementation of a transit
% compartment model for describing drug absorption in pharmacokinetic
% studies. J Pharmacokinet Pharmacodyn 2007; 34:711-26>
Ktr = Th6 + (rate-345)*Th7; %min^-1
n = 6;
MTT = n+1/Ktr;  % mean transit time

% 6th order lag specified as a transfer function and converted to state
% space form
pump = ss(tf(Ktr,[1 Ktr])^n);   
sys = mam2ss(clearance,volume);
sys.inputd = LAG;
sys = sys*pump;

end

%% mam2ss
% Returns a state space model from vectors of clearance
% and compartment volumes. All states are drug amount
% in the compartment; to get concentration, divide by compartment
% volume in the C matrix. Note that this code can handle an
% arbitrary number of compartments, but assumes the first compartment is
% the central compartment (mam is short for mamillary). 

function [sys] = mam2ss(clearance,volume)

n = length(clearance);
k1 = [clearance./volume(1)];
k2 = [clearance(2:n)./volume(2:n)];
A = [-sum(k1) k2;transpose(k1(2:n)) -diag(k2)];
B = [1;zeros(n-1,1)];           %Drug enters central compartment
C = [1/volume(1) zeros(1,n-1)]; %Central compartment concentration.

C = C./1000; % All compartment volumes are in L, convert to ml.
D=0;

sys = ss(A,B,C,D);
end

%% superimposePlot
% Create an axis with normalized coordinates into which we can draw the
% PNG. We fill the entire figure with the PNG at its native size. 
% We add a second axis above the image, make it transparent, and use this
% to get the coordinates.


function superimposePlot(pngName, first)

axes('units','normalized','position',[0 0 1 1]);
imshow(pngName);
truesize;
axis image;
hb = axes('units','normalized','position',[0 0 1 1], 'Layer', 'top');
axis off;
set(gca,'Color', 'none');

% name for a .mat file to save the rect
saveFname = regexprep(pngName,'\.png','Rect');

if first
    % The user drags a rectangle from the top left to bottom right of the
    % axis in the PNG. We transform this to the coordinate space of the hb
    % axis
    rect = getrect(hb);
    axlim = axis(hb);
    axwidth = diff(axlim(1:2));
    axheight = diff(axlim(3:4));
	rect(1) = (rect(1) - axlim(1)) / axwidth;
	rect(2) = (rect(2) - axlim(3)) / axheight;
	rect(3) = rect(3) / axwidth;
	rect(4) = rect(4)/ axheight;
    save(saveFname, 'rect');
else
    S=load(saveFname, 'rect');
    rect = S.rect;
end
axes('position', rect, 'Layer', 'top');

end