Code covered by the BSD License  

Highlights from
RGDS_Practical_Guide

RGDS_Practical_Guide

by

 

15 Mar 2013 (Updated )

MATLAB routines for the book: "Development of Innovative Drugs via Modeling with MATLAB".

concmod(regimenType, numSubjects)
function concmod(regimenType, numSubjects)
%CONCMOD Model Tumor Growth
%   CONCMOD(REGIMENTYPE, NUMSUBJECTS) models tumor growth and creates
%   graphs of the tumor growth and epidermis over time. |REGIMENTYPE|
%   must be either 'intermittent' or 'continuous'. |NUMSUBJECTS| is
%   the number of subjects.
%
%   Example:
%   Model tumor growth with an intermittent regimen and 10 subjects 
%   concmod('intermittent',10)

% Check input arguments and provide default values if necessary
error(nargchk(0, 2, nargin));       %#ok<*NCHKN>

if nargin < 2
    numSubjects = 10;
end
if nargin < 1
    regimenType = 'intermittent';
end

% Reset the random number generator to the default
rng default; rand(100);

% Calculate dosing times and amount based on regimen
[doseTimes, doseAmount] = doseSchedule(regimenType);

% Set up figures for plotting results
tumorFigure = figure;
tumorAxes   = axes; set(tumorAxes,'FontSize',15)
xlabel('Time [h]')
ylabel('Number of tumor cells (relative to baseline)')
title('Tumor Growth')
xlim([0, doseTimes(end)]); hold on; grid on; 
epidermisFigure = figure;
epidermisAxes   = axes; set(epidermisAxes,'FontSize',15)
xlabel('Time [h]')
ylabel('% Change (relative to baseline)')
title('Epidermis')
xlim([0, doseTimes(end)]); hold on; grid on; ylim([50, 100]);

set([tumorFigure; epidermisFigure], ...
    'units', 'normalized', ...
    {'Position'}, {[0.1, 0.5, 0.3, 0.4]; [0.6, 0.5, 0.3, 0.4]});

% Simulate system for each subject
for subjectID = 1:numSubjects
    
    % Initialize parameters and set initial values
    p  = initializeParams;
    y0 = [doseAmount, p.c10, p.c20, p.n0, p.pc0, p.dc0, p.sc0];
    
    % Allocate variables to store results
    timePoints  = [];
    tumorGrowth = [];
    epidermis   = [];
    
    % Simulate system for each treatment period
    for dose = 1:(length(doseTimes)-1)
        
        % Set time interval for this treatment period
        tspan = [doseTimes(dose), doseTimes(dose+1)];  

        % Call Runge-Kutta method
        [t,y] = ode45(@derivatives, tspan, y0, [], p);  

        % Record values for plotting
        timePoints  = [timePoints ; t];                     %#ok             
        tumorGrowth = [tumorGrowth; y(:,4)/p.n0];           %#ok      
        epidermis   = [epidermis  ; 100*y(:,7)/p.sc0];      %#ok 

        % Reset initial values for next treatment period 
        % and add next dose
        y0    = y(end,:);
        y0(1) = y0(1) + doseAmount;  
    end
    
    % Plot results for this subject
    plot(tumorAxes, timePoints, tumorGrowth, 'Color', 'black')
    plot(epidermisAxes, timePoints, epidermis, 'Color', 'black')
    drawnow
end

% save graphs as TIFF file 
print(tumorFigure, '-r900', '-dtiff',      ...
    ['tumor', '_', regimenType])
print(epidermisFigure, '-r900', '-dtiff',  ...
    ['epidermis', '_', regimenType])

end


function dydt = derivatives( ~, y, p)
%DERIVATIVES Compute the right-hand side of the ODE.
%   DYDT = DERIVATIVES(T, Y, P) calculates |DYDT|, the right-hand
%   side of the ODE model, at points defined by the vector of
%   dependent variables |Y|, and time |T|, and with parameters |P|.

amount  = y(1);     % drug amount
c1      = y(2);     % parent drug concentration
c2      = y(3);     % active metabolite concentration
n       = y(4);     % tumor growth
pc      = y(5);     % cells in proliferative compartment
dc      = y(6);     % cells in differentiated compartment
sc      = y(7);     % cells in stratum corneum compartment 

