ToxFit.m 4.22 KB
Newer Older
Gupta's avatar
Gupta committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
function ToxFit(input_file,iterations,start_time)
% ToxFit calculates IC50 and growth parameters for a given set of
% measurements.
% Usage: [IC50,growth] = ToxFit(input_file,iterations,start_time)
% input:
%   input_file: xlsx file giving time points in the first column and test 
%               concentrations in the first row.
%   iterations: gives the number of bootstrap iterations
%   start_time: gives the time point at which the fitting should start
% output:
%   IC50: gives the compound's concentration needed to half the cell growth
%   growth: gives the estimated growth rate of the cells during the experiment

global test_concs increment_start_time
RandStream.setGlobalStream(RandStream('mt19937ar','seed',12345));

% loading input data
Data = xlsread(input_file);

% test concentations must be given in the first row of the input file
test_concs = Data(1,2:end);
M = length(test_concs);

% time points for all measurements must be given in the first column of the input file
time_points = Data(2:end,1);
index_timept = find(time_points>start_time,1);
increment_start_time = time_points(index_timept);

% X_cell_data is a matrix holding all measured data. This matrix must have
% a size of length(time_points)*length(test_concs)
X_cell_data = Data(2:end,2:end);


%% Bootstrapping and Data fitting

x_boot = zeros(iterations,M+2);

for iter = 1:iterations
    
    % preparing x and y vectors for data fitting
    xdata = [];
    ydata = [];
    for i=1:M
        xnew = time_points(index_timept:end);
        ynew = (X_cell_data(index_timept:end,i));        
        % generate a random sample of rows (with replacement)
        samp = ceil(rand(length(ynew),1) .* length(ynew));        
        % sample rows from x and y data
        xdata = [xdata xnew(samp,:)];
        ydata = [ydata ynew(samp,:)];
    end
    
    % preparing random initial values for the parameters: IC50, growth rate, I_m (for m=1:M)
    x0 = [(5000)*rand(1,1), rand(1,1), 100*rand(1,M)];
    
    % defining lower and upper bound for the parameter vector
    lb = zeros(size(x0));
    ub = inf*ones(size(x0));
    
    % preparing fitting options
    options = optimset('FunValCheck','on','MaxFunEvals',20000); 
    
    % lsqcurvefit --> Fit a model on the sample in order to estimate the parameters obtained from this bootstrap iterations. The parameters obtained in each iteration are put in a row of x_boot
    [x_boot(iter,:),~,~,~] = lsqcurvefit(@kinetic_function,x0,xdata,ydata,lb,ub,options);
end

% Compute 95% confidence interval for the estimated paramters
x_median = median(x_boot,1);
x_95CI = prctile(x_boot,[5 95]);

% final return values
IC50 = x_median(1);
growth = x_median(2);

%% Ploting results

figure(), clf, hold on

% plotting experimental data
cmap = colormap(jet); cmap=cmap(20:end,:);
c = cmap(round(linspace(0,1,M)*(size(cmap,1)-1)+1),:);
L = {};
for i=1:M
    plot(time_points,X_cell_data(:,i),'.','MarkerSize',10,'Color',c(i,:));
    L{i} =  sprintf('%0.5g \\muM',test_concs(i));
end
legend(L,'Location','NorthWest');
title(input_file);
xlabel('Time (h)','FontSize',12); ylabel('Normalized cell index','FontSize',12);set(gca,'FontSize',12);

% plotting estimated data
t_plot = [25:.1:100]'*ones(1,M);
plot(t_plot,kinetic_function(x_median,t_plot),'c--','LineWidth',2)
t_plot = time_points(index_timept:end)*ones(1,M);
plot(t_plot,kinetic_function(x_median,t_plot),'b-','LineWidth',2)


disp(['IC50=',num2str(IC50), ', Growth rate=',num2str(growth)])

% function Result = kinetic_function(x,x_data)
% global test_concs increment_start_time;
% 
% [n,m] = size(x_data);
% Result = zeros(n,m);
% for i = 1:m
% %     Result(:,i) = RHS_once(x,test_concs,i,0,x_data(:,i));
%     Result(:,i) = x(2+i)*exp((x_data(:,i)-increment_start_time)*x(2)*(x(1)/(x(1)+test_concs(i))));
%     
% end

function Result = kinetic_function(x,x_data)
% Result = kinetic_function(x,x_data)
% x:     parameter vector holding [IC50, growth rate, I_m] for m=1:M
% xdata: matrix containing time points; every row corresponds to a
%        different test concentation
global test_concs increment_start_time;

[n,m] = size(x_data);
Result = zeros(n,m);
for i = 1:m
%     Result(:,i) = RHS_once(x,test_concs,i,0,x_data(:,i));
    Result(:,i) = x(2+i)*exp((x_data(:,i)-increment_start_time)*x(2)*(x(1)/(x(1)+test_concs(i))));
    
end