% Supplementary File 1: MATLAB code for % The Effect of Combination Therapy on the Evolution of Resistance in % Cancer: The Theoretical Merits of Synergism and Antagonism % Copyright: Elysia C. Saputra, September 26, 2017 % NOTE: to run on MATLAB, change file extension from .txt to .m % This script runs a Monte Carlo simulation of the stochastic evolutionary % tumor cell population model under dual combination treatment, % according to either the Constant-Dose or the Constant-Efficacy Method, % where all combination treatments are dosed symmetrically. T_lethal (time % until a lethal tumour burden) and T_RR (time until fully resistant RR % cells emerge) are determined from the stochastic simulations. The % probability distribution of T_RR = t from the stochastic simulations is % then compared to the probability distribution curve of T_RR generated % from the theoretical approximation of T_RR. clear;clc; %% ======================================================================= % Part 1: Stochastic evolutionary tumour cell population model for % simulating Constant-Dose Method or Constant-Efficacy Method % -------- % Initializing simulation parameters MC = 1000; % number of experiments for Monte Carlo simulation runTime = 350; % maximum run time % Parameters for population growth dynamics p = 0.95; % proliferation rate of all subpopulations in the absence of treatment u = 10^(-6); % phenotype alteration rate a_b = 0.12; % basal death rate a_c = 0.2537; % scaling factor for drug effect on apoptosis PP_c = 1; % scaling factor for drug effect on proliferation N_init = [10^9 0 0 0]; % number of cells in each subpopulation at the % start of each evolution cycle % N_init(1) = initial number of SS cells % N_init(2) = initial number of RS cells % N_init(3) = initial number of SR cells % N_init(4) = initial number of RR cells % Parameters for treatment A50 = 10; % ED50 of drug A B50 = 10; % ED50 of drug B E_max = 1; % maximum value of PPi R = A50/B50; % potency ratio constant alphaRef = 0; % alpha of reference treatment doseRef = 11; % dose of reference treatment treatment = 1; % treatment = 1 for Constant-Dose Method % treatment = 2 for Constant-Efficacy Method alpha = 0; % alpha parameter of the treatment to be evaluated % -------- % Defining dose limits if treatment == 1 % Constant-Dose Method doseA = doseRef; % dose of drug A administered doseB = doseRef; % dose of drug B administered elseif treatment == 2 % Constant-Efficacy Method % Equivalent dose of drug A with the same effect as the % reference drug combination D_A = doseRef + doseRef*R + alphaRef*doseRef^2/A50; % Calculate the dose of combination treatment to be simulated syms dA solx = solve(D_A == dA*(1+R) + alpha*dA^2/B50, dA); if alpha >= 0 doseA = max(eval(solx)); doseB = doseA; else doseA = min(eval(solx)); doseB = doseA; end end % Calculating drug effect parameters (DEP) % Dose sensed by cell: % SS = doseA + doseB % RS = 0.1*doseA + doseB % SR = doseA + 0.1*doseB % RR = 0.1*doseA + 0.1*doseB S_A = [1 0.1 1 0.1]; % sensitivity of each subpopulation to drug A, for [SS RS SR RR] S_B = [1 1 0.1 0.1]; % sensitivity of each subpopulation to drug B, for [SS RS SR RR] % Calculating normalized drug effect % D_A_i = equivalent dose of drug A effectively sensed by subpopulation i D_A_i = S_A*doseA + S_B*doseB*R + alpha*S_A.*S_B*doseA*doseB/B50; % E_i = normalized drug effect for subpopulation i E_i = D_A_i./(A50+D_A_i); % Calculating drug effect parameters % PP_i = proliferation probability of subpopulation i PP_i = 1 - PP_c*E_i; % a_i = drug-induced death rate of subpopulation i a_i = a_c*E_i; % effective growth rate k_i = (1+p*PP_i).*(1-a_b-a_i); % effective growth rate of subpopulation i % -------- % Simulation model t = 0:runTime; % time (unit: generation) % [length(t) x n] matrices storing the growth dynamics of each % subpopulation for each of the n experiments in the Monte Carlo run SS = zeros(runTime+1, MC); RS = zeros(runTime+1, MC); SR = zeros(runTime+1, MC); RR = zeros(runTime+1, MC); total = zeros(runTime+1, MC); T_RR = NaN(1, MC); % vector storing T_RR for each of the n experiments % in the Monte Carlo run T_lethal = NaN(1, MC); % vector storing T_lethal for each of the n experiments % in the Monte Carlo run % To improve the efficiency of the simulation, pre-generate a matrix of % Poisson numbers with mean 1000, to be used for approximating Poisson % numbers with large means (number of cells > 25/u) in the simulation, % using the standard normal distribution gen_poissMat = poissrnd(1000, 12, runTime*MC); % generate Poisson numbers % with mean 1000 poissMat = (gen_poissMat - 1000)./sqrt(1000); % converting gen_poissMat % into a standard normal distribution for v = 1:MC % for each experiment in the Monte Carlo simulation Nstart = N_init; % re-initialize starting populuation Ncompiled = zeros(length(t),4); % matrix storing the growth dynamics of % the 4 subpopulations over time (column 1: SS, 2: RS, 3: SR, 4: RR) Ncompiled(1,:) = Nstart; dN = zeros(length(Nstart),length(Nstart)); % matrix containing the % number of cells transitioning between phenotype states for i = 1:runTime if Nstart(1) <= 25/u % Transformations of SS dN(1,2) = poissrnd(u*Nstart(1)); dN(1,3) = poissrnd(u*Nstart(1)); % Proliferation and death of SS cells Nstart(1) = poissrnd(Nstart(1)*k_i(1)); else % Transformations of SS dN(1,2) = round(poissMat(1,(v-1)*runTime+i)*sqrt(u*Nstart(1)) ... + u*Nstart(1)); dN(1,3) = round(poissMat(2,(v-1)*runTime+i)*sqrt(u*Nstart(1)) ... + u*Nstart(1)); % Proliferation and death of SS cells Nstart(1) = round(poissMat(3,(v-1)*runTime+i)*sqrt(Nstart(1)*k_i(1)) ... + Nstart(1)*k_i(1)); end if Nstart(2) <= 25/u % Transformations of RS dN(2,1) = poissrnd(u*Nstart(2)); dN(2,4) = poissrnd(u*Nstart(2)); % Proliferation and death of RS cells Nstart(2) = poissrnd(Nstart(2)*k_i(2)); else % Transformations of RS dN(2,1) = round(poissMat(4,(v-1)*runTime+i)*sqrt(u*Nstart(2)) ... + u*Nstart(2)); dN(2,4) = round(poissMat(5,(v-1)*runTime+i)*sqrt(u*Nstart(2)) ... + u*Nstart(2)); % Proliferation and death of RS cells Nstart(2) = round(poissMat(6,(v-1)*runTime+i)*sqrt(Nstart(2)*k_i(2)) ... + Nstart(2)*k_i(2)); end if Nstart(3) <= 25/u % Transformations of SR dN(3,1) = poissrnd(u*Nstart(3)); dN(3,4) = poissrnd(u*Nstart(3)); % Proliferation and death of SR cells Nstart(3) = poissrnd(Nstart(3)*k_i(3)); else % Transformations of SR dN(3,1) = round(poissMat(7,(v-1)*runTime+i)*sqrt(u*Nstart(3)) ... + u*Nstart(3)); dN(3,4) = round(poissMat(8,(v-1)*runTime+i)*sqrt(u*Nstart(3)) ... + u*Nstart(3)); % Proliferation and death of SR cells Nstart(3) = round(poissMat(9,(v-1)*runTime+i)*sqrt(Nstart(3)*k_i(3))... + Nstart(3)*k_i(3)); end if Nstart(4) <= 25/u % Transformations of RR dN(4,2) = poissrnd(u*Nstart(4)); dN(4,3) = poissrnd(u*Nstart(4)); % Proliferation and death of RR cells Nstart(4) = poissrnd(Nstart(4)*k_i(4)); else % Transformations of RR dN(4,2) = round(poissMat(10,(v-1)*runTime+i)*sqrt(u*Nstart(4)) ... + u*Nstart(4)); dN(4,3) = round(poissMat(11,(v-1)*runTime+i)*sqrt(u*Nstart(4)) ... + u*Nstart(4)); % Proliferation and death of RR cells Nstart(4) = round(poissMat(12,(v-1)*runTime+i)*sqrt(Nstart(4)*k_i(4)) ... + Nstart(4)*k_i(4)); end % Adding and substracting transforming cells for h = 1:length(Nstart) Nstart(h) = Nstart(h) - sum(dN(h,:)) + sum(dN(:,h)); end Ncompiled(i+1,:) = Nstart; % start the next iteration with the current number of cells end SS(:,v) = Ncompiled(:,1); RS(:,v) = Ncompiled(:,2); SR(:,v) = Ncompiled(:,3); RR(:,v) = Ncompiled(:,4); total(:,v) = SS(:,v)+RS(:,v)+SR(:,v)+RR(:,v); % Find the time when RR cells first appear if sum(find(RR(:,v) >= 1)) ~= 0 T_RR(v) = find(RR(:,v) >= 1,1,'first')-1; end % Find the time when the tumor reaches a lethal tumor burden of 10^12 % cells if sum(find(total(:,v) >= 10^12)) ~= 0 T_lethal(v) = find(total(:,v) >= 10^12,1,'first')-1; end end % -------- % Calculate mean and standard deviation of T_lethal and T_RR % T_lethal validExp_ind_lethal = isnan(T_lethal) == 0; % if T_lethal = NaN (because % T_lethal > runTime), exclude the experiment nValidExp_lethal = sum(validExp_ind_lethal); % number of valid experiments T_lethal_valid = T_lethal(validExp_ind_lethal); % valid T_lethal data points T_lethal_mean = mean(T_lethal_valid); % average T_lethal T_lethal_std = std(T_lethal_valid); % standard deviation of T_lethal % T_RR validExp_ind_RR = isnan(T_RR) == 0; % exclude experiments for which T_RR is % not valid nValidExp_RR = sum(validExp_ind_RR); % number of valid experiments T_RR_valid = T_RR(validExp_ind_RR); % valid T_RR data points T_RR_mean = mean(T_RR_valid); % average T_RR T_RR_std = std(T_RR_valid); % standard deviation of T_RR T_RR_med = median(T_RR_valid); fprintf('T_lethal (simulated) = %.2f +- %.2f, (n = %.2f)\n', ... T_lethal_mean, T_lethal_std, nValidExp_lethal); fprintf('T_RR (simulated) = %.2f +- %.2f, (n = %.2f)\n', ... T_RR_mean, T_RR_std, nValidExp_RR); % -------- % Plot growth dynamics % Plot the average dynamics of 10% of experiments around the median T_RR_plot = find(T_RR >= prctile(T_RR_valid,45) & ... T_RR <= prctile(T_RR_valid,55)); SS_plot = SS(:, T_RR_plot); RS_plot = RS(:, T_RR_plot); SR_plot = SR(:, T_RR_plot); RR_plot = RR(:, T_RR_plot); total_plot = total(:, T_RR_plot); SS_ave = mean(SS_plot,2); RS_ave = mean(RS_plot,2); SR_ave = mean(SR_plot,2); RR_ave = mean(RR_plot,2); total_ave = mean(total_plot,2); % Expressing the average number of cells as logistic numbers lgSS = log10(SS_ave); lgRS = log10(RS_ave); lgSR = log10(SR_ave); lgRR = log10(RR_ave); lgT = log10(total_ave); % Plot only data points where total number of cells < 10^12 if max(lgT > 12) max_ind = find(lgT>=12,1,'first'); else max_ind = length(t); end t_plot = t(1:max_ind); lgSS = lgSS(1:max_ind); lgRS = lgRS(1:max_ind); lgSR = lgSR(1:max_ind); lgRR = lgRR(1:max_ind); lgT = lgT(1:max_ind); figure(1); clf; hold on; pl = plot(t_plot,lgSS,'b',t_plot,lgRS,'m',t_plot,lgSR,'r',... t_plot,lgRR,'c','LineWidth',1.5); pl2 = plot(t_plot,lgT,'LineWidth',1.5); plot([0 runTime],[9 9],'k--',[t_plot(end) t_plot(end)],[0 12],'k--',... 'LineWidth',1); axis([0 runTime 0 12]); title(strcat('Simulation model, T_R_R= ',num2str(T_RR_mean),... ', no. of experiments= ',num2str(nValidExp_RR))); xlabel('Time (generation)'); ylabel('log Cell number') legend('SS','RS','SR','RR','Total','Location','SouthOutside'); c = [1 0.8 0]; % setting the color of the 'Total' curve pl2.Color = c; %% =========================================================================== % Part 2: Comparing the calculated T_RR from the stochastic model and the % theoretical approximation of T_RR N_SS_0 = 10^9; % initial number of SS cells at t = 0 probDist_T_RR = zeros(runTime+1,1); % P(T_RR = t); probability distribution % curve of T_RR % Calculating the probability that T_RR = t for each time point in the % range observed for i = 1:runTime r = 2*u^2/(k_i(2) - k_i(1))*N_SS_0; probDist_T_RR(i+1) = exp(-r*((1-k_i(2)^(i-1))/(1-k_i(2))... - (1-k_i(1)^(i-1))/(1-k_i(1))))*(1-exp(-r*(k_i(2)^(i-1)-k_i(1)^(i-1)))); end % Calculate the mean T_RR computed by the theoretical model T_RR_mean_the = sum(probDist_T_RR.*t'); % mean T_RR (theoretical) percentDiscrepancy_T_RR = abs(T_RR_mean_the-T_RR_mean)/T_RR_mean*100; % percent discrepancy of theoretical mean T_RR compared to simulated mean % T_RR fprintf('T_RR (theoretical) = %.2f\n', T_RR_mean_the); fprintf('Percent discrepancy of theoretical mean T_RR from the simulation = %.2f percent\n',... percentDiscrepancy_T_RR); dist_T_RR_the = round(probDist_T_RR*MC); T_RR_the_list = zeros(sum(dist_T_RR_the),1); ind_counter = 0; T_RR_the_list(1:dist_T_RR_the(1)) = t(1); ind_counter = ind_counter + dist_T_RR_the(1); for h = 2:length(dist_T_RR_the) T_RR_the_list(ind_counter+(1:dist_T_RR_the(h))) = t(h); ind_counter = ind_counter + dist_T_RR_the(h); end T_RR_med_the = median(T_RR_the_list); percentDiscrepancy_T_RR_med = abs(T_RR_med_the-T_RR_med)/T_RR_med*100; fprintf('median T_RR (simulated) = %.2f\n', T_RR_med); fprintf('median T_RR (theoretical) = %.2f\n', T_RR_med_the); fprintf('Percent discrepancy of theoretical median T_RR from the simulation = %.2f percent\n',... percentDiscrepancy_T_RR_med); % -------- % Plot the histogram of T_RR (simulated) and the probability distribution % curve of T_RR (theoretical) figure(2); clf; hold on; h = histogram(T_RR_valid,round(nValidExp_RR/3)); plot(t,probDist_T_RR, 'k--', 'LineWidth',1.5); xlim([0 200]);xlabel('T_R_R (generation)'); ylabel('Probability of RR emerging'); title('Probability distribution of T_R_R'); legend('Simulation','Approximation'); h.Normalization = 'probability'; h.FaceColor = 'k'; h.EdgeColor = 'k';