function BE20109F13_fitforWB
% 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+(X-IC50))
% 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 point
% that we have to capture the change in clope 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. 



%Input your data below starting with the lowest concentration of inhibitor
%This should be the normalized densitometry data calculated from your WB 
%This example is ([pEGFR]/[EGFR])/[GAPDH]
y=[0.54,0.69,0.57,0.29,0.05,0.03];
Erlotinib_vals = (-4:1:1);%These values span the range of Erlotinib concentration that was used. 
%To avoid errors associated with handling zeros, we are instead using
%0.0001 as our bottom value
maxy = max(y); %Largest data point for an initial 'top of curve' guess
miny = min(y); %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,@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 the fitted IC50 and confidence intervals from Log values to uM
IC50ci = 10.^IC50_confint;
IC50 = 10.^params_nlinfit(1);

% Display result in command window
sprintf('Estimated IC50 from nlinfit is: %f uM', IC50)
sprintf('The lower 95 percent CI bound is: %f uM',(IC50ci(1,1)))
sprintf('The upper 95 percent CI bound is: %f uM',(IC50ci(1,2)))
%sprintf('Estinated Hill coefficient from nlinfit is: %f',params_nlinfit(2))

% 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,'ro',Erlotinib_vals,y_check,'b')
title('Non-linear fit of Inhibitor Dose-Response')
xlabel('Log [Erlotinib] (uM)')
ylabel('(pEGFR/EGFR)/GAPDH')
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')


%% 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)));