%Applicando un impulso di clamping (partendo da Vr), calcolare e disegnare l'andamento di: alpha, beta, n, m, h, costanti di tempo, 
%correnti ioniche e
% correnti capacitive 

%Impulso: potenzziale di membrana = VI per una durata di 20 msec
%condizioni: assone gigante di calamaro T=6.3 C

clear variables 
close all
clc

%% costanti

GNaM=120; %[mS/cm2]
GKM=36; %[mS/cm2]
GL=0.3; %[mS/cm2]
CM=1; %[uF/cm2]
Vr=-65; %[mv]
VNa=50; %[mv]
VK=-77; %[mv]
VL=-54.387; %[mv]



%% condizioni iniziali : le calcolo come alpha(V=Vr)/(alpha(V=Vr)+beta(V=Vr))

V0=Vr;
dv=0; %dv=Vm-Vr differenza rispetto al potenziale di riposo
an0=(0.01*(10-(dv)))/(exp((10-dv)/10)-1);
bn0=0.125/(exp(dv/80));
am0=(0.1*(25-(dv)))/(exp((25-dv)/10)-1); 
bm0=4/(exp(dv/18));
ah0=0.07/(exp(dv/20));
bh0=1/(exp((30-dv)/10)+1); 

n0=an0/(an0+bn0);
m0=am0/(am0+bm0);
h0=ah0/(ah0+bh0);


%% inizializzo le variabili con la giusta dimensione 

tmax=100; % [msec] Durata della simulazione
dt=0.001;% [ms] intervallo tra due passi della della simulazione 
steps=round(tmax/dt);
t=(0:steps-1)*dt;

  Vm = Vr*ones(steps,1);
  m = zeros(steps,1);
  h = zeros(steps,1);
  n = zeros(steps,1);
  GK=zeros(steps,1);
  GNa=zeros(steps,1);
  
 
   
%% impulso clamp 

 VI=20; %valore impulso Interessante a VNa
  durata=20; %durata impulso in msec
  ncampioni=round(durata/dt); %numero di passi o campioni relativi alla durata dell'impulso
  Vm(ncampioni:ncampioni+ncampioni)=VI;
  %figure
  %plot(t,Vm);

  
%%  setto le condizioni iniziali

  
  
  m(1)=m0;
  n(1)=n0;
  h(1)=h0;
  
%% Metodo Eulero: soluzione di Im=Ic+Ii -> impulso non propagante (space clamp) stimolato con una corrente esterna
%   

for i=1:steps-1
    
    
      
      %conduttanze
      GK(i)=GKM*(n(i)^4);
      GNa(i)=GNaM*(m(i)^3)*h(i);
      
      
      %correnti
      IK(i)=GK(i)*(Vm(i)-VK); %K
      INa(i)=GNa(i)*(Vm(i)-VNa); %Na
      IL(i)=GL*(Vm(i)-VL); %L
      Ii(i)=IK(i)+INa(i)+ IL(i);
            
      %Vm(i+1)=Vm(i)+dt*( I(i)- Ii(i))/CM;
      
      
      dv=Vm(i)-Vr; % cambia Vm -> cambiano i coefficienti alpha e beta
      
      an(i)=(0.01*(10-(dv)))/(exp((10-dv)/10)-1);
      bn(i)=0.125/(exp(dv/80));
      am(i)=(0.1*(25-(dv)))/(exp((25-dv)/10)-1); 
      bm(i)=4/(exp(dv/18));
      ah(i)=0.07/(exp(dv/20));
      bh(i)=1/(exp((30-dv)/10)+1); 
      
    

      % rquazioni differenziali coefficienti alpha e beta
      n(i+1)= n(i)+dt*(an(i)*(1-n(i))-bn(i)*n(i));
      m(i+1)= m(i)+dt*(am(i)*(1-m(i))-bm(i)*m(i));
      h(i+1)= h(i)+dt*(ah(i)*(1-h(i))-bh(i)*h(i));
      
     
    
      
end
  
 
%% Grafico i risultati


%N part
figure
subplot(2,1,1)
plot(t,Vm);
ylabel('[mV]');
xlabel('t [msec]');
legend('Vm');
ylim([-90 40]);
subplot(2,1,2)
plot(t,n);
hold on
plot(t(1:end-1),an);
plot(t(1:end-1),bn)
plot(t,n.^4);
legend('n','\alpha_n','\beta_n','n^4');
xlabel('t [msec]');

figure
plot(t(1:end-1),1./(an+bn));
legend('\tau_n');
xlabel('t [msec]');

%M,H art
figure
subplot(3,1,1)
plot(t,Vm);
ylabel('[mV]');
xlabel('t [msec]');
legend('Vm');
ylim([-90 40]);
subplot(3,1,2)
plot(t(1:end-1),am);
hold on
plot(t(1:end-1),bm)
legend('\alpha_m','\beta_m');
xlabel('t [msec]');
subplot(3,1,3)
plot(t(1:end-1),ah);
hold on
plot(t(1:end-1),bh)
legend('\alpha_h','\beta_h');
xlabel('t [msec]');

figure
subplot(3,1,1)
plot(t,Vm);
ylabel('[mV]');
xlabel('t [msec]');
legend('Vm');
ylim([-90 40])
subplot(3,1,2)
plot(t,m);
hold on
plot(t,m.^3);
plot(t,h);
legend('m','m^3','h');
xlabel('t [msec]');
subplot(3,1,3)
plot(t,(m.^3).*h);
legend('m^3h');
xlabel('t [msec]');

figure 
plot(t(1:end-1),Ii)
hold on
plot(t(1:end-1),IK)
plot(t(1:end-1),INa)
legend('Ii','Ik','INa');
%plot(t(1:end-1),CM*diff(Vm)/dt);


 
