%
%  hefib_ampanmda.m   Steve Cox, 2/1/09
%
% solve the active sealed fiber problem with an AMPA/NMDA synapse on a spine
%
% usage:   hefib_ampanmda(a,ell,N,dt,Tfin,Ns,pinc)
%
%		a = fiber radius (cm)
%		ell = fiber length (cm)
%       N = number of compartments
%       dt = timestep (ms)
%       Tfin = end time (ms)
%       Ns = spine compartment number(s)
%       pinc = number of time steps between plot3 calls
%
% example:	hefib_ampanmda(2e-4,.5,400,.02,14.9,120,10)
%

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

close all
Cm = 1;		% micro F / cm^2
R2 = 0.1;		% 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

VL = -60;       % leak potential  mV
GL = 0.3;      % mS/(cm)^2
tau = Cm/GL;    % membrane time constant
lam2 = a/2/R2/GL;   % square of space constant

% spine and synapse parameters

ass = 5e-5;     % cm
ellss = 1e-4;   % cm
Ash = 5e-8;     % cm
gss = pi*ass^2/R2/ellss;
km = 2*0.19;      % off rate, 1/ms
kp = 1.1;       % on rate, 1/(mM ms)
T0 = 1;          % transmitter conc, mM
r0 = T0/(T0+km/kp);  % effective fraction of open receptors
t1 = 4;         % start neurotran release
t2 = 5;         % stop neurotran release
gbarsyn = 0*1.3*9e4/55;  % mS/cm^2
Vsyn = 0;       % synaptic reversal potential
r = 0;          % time course of synaptic conductance
T1 = 0;
k1 = kp*T1+km;
rN = 0;
kmN = 4*0.0066;
kpN = 0.072;
KdN = kmN/kpN;
r0N = T0/(T0+KdN);
gbarsynN = 2e2;
k1N = kpN*T1+kmN;

e = ones(N,1);
v = VL*e;		% initial conditions
w = VL;

vK = -72;        % mV
vNa = 55;      % mV
vCl = -60;        % mV
VL = vCl;

GK = 40 - 20*(x'>0.05).*(x'<0.1);        % mS/(cm)^2
GNa = 44 + 560*(x'>0.05).*(x'<0.1);      % mS/(cm)^2
GL = 0.3;      % mS/(cm)^2
GK = GK/GL;
GNa = GNa/GL;
tau = Cm/GL;
lam2 = a/2/R2/GL;

Vr = -62.4;
e = ones(N,1);
v = Vr*e;		% initial conditions
w = v(1);
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)));

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;
plot3(x,t*ones(N,1),v)
hold on

Nt = 1+round(Tfin/dt);
tim = zeros(Nt,1);
GsynA = tim;
GsynN = tim;
W = tim;
W(1) = VL;
Iampa = tim;
Inmda = tim;

while (t <= Tfin)
    
      t = t + dt;
      tm = mod(t,20);
      T2 = T0*(tm>t1)*(tm<t2);
      k2 = kp*T2+km;
      k2N = kpN*T2+kmN;
      
      r = ((2/dt-k1)*r + kp*(T1+T2))/((2/dt)+k2);
      gsynA = gbarsyn*r0*r;
      
      rN = ((2/dt-k1N)*rN + kpN*(T1+T2))/((2/dt)+k2N);
      gsynN = gbarsynN*r0N*rN;
      gsynN = gsynN./(1+2*exp(-w*0.062)/3.57);
      
      gsyn = gsynA + gsynN;
      
      T1 = T2;
      k1 = k2;
      k1N = k2N;
        
      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)) );

      cNa = GNa .* m.^3 .* h;
      cK = GK .* n.^4;

      ed = 1 + dt*(1 + cNa + cK)/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 + vCl)/tau;

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

      v = B \ rhs;
      
      tcnt = tcnt + 1;
      tim(tcnt) = t;
      GsynA(tcnt) = gsynA;
      GsynN(tcnt) = gsynN;
      W(tcnt) = w;
      Iampa(tcnt) = gsynA*Ash*(w-Vsyn);
      Inmda(tcnt) = gsynN*Ash*(w-Vsyn);

      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)
hold off

figure(2)
subplot(3,1,1)
plot(tim,GsynA,'linewidth',2)
hold on
plot(tim,100*GsynN,'--','linewidth',2)
ylabel('(mS/cm^2)','fontsize',14)
legend('g_{ampa}','100g_{nmda}','location','best')

subplot(3,1,2)
plot(tim,W,'linewidth',2)
ylim([-70 0])
ylabel('w (mV)','fontsize',14)

subplot(3,1,3)
plot(tim,Iampa*1e3,'linewidth',2)
hold on
plot(tim,Inmda*1e6,'--','linewidth',2)
legend('I_{ampa}','1000I_{nmda}','location','best')
xlabel('t  (ms)','fontsize',14)
ylabel('(nA)','fontsize',14)

figure(3)
plot(x,v,'linewidth',2)
xlabel('x  (cm)','fontsize',16)
ylabel('v  (mV)','fontsize',16)

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