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',14)
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',14)
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 |Y| and
% |T|, and with parameters |P|. |Y| stands for the dependent
% variables.
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.kCkill50, 1);
dndt = (p.k0*log(p.nn00/n) - (p.k + kkill))*n;
% Toxicity Model
tox = hillEffect(p.toxFactor*c2, 0, p.toxMax, p.toxC50, 1);
br1 = p.br0*(1-tox); % 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.)
kCkill50 = 100; % concentration linked to 50% of max effect
kCkill50CV = 1.5;
params.kCkill50 = variability('varlog', ...
kCkill50, kCkill50CV*kCkill50, 1);
% Toxicity parameters (epidermis)
params.toxFactor = 8.0; % toxicity factor
params.toxMax = 0.8; % max effect (Hill Eq.)
toxC50 = 10.0; % concentration linked to 50% of max effect
toxC50CV = 0.8;
params.toxC50 = variability('varlog', ...
toxC50, toxC50CV*toxC50, 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 lognormally distribution
% |S| - standard deviation of the lognormally 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