% PK Model 
dAmountdt = -p.k01*amount;
dc1dt     = p.k01/p.V1*amount - (p.CL12/p.V1 + p.CL10/p.V1)*c1;   
dc2dt     = p.CL12/p.V2*c1 - p.CL20/p.V2*c2;      

% Efficacy Model
kkill = hillEffect(p.tumFactor*c2, 0, p.kkillMax, p.kkill50, 1); 
dndt  = (p.k0*log(p.nn00/n) - (p.k + kkill))*n; 

% Toxicity Model
DS    = hillEffect(p.skFactor*c2, 0, p.DSMax, p.DS50, 1);    
br1   = p.br0*(1-DS);   % birth rate;
dpcdt = br1 - p.kpc*pc;    
ddcdt = p.kpc*pc - p.kdc*dc;
dscdt = p.kdc*dc - p.ksc*sc;

% Derivatives vector of ODE system
dydt = [ dAmountdt;     % drug amount    
         dc1dt;         % parent drug concentration
         dc2dt;         % active metabolite concentration
         dndt;          % tumor growth
         dpcdt;         % changes in proliferative compartment    
         ddcdt;         % changes in differentiated compartment
         dscdt;    ];   % changes in stratum corneum compartment   
end


function [doseTimes, doseAmount] = doseSchedule(regimenType)
%DOSESCHEDULE Creates a vector of dosing times and amounts of drug
%   [DOSETIMES, DOSEAMOUNT] = DOSESCHEDULE(REGIMENTYPE) returns 
%   a vector |DOSETIMES| of dosing times and a scalar |DOSEAMOUNT| 
%   of the dosing amount. |REGIMENTYPE| must be either 'intermittent' 
%   or 'continuous'.

treatmentWeeks  = 12;
cycleWeeks      = 3;
daysInWeek      = 7;
hoursInDay      = 24;
initialTime     = 0;
% treatment time [hours]
endTime = treatmentWeeks * daysInWeek * hoursInDay;
% cycle time [hours]
cycleTime = cycleWeeks * daysInWeek * hoursInDay;

switch regimenType
    case 'intermittent' 
        % intermittent treatment
        dailyDoses = 2;         % BID
        doseAmount = 1500;      % 1 dose in [mg]
        timeOnDrug = cycleTime - daysInWeek * hoursInDay;  
    case 'continuous'
        % continuous (over 12 weeks) treatment
        dailyDoses = 2;         % BID
        doseAmount = 1000;      % 1 dose in [mg]
        timeOnDrug = cycleTime - hoursInDay/dailyDoses; 
    % case 'other' % placeholder for other regimen(s)
    otherwise
        messageID = 'Book:concmod:UnknownRegimen';
        messageStr = ['Unknown regimen ''%s''. The regimen must ',...
            'be either ''intermittent'' or ''continuous''.'];
        error(messageID, messageStr, regimenType);
end 

% create vector of treatment times
initialCycleTimes = ... 
    initialTime : cycleTime : (endTime - cycleTime);
doseTimesWithinCycle = ...
    initialTime : hoursInDay/dailyDoses : timeOnDrug; 
