%
%  heafibAampa2nmda.m   Steve Cox, 2/1/07
%
% solve the active I_A fiber with AMPA & NMDA on spine
%
% usage:	heafibAampa2nmda(a,ell,N,dt,Tfin,Ns,pinc)
%
%		a = fiber radius (cm)
%		ell = fiber length (cm)
%             	N = number of compartments
%               dt = timestep (ms)
%		t1 = start of current pulse (ms)
%		t2 = end of current pulse (ms)
%               T = stopping time (ms)
%		Iapp = size of current pulse (micro amps)
%               pinc = number of time steps between plots
%                      if pinc > 0 then 3D full
%                      if pinc < 0 then 2D movie
%
% example:	heafibAampa2nmda(.0002,.5,200,.03,12,[50 150]',10)
%

function heafibampa2nmda(a,ell,N,dt,Tfin,Ns,pinc)

%close all
Cm = 1;		% micro F / cm^2
R2 = 0.034;		% k Ohm cm
dx = ell/N;		% patch length
A = 2*pi*a*dx;		% patch surface area
x = dx/2:dx:ell-dx/2;	% vector of patch midpoints

vK = -72;        % mV
vA = -75;
vNa = 55;      % mV
vCl = -47;        % mV
VL = vCl;

GK = 20;        % mS/(cm)^2
GNa = 120;      % mS/(cm)^2
GA = 47.7*(.2+2*x')/2.3;
GL = 0.3;      % mS/(cm)^2
GK = GK/GL;
GA = GA/GL;
GNa = GNa/GL;
tau = Cm/GL;
lam2 = a/2/R2/GL;

Vr = -61;
e = ones(N,1);
v = Vr*e;		% initial conditions
n = an(v(1))*e/(an(v(1))+bn(v(1)));
m = am(v(1))*e/(am(v(1))+bm(v(1)));
h = ah(v(1))*e/(ah(v(1))+bh(v(1)));
a = ainf(v(1))*e;
b = binf(v(1))*e;

B = spdiags([-e 2*e -e], -1:1, N, N)/dx/dx;
B(1,1) = 1/dx/dx;
B(N,N) = 1/dx/dx;
B = dt*lam2*B/tau;
dB = diag(B);

t = 0;
tcnt = 1;
subplot(2,2,1)
plot3(x,t*ones(N,1),v)
hold on
vcnt = 1;

% spine and synapse parameters

ass = 1e-5;
ellss = 1e-4;
Ash = 5e-8;
gss = pi*ass^2/R2/ellss;

t1 = 2;
t2 = 3;
T0 = 1;

kmA = 0.19;
kpA = 1.1;
KdA = kmA/kpA;
r0A = T0/(T0+KdA);
kA = kpA*T0+kmA;
gbarsynA = 3e3;

kmN = 0.0066;
kpN = 0.072;
KdN = kmN/kpN;
r0N = T0/(T0+KdN);
kN = kpN*T0+kmN;
gbarsynN = 3e2;

Vsyn = 0;

w = Vr*ones(size(Ns));

while (t <= Tfin)

      t = t + dt;

      gsynA = gbarsynA*(t>t1)*r0A*(exp(kA*min(0,t2-t))-exp(kA*(t1-t)));

      Iampa(:,tcnt) = Ash*gsynA*(w-Vsyn);

      gsynN = gbarsynN*(t>t1)*r0N*(exp(kN*min(0,t2-t))-exp(kN*(t1-t)));
      gsynN = gsynN./(1+2*exp(-w*0.062)/3.57);

      Inmda(:,tcnt) = Ash*gsynN.*(w-Vsyn);

      gsyn = gsynA + gsynN;

      bot = tau/dt + 1 + gsyn/GL + gss/Ash/GL;

      w = (tau*w/dt + gsyn*Vsyn/GL + gss*v(Ns)/GL/Ash + VL)./bot;

      ant = an(v);
      amt = am(v);
      aht = ah(v);
      n = ( n + dt*ant )./(1 + dt*(ant+bn(v)) );
      m = ( m + dt*amt )./(1 + dt*(amt+bm(v)) );
      h = ( h + dt*aht )./(1 + dt*(aht+bh(v)) );
      a = (a.*taua(v) + dt*ainf(v))./(dt + taua(v));
      b = (b.*taub(v) + dt*binf(v))./(dt + taub(v));

      cNa = GNa * m.^3 .* h;
      cK = GK * n.^4;
      cA = GA .* a.^3 .* b;

      ed = 1 + dt*(1 + cNa + cK + cA)/tau;
      ed(Ns) = ed(Ns) + dt*gss/GL/A/tau;

      B(1:N+1:end) = dB + ed;	      % update the diagonal

      rhs = v + dt*(cNa*vNa + cK*vK + cA*vA + vCl)/tau;

      rhs(Ns) = rhs(Ns) + dt*(gss*w/GL/A)/tau;

      v = B \ rhs;

      W(:,tcnt) = w;

      tcnt = tcnt + 1;
      
      if mod(tcnt,pinc) == 0
         plot3(x,t*ones(N,1),v)
      end

end

xlabel('x  (cm)','fontsize',16)
ylabel('t  (ms)','fontsize',16)
zlabel('v  (mV)','fontsize',16)
axis([0 ell 0 Tfin -80 60])

hold off

tim = linspace(0,Tfin,tcnt-1);

subplot(2,2,2)
plot(tim,1e3*Iampa(1,:))
hold on
plot(tim,1e3*Iampa(2,:),'--')
ax = axis;
ax(2) = 12;
axis(ax)
xlabel('t  (ms)','fontsize',16)
ylabel('I_{AMPA}  (nA)','fontsize',16)
hold off

subplot(2,2,3)
plot(tim,1e3*Inmda(1,:))
hold on
plot(tim,1e3*Inmda(2,:),'--')
ax = axis;
ax(2) = 12;
axis(ax)
xlabel('t  (ms)','fontsize',16)
ylabel('I_{NMDA}  (nA)','fontsize',16)
hold off

subplot(2,2,4)
plot(tim,W(1,:))
hold on
plot(tim,W(2,:),'--')
ax = axis;
ax(2) = 12;
axis(ax)
xlabel('t  (ms)','fontsize',16)
ylabel('w  (mV)','fontsize',16)
hold off

return

function val = an(v)
val = .02*(v+45.7)./(1-exp(-(v+45.7)/10));

function val = bn(v)
val = .25*exp(-0.0125*(v+55.7));

function val = am(v)
val = 0.38*(v+29.7)./(1-exp(-(v+29.7)/10));

function val = bm(v)
val = 15.2*exp(-0.0556*(v+54.7));

function val = ah(v)
val = 0.266*exp(-0.05*(v+48));

function val = bh(v)
val = 3.8./(1+exp(-(v+18)/10));

function val = ainf(v)
val = 0.0761*exp(0.0314*(v+94.22))./(1+exp(0.0346*(v+1.17)));
val = val.^(1/3);

function val = taua(v)
val = 0.3632 + 1.158./(1+exp(0.0497*(v+55.96)));

function val = binf(v)
val = 1./(1+exp(0.0688*(v+53.3)));
val = val.^4;

function val = taub(v)
val = 1.24 + 2.678./(1+exp(0.0624*(v+50)));