%SIMULAZIONE DEL POTENZIALE DI AZIONE 
% calcolare l'andamento del potenziale di membrana quando viene applicato
% un impulso di corrente di carica data in condizioni di space clamp (no
% propagazione lungo x)
%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))

%condizioni iniziale n_inf(V=Vr) h_inf(V=Vr)  m_inf(V=Vr) 
V0=Vr;
dv=0;

%utilizzo le formule ricavate da HH
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); %0.3177
m0=am0/(am0+bm0); %0.0529
h0=ah0/(ah0+bh0); %0.5961


%% variabili di integrazione

%NB: i coefficienti di vel di transizione sono in 1/msec!

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

%% inizializzo le variabili con la giusta dimensione 

  Vm = zeros(steps,1);
  m = zeros(steps,1);
  h = zeros(steps,1);
  n = zeros(steps,1);
  GK=zeros(steps,1);
  GNa=zeros(steps,1);
  IK=zeros(steps,1);   %corrente potassio
  INa=zeros(steps,1);  %corrente sodio 
  IL=zeros(steps,1);   %corrente leakege
  Ii=zeros(steps,1);   %corrente ionica 
  
  %ingresso forzante
  I=zeros(steps,1);    %corrente di membrana che impongo dall'esterno (ingresso forzante) 
  
 
%%  Corrente di stimolazione

%si fissa la "carica" dello stimolo, se poi fissiamo la durata otteniamo la
%corrente Is (carica integrale della corrente nel tempo) 
%nota: contanto intensità e durata (intensità corrente) x (durata) = carica
% sperimentando si ottiene una soglia sui  6.50 C/cm2 
  
stimolo=20e-9; %C/cm2  
%stimolo=6e-9; %C/cm2
%stimolo=7e-9; %C/cm2

durata=0.1; %msec
IS=(stimolo)/(durata*1e-3); %corrente in A/cm2
ncampioni=round(durata/dt); %numero di passi (campioni)  relativi alla durata dell'impulso
IS=IS*10^6;% [uA/cm2]
I(ncampioni:ncampioni+ncampioni)=IS;
% figure
% plot(t,I);
  
%% refrattario (calcoli fatti per stimolo di 20e-9)
 % N=30; %8-37 assoluto;  N=60 o 80  realativo il solito stimolo di prima non basta, uno pi? alto si (4Is) (Is) non ce la fa nel relativo h deve essere sufficientemente alto, vedi figura finale
 
 %N=30;
 %I(N*ncampioni:N*ncampioni+ncampioni)=4*IS; %REFRATTARIO: con  stimolo (Riperto lo stimolo mentre la membrana ? attivata)
  
  
  
%% firing periodico
%I(1:end)=10; %per firing periodico -> piccola corrente costante
  
%%  setto le condizioni iniziali

  
  Vm(1)=V0;
  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
      
      %#1 calcolo delle conduttanze
      GK(i)=GKM*(n(i)^4);
      GNa(i)=GNaM*(m(i)^3)*h(i);
    
      
      %#2 calcolo le correnti ioniche 
      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); %corrente ionica
      
      %#3 applico Eulero per aggiornare il potenziale di membrana al passo
      %successivo 
            
      Vm(i+1)=Vm(i)+dt*( I(i)- Ii(i))/CM;
      
      %#4 calcolare alpha e beta 
      
      dv=Vm(i)-Vr; % V'
      
      an=(0.01*(10-(dv)))/(exp((10-dv)/10)-1);
      bn=0.125/(exp(dv/80));
      am=(0.1*(25-(dv)))/(exp((25-dv)/10)-1); 
      bm=4/(exp(dv/18));
      ah=0.07/(exp(dv/20));
      bh=1/(exp((30-dv)/10)+1); 

      %#5 applico Eulero per calcolare n m e h al passo successivo  
      n(i+1)= n(i)+dt*(an*(1-n(i))-bn*n(i));
      m(i+1)= m(i)+dt*(am*(1-m(i))-bm*m(i));
      h(i+1)= h(i)+dt*(ah*(1-h(i))-bh*h(i));
      
    
      
end
  
 
%% Grafico la corrente di stimolo e il potenziale di azione

figure
plot(t,I);
xlabel('t [msec]');
legend('IS');
ylabel('[\mu A/cm^2]'); %mS*mV

figure
plot(t,Vm);
ylabel('[mV]');
xlabel('t [msec]');
hold on
plot(t,ones(size(Vm))*Vr,'-.');
legend('Vm','Vr');
ylim([-90 40]);


%% Grafico  conduttanze

figure
subplot(3,1,1)
plot(t,GK);
hold on
plot(t,GNa);
scaleV=40/100;
plot(t,(Vm-Vr)*scaleV,'k--');
legend('GK','GNa','\Delta Vm');
xlabel('t [msec]');
ylabel('mS/cm^2');
subplot(3,1,2)
plot(t,h)
hold on
plot(t,m)
xlabel('t [msec]');
legend('h','m');
subplot(3,1,3)
plot(t,n)
xlabel('t [msec]');
legend('n');

INa=GNa.*(Vm-VNa);
IK=GK.*(Vm-VK);
figure
plot(t,IK)
hold on
plot(t,INa);
legend('INa','IK');


%% refrattario

%  figure
%  subplot(2,1,1)
%  plot(t,h)
%  legend('h')
%  xlabel('t [msec]');
%  subplot(2,1,2)
%  plot(t,I)
%  xlabel('t [msec]');