function AmpVsPol
% This function simulates the steady state Levchenko and Iglesias model [1] 
% on ellipses with increasing polarity.
%
%
% Matt Onsum
% 
% Reference: 
% [1] Levchenko A., Iglesias, P.A. Models of Eukaryotic Gradient Sensing: 
% Applications to Chemotaxis of Amoebae and Neutrophils, Biophys J, 2002 
% 82: p 50-63 
%
%

close all
%As defined in text
polarity = [1:.1:2];

n_p = length(polarity);

Amp = zeros(n_p,1);

for ind_p = 1:n_p
   
    %Define ellipse with polarity(ind_p)
    %and put into polar coordinates
    p = polarity(ind_p);
    a = sqrt(2/(p^2+1));
    b = p*a;
    theta = linspace(-pi,pi);
    x = a*cos(theta);
    y = b*sin(theta);
    r = sqrt(x.^2 + y.^2);
    
    figure(1)
    hold on
    plot(x,y,'LineWidth',2);
    axis equal
    title('Set of Ellipses');
    xlabel('X','Fontsize',12);
    ylabel('Y','Fontsize',12);
    
    %Define Linear Field of attractant
    %L(x,y) = a + by;
    c = 1;
    d = .4;
  
    %Apply to ellipse
    G1 = d*r.*sin(theta)+c;
  
   % Solve LI Model with above input
    
    
   
   
   %Linear Amplifier-- Epsilon =0
    s0 = Iglesias02(r,c,d,0);
    
   %Nonlinear Amplifier-- Epsilon =.45
    s1 = Iglesias02(r,c,d,0.45);
    
    
 
    %Calculate Amp linear
    [P0,S] = polyfit(G1,s0,1);   
    alpha = P0(1);
    Amp0(ind_p)= alpha*min(G1)/(P0(1)*min(G1)+P0(2));
    
    %Calculate Amp, nonlinear
    [P1,S] = polyfit(G1,s1,1);   
    alpha = P1(1);
    Amp1(ind_p)= alpha*min(G1)/(P1(1)*min(G1) + P1(2));






end

 figure(2); 
 plot(polarity,Amp0/min(Amp0),'LineWidth',2);
 hold on
 plot(polarity,Amp1/min(Amp1),'g','LineWidth',2);
axis([1 2 0 2.5])
xlabel('Polarity','FontSize',12);
ylabel('Amplification', 'FontSize', 12);


function [r1] = Iglesias02(R,a,b,epsilon);
% Steady State Solution to L&I's model
% Reference

%1.0 Biochem Parameters-%
% 1.1 Activator Dynamics
ka  = 1;        %S-[+ka]-> A 
k_a = 1;        %S<-[-k_a]- A
KDa = ka/k_a;

% 1.2 Inhibitor Dynamics
ki = 0.05;      %S -[+ki]-> I
k_i = 0.05;     %S <-[-k_i]-I
KDi = ki/k_i;
D = 20;        %Diffusion of I

% 1.3 Lumped Parameters
alpha = 1;
beta = 1;
gamma = 1;

%2.0 Chemoattractant about ellipse
theta = linspace(-pi,pi,length(R));
s = a + b*R.*sin(theta);

%3.0 Steady State Values of Activator and Inhibitor in polar coordinates
A = KDa*(s);
I = KDi*a + ki*inv(k_i+D)*b*sin(theta);
F = A./I;

%4.0 distribution of response element r1*
r1 = alpha*gamma*F./(beta*(1-alpha*epsilon*F));