function BE20109F13_singlefit
% 20.109 Matlab tutorial: nonlinear curve fitting
% Thanks to Anthony Soltis (BE PhD student) for sharing some code to build from.

clc; clear all; close all
%% Curve fit using 'nlinfit' function

% TO CALL nlinfit:
% [param_estimate,residuals] = nlinfit(xdata,ydata,fitfunc,initial_guesses)
% nlinfit is a non-linear curve fit algorithm that arrives at the best
% answer by minimizing the distance between the experimental data and the
% output from the model (a least-squares method).
%
% INPUTS:
%   xdata - Erlotinib Concentration in uM
%   ydata - Densitometry or viability data matrix named 'data'
%   fitfunc - function definining relationship between xdata and ydata;
%             called with function handle (use @ symbol)
%   initial_guesses - vector of initial guesses for parameter fits
%
% FUNCTION USED for fitfunc:
% Y = min+(max-min)/(1+10^(X-LogIC50))
% This is a modified 4 parameter logistic fit equation. It models the data
% as a sigmoidal curve spanning from your minimum to maximum values. We are
% assuming a hill coefficient slope = 1 because the number of data points
% that we have to capture the change in slope is relatively small (only one
% or two points will fall in the transition zone).
%
% OUTPUTS
%   param_estimate - vector of final parameter estimates
%   Jacobian - this matrix will be used when determining the 95% confidence
%   intervals of our prediction
%   residuals - vector of residuals between data points and fit function
%               using fit parameters. 


% Ok -- let's get started:
% Assign data for fitting
d = importdata('data.xls');
text = d.textdata;
data = d.data';


%Unlike the singlefit file, this code will import your averaged data and
%loop through all six of the InhibitorX concentrations to find a IC50 value
for i = 1:length(data(:,1))
y(i,:) = data(i,:); 
Erlotinib_vals = [-4, -3, -2, -1, 0, 1];
%These values span the range of Erlotinib concentration that was used. 
%To avoid errors associated with handling zeros, we are using
%0.0001 as our bottom value instead of zero

maxy = max(y(i,:)); %Largest data point for an initial 'top of curve' guess
miny = min(y(i,:)); %Smallest data point for an initial 'bottom of curve' guess

% Call the fitting function
initial_guesses = [-1 miny maxy]; % Provide initial estimate(s) for parameter(s); here we've started with 0.1 uM and our largest/smallest values
[params_nlinfit,residuals,Jacobian] = nlinfit(Erlotinib_vals,y(i,:),@fitfunc,initial_guesses); %Passing variables to our function

%Calculate the 95% confidence intervals
%For more information about how nlparci calculates the 95% confidence
%intervals see: http://www.mathworks.com/support/solutions/en/data/1-194W2/index.html?product=ML&solution=1-194W2
IC50_confint = nlparci(params_nlinfit,residuals,'jacobian',Jacobian);

%Convert IC50 values back to uM from the Log value
IC50ci(i,:) = 10.^IC50_confint(1,:);
IC50(i,1) = 10.^params_nlinfit(1);


% Check fit visually
y_check = fitfunc(params_nlinfit,Erlotinib_vals); % Call fitting function with fit param
figure()
subplot(1,2,1)
plot(Erlotinib_vals,y(i,:),'ro',Erlotinib_vals,y_check,'b')
title('Non-linear fit of Inhibitor Dose-Response')
xlabel('Log[Erlotinib] (uM)')
ylabel([text(i)])
legend('Data','Fit',2)
subplot(1,2,2)
plot(1:length(y_check),residuals,'o',...
     1:length(y_check),zeros(size(Erlotinib_vals)),'k--')
title('Residual Plot')
end


% Put it all together in an output matrix where column 1 is IC50, col 2 is
% lower 95% CI and col 3 is upper 95% CI
output(:,1) = IC50;
output(:,2) = IC50ci(:,1);
output(:,3) = IC50ci(:,2);
display('The model_out exel file contains IC50 (col 1) and 95% confidence interval of fit in uM')
xlswrite('model_out',output)


%% Fitting function 
function yhat = fitfunc(k,x)
% INPUTS:
%   k - vector of parameters for fitting:
%   x - independent variable data (drug concentration!)
% 
% OUTPUTS:
%   yhat - value of function at points 'x' given parameters 'k' 

LogIC50 = k(1); 
miny = k(2);
maxy = k(3);
yhat = miny+((maxy-miny)./(1+10.^(x-LogIC50)));