%%
%CHE 216 - Enzyme Kinetics - Michaelis-Menten
%

% Introduction to M-M equation
% v = v_m [S]/(k_s + [S])

%vm is the highest possible rate, when all the enzyme 
%is binding the substrate all the time
%
%k_s is the saturation constant, the 
%concentration at half maximum velocity
%

figure(1);
vm = 1;
ks = 1;
Sall = linspace(0,5,1000);
vall = vm*Sall./(ks+Sall);
plot(Sall,vall,'-');
xlabel('S (substrate concentration)')
ylabel('v (rate of product)')
title(['Michaelis-Menten kinetics   ' ...
    'vm = ' num2str(vm),', ks = ' num2str(ks)])

%% Let's look at some data
S = [1.3 ,  1.8,    3,   4.5,     6,    8,    9];
v = [0.1, 0.13, 0.22, 0.275, 0.335, 0.35, 0.36];
%record the number of data points
sizedata = size(S);
numdata = sizedata(2);

figure(2);
plot(S,v,'sr');
xlabel('S (substrate concentration)')
ylabel('v (rate of product)')
title('the raw data');

%% transform to linear, and plot
%y_i = a    +     b * x_i
%1/v = 1/vm + ks/vm * 1/[S];

oS = S.^(-1);
ov = v.^(-1);

figure(3);
plot(oS,ov,'sr')

%% three ways of setting up the equations
% first, lets look at the sums
one = ones(numdata,1);

Sn  = sum(one)
SX  = sum(oS)
SXX = sum(oS.^2)
SY  = sum(ov)
SXY = sum(oS.*ov)

%% Now let's look at the normal equations

Z = [one oS']

% let's see what Z'*Z is . . . 
Z'*Z
Z'*ov'

% It's giving us the same equations!  That's what Z'Z is -- it's just
% bookkeeping

% Z'*Z A = Z'Y' 
% A = inv(Z'*Z) * Z'Y'
% or A = Z'*Z \ Z'Y' 

A = Z'*Z \ Z'*ov';

vm = 1/A(1)
ks = A(2)*vm

%% Now, use matlab tools

[p, s] = polyfit(oS,ov,1)

vm_poly = 1/p(2)
ks_poly = p(1)*vm_poly

%% now plot 
Sline = linspace(0,9,100);

vline = (vm.*Sline) ./ (ks +Sline);

% compute r2 of the linear fit

St = var(ov)
ovmod = 1/vm + ks/vm*oS;
Sr = sum((ovmod-ov).^2)
r2 = (St-Sr)/St

numpoints = size(S);
s = sqrt(Sr/(numdata-2))

%error analysis 

M = inv(Z'*Z);
ErrA(1) = sqrt(M(1,1)*s^2);
ErrA(2) = sqrt(M(2,2)*s^2);

%now we need to propagate the error -- vm = 1 /A(1), so the
%error in vm = 1/(vm)^2 * ErrA(1).  ks = A(2) * vm, so the error in ks is
%dks/ks = dA(2)/A(2) + dvm/vm 

MErrA(1) = ErrA(1)*vm^-2;
MErrA(2) = ks * (ErrA(2)/ks + MErrA(1)/vm);

% compute the sum of the residuals for the untransformed curve

vmod = vm*S./(ks+S);
Sr = sum((v-vmod).^2)

figure(4);
hold on;
plot(Sline,vline);
plot(S,v,'sr');
xlabel('S (substrate concentration)')
ylabel('v (rate of product)')
title(['Linear fit; Sr = ' num2str(Sr) ...
    ', vm = ' num2str(vm) ' +/- ' num2str(MErrA(1)) ...
    ', ks = ' num2str(ks) ' +/- ' num2str(MErrA(2))])

hold off;

%% Now, nonlinear fits

% compute the derivatives with respect to the parameters
%
% dv/dvm = S/(ks + S)
%
% dv/dks = -vm*S*(ks+S)^-2
%
% vm = a(1)
% ks = a(2)

maxsteps = 100;
a = [1 8];
tol = 10^-6;

for i=1:maxsteps
% rate at the given values 
   fv = a(1)*S./(a(2)+S);
% compute gradient with respect to the parameters
   dvdvm = S./(a(2)+S);
   dvdks = -(a(1)*S).*(a(2)+S).^-2;
% compute the normal equation matrix -- here' it consists of the derivative 
% evaluated at each data point, not the function.
   Z = [ dvdvm' dvdks'];

% Solve for da, the increments
% AX = B solved by X = A \ B; could do this by inverse as well.
   da = Z'*Z \ Z'*(v'-fv');

   err = norm(da)
   a = a + da'
   if (err < tol) 
       break   
   end
end

vm = a(1)
ks = a(2)

%% Now plot the nonlinear fit
Sline = linspace(0,9,100);
vline = (vm.*Sline) ./ (ks +Sline);

% compute the sum of the residuals for the curve from the nonlinear fit

vmod = vm*S./(ks+S);
Sr = sum((v-vmod).^2)

figure(5);
hold on;
plot(Sline,vline);
plot(S,v,'sr');
xlabel('S (substrate concentration)')
ylabel('v (rate of product)')
title(['Nonlinear fit; Sr = ' num2str(Sr) ...
    ', vm = ' num2str(vm),', ks = ' num2str(ks)])
hold off;