doseTimes = reshape(...
    repmat(doseTimesWithinCycle', 1, length(initialCycleTimes)) + ...
    repmat(initialCycleTimes, length(doseTimesWithinCycle), 1), ...
    1, length(initialCycleTimes) * length(doseTimesWithinCycle));
doseTimes = [doseTimes, endTime]; % add the end time of treatment
    
end


function params = initializeParams
%INITIALIZEPARAMS create initial values for model parameters.
%   PARAMS = INITIALIZEPARAMS returns a structure |PARAMS| containing
%   initial values for model parameters, including variability when
%   necessary.

% PK parameters 
pkCV = 0.3;
lim1 = 0.3;
lim2 = 5.0;

k01  = 0.7;     % first order elimination 
params.k01 = variability('varnorm', k01, k01*pkCV, 1, lim1);
CL12 = 10;      % metabolite clearances
params.CL12 = variability('varnorm', CL12, CL12*pkCV, 1, lim2);
CL10 = 80;      % parent drug clearances
params.CL10 = variability('varnorm', CL10, CL10*pkCV, 1, lim2);
CL20 = 60;      % metabolite elimination clearances
params.CL20 = variability('varnorm', CL20, CL20*pkCV, 1, lim2); 
V1  = 30;       % metabolite volumes of distribution
params.V1 = variability('varnorm', V1, V1*pkCV, 1, lim2);
V2  = 150;      % metabolite elimination volumes of distribution
params.V2 = variability('varnorm', V2, V2*pkCV, 1, lim2); 
params.c10 = 0;     % initial value: parent drug concentration
params.c20 = 0;     % initial value: active metabolite concentration   

% Efficacy parameters (tumor growth parameters)
params.k0 = 4.2e-5;       % doubling time 105 days
params.k = 0;             % natural elimination rate
params.tumFactor = 12;    % tumor factor   
params.nn00 = 1e12;       % tumor growth limiting value
params.n0 = 1e9;          % tumor growth initial value  
params.kkillMax = 0.05;   % max effect (Hill Eq.)     
kkill50 = 100;       % concentration linked to 50% of max effect
kkill50CV = 1.5;  
params.kkill50 = variability('varlog', ...
    kkill50, kkill50CV*kkill50, 1); 

% Toxicity parameters (epidermis)
params.skFactor = 8.0;   % toxicity factor
params.DSMax = 0.8;      % max effect (Hill Eq.)
DS50  = 10.0;       % concentration linked to 50% of max effect 
DS50CV = 0.8;
params.DS50 = variability('varlog', ...
    DS50, DS50CV*DS50, 1);

% proliferative compartment (pc)
ttpcd = 6;                  % cell cycle time in pc (days)      
ttpc = ttpcd*24;            % cell cycle time in pc (hours)
params.kpc = 1/ttpc;        % app. elimination rate at steady state 
growthFactor = 1.0;         % growth factor in pc 
pcTot = 27000;              % number of cells in pc
params.pc0 = growthFactor*pcTot;    % cells in pc at time 0
params.br0 = params.pc0/ttpc;       % initial birth rate

% differentiated compartment (dc)
ttdcd = 9;                  % transit time in dc (days)
ttdc = ttdcd*24;            % transit time in dc (hours)
params.kdc = 1/ttdc;        % app. elimination rate at steady state 
params.dc0 = params.br0*ttdc;     % cells in dc at time 0  

% stratum corneum compartment (sc)
ttscd = 7;                 % transit time in sc (days)
ttsc = ttscd*24;           % transit time in sc (hours)
params.ksc = 1/ttsc;       % app. elimination rate at steady state 
params.sc0 = params.br0*ttsc;     % cells in sc at time 0   

end


function x = variability(distribution, m, s, num, lim) 
%VARIABILITY Generate variability.
%   X = VARIABILITY(DISTRIBUTION, M, S, NUM, LIM) returns a vector 
%   of random numbers for variability, specified by |DISTRIBUTION|   
%   with the mean value |M| and standard deviation |S|. |NUM| is the
%   length of the vector |X|.
%   |DISTRIBUTION| can be 'varlog' or 'varnorm':  
%   'varlog' generates random numbers |X| ~ logN(|M|,|S|).
%   |M| stands for mean value of the lognormal distribution 
%   |S| - standard deviation of the lognormal distribution
%   'varnorm' generates normally distributed random numbers |X|,  
%   where |X| ~ N(|M|,|S|), right-censored by |LIM|. 

switch distribution
    case 'varlog'
        %lognormal distribution
        mu    = log(m^2/sqrt(m^2 + s^2));   % where: N(|mu|,|sigma|)
        sigma = sqrt(log(1 + (s/m)^2));     % where: N(|mu|,|sigma|) 
        x     = lognrnd(mu,sigma, 1, num);
    case 'varnorm'
        %lognormal distribution
        mu    = m;
        sigma = s;
        x     = max(normrnd(mu, sigma, 1, num), lim);
    otherwise 
        % placeholder for other distributions
end

end


function E = hillEffect(c, E0, Emax, EC50, n)
%HILLEFFECT Compute drug effect based on the Hill equation.
%   E = HILLEFFECT(C, E0, EMAX, EC50, N) computes drug effect |E|, 
%   based on the Hill equation as a function of concentration |C|,
%   and with the following concentration-response parameters:
%   |E0| - baseline response, |EMAX| - maximum effect, 
%   |EC50| - concentration related to 50% of max. effect, 
%   and |N| - Hill coefficient of sigmoidicity.


E = E0 + Emax.*c.^n./(EC50+c.^n); 

end

Contact us