clear;
close;

n=100; % number of discretized spatial points
L=1;   % spatial length
dx=L/n; % incremental change in x coordinate
dt=0.05;   % incremental change in time
T=1000;   % total integration time
lam0=0.1; % constant background level of ligand
lamA=0.9; % amplitude of ligand peak
tplot=1; % plot frequency
A=3;   % magnitude of peak at trailing edge (choose 0 for uniform distribution)
%rom Table 1 (Narang et al.)

kappaeps=0.02;
kappa_f=1.1;
kappa_i=0.1;
kappa_p=0.1;
eps_d=0.1;
psi_p=0.01;
psi_i=0.01;
delta_p=0.0001;
delta_r=0.0001;

% initial conditions

pold=0.0991*ones(1,n);
iold=0.0991;

sigma=0.5*n;
for x=1:n,
    lam(x)=lam0*(1+lamA*exp(-(x-n/2)^2/sigma));
    lamold(x)=lam0*(1+lamA*exp(-(x-n/2)^2/sigma));
end;
rhoold=(1+1/lam0+eps_d)./(1+1./lam+eps_d);

% plot initial conditions

xvec=linspace(0,1,n);
subplot(2,2,1);
plot(xvec,lam,'bo');
title('cAMP GRADIENT');
subplot(2,2,2);
plot(xvec,rhoold,'go');
axis([0 1 0 2]);
title('ACTIVE RECEPTORS');
subplot(2,2,3);
plot(xvec,pold,'ro');
title('PIP2s (ACTIVATOR)');
axis([0 1 0 1]);
subplot(2,2,4);
ivec=iold*ones(n,1);
plot(xvec,ivec,'yo');
title('IP3 (INHIBITOR)');
axis([0 1 0 0.2]);
drawnow;
        
input('Hit a key to start simulation');
close;
figure;


% start of iterative time loop

for timestep=1:(T/dt);
    
    t=timestep*dt;

    if t==50
        input('Hit a key to apply ligand distribution at trailing edge');
        close;
        figure;
        for x=1:n/2,
            lam(x)=lam0*(1+A*exp(-(x)^2/sigma));
            lamnew(x)=lam(x);
        end;
        for x=n/2+1:n,
            lam(x)=lam0*(1+A*exp(-(x-n)^2/sigma));
            lamnew(x)=lam(x);
        end;
        
        rhonew=rhoold.*(1+1./lamold+eps_d)./(1+1./lamnew+eps_d);
        rhoold=rhonew;
    end;
    
    % equations for x=1, incorporating periodic boundary conditions 
    drhodt(1)=(1-rhoold(1))*kappaeps/(1+1/lam(1)+eps_d) + ...
        delta_r*(rhoold(n)-2*rhoold(1)+rhoold(2))/dx^2;
    dpdt(1)=kappa_f*rhoold(1)*pold(1)^2*(1-pold(1))-pold(1)*iold + ...
        psi_p-kappa_p*pold(1)+ delta_p*(pold(n)-2*pold(1)+pold(2))/dx^2;
    
    % equation for internal points (from 2 to n-1)
    for x=2:n-1,
        drhodt(x)=(1-rhoold(x))*kappaeps/(1+1/lam(x)+eps_d) + ...
            delta_r*(rhoold(x-1)-2*rhoold(x)+rhoold(x+1))/dx^2;
        dpdt(x)=kappa_f*rhoold(x)*pold(x)^2*(1-pold(x))-pold(x)*iold + ...
            psi_p-kappa_p*pold(x)+ delta_p*(pold(x-1)-2*pold(x)+pold(x+1))/dx^2;
    end;

    % equations for x=n, incorporating periodic boundary conditions 
    drhodt(n)=(1-rhoold(1))*kappaeps/(1+1/lam(1)+eps_d) + ...
        delta_r*(rhoold(n-1)-2*rhoold(n)+rhoold(1))/dx^2;
    dpdt(n)=kappa_f*rhoold(n)*pold(n)^2*(1-pold(n))-pold(n)*iold + ...
        psi_p-kappa_p*pold(n-1)+ delta_p*(pold(n-1)-2*pold(n)+pold(1))/dx^2;

    % calculate di/dt
    didt=kappa_f*mean(rhoold)*mean(pold)^2*mean(1-pold)-mean(pold)*iold + psi_i-kappa_i*iold;

    % update the new values of rho, p, and i.
    for x=1:n,
        rhonew(x,timestep)=rhoold(x)+drhodt(x)*dt;
        pnew(x,timestep)=pold(x)+dpdt(x)*dt;
        rhoold(x)=rhoold(x)+drhodt(x)*dt;
        pold(x)=pold(x)+dpdt(x)*dt;
    end;
    inew(timestep)=iold+didt*dt;
    iold=iold+didt*dt;
    
    if rem(t,tplot)==0
        xvec=linspace(0,1,n);
        subplot(2,2,1);
        plot(xvec,lam,'bo');
        title('cAMP GRADIENT');
        subplot(2,2,2);
        plot(xvec,rhoold,'go');
        title('ACTIVE RECEPTORS');
        axis([0 1 0 2]);
        subplot(2,2,3);
        plot(xvec,pold,'ro');
        title('PIP2s (ACTIVATOR)');
        axis([0 1 0 1]);
        subplot(2,2,4);
        ivec=iold*ones(n,1);
        plot(xvec,ivec,'yo');
        title('IP3 (INHIBITOR)');
        axis([0 1 0 0.2]);
        drawnow;
        t
        iold
    end;   
end;