%SIMULAZIONE DEL POTENZIALE DI AZIONE 
% calcolare l'andamento del potenziale di membrana quando viene applicato
% un impulso di corrente di carica data 
%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;

%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

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


l=20; %cm 
dx=0.1;
stepx=l/dx; 
x=(0:stepx-1)*dx;

a=0.0238; %cm
R=34.5/1000; %Resistività intracellulare kohm*cm

%% inizializzo le variabili con la giusta dimensione 

Vm = ones(steps,stepx)*V0;
m = ones(steps,stepx)*m0;
h = ones(steps,stepx)*h0;
n = ones(steps,stepx)*n0;
%GK=zeros(steps,stepx);
%GNa=zeros(steps,stepx);
%GTot=zeros(steps,stepx); %conduttanza totale
%IK=zeros(steps,stepx);   %corrente potassio
%INa=zeros(steps,stepx);  %corrente sodio 
%IL=zeros(steps,stepx);   %corrente leakege
%Ii=zeros(steps,stepx);   %corrente ionica 
I=zeros(steps,1);    %corrente di membrana che impongo dall'esterno (ingresso forzante) 
IM=zeros(steps,stepx);

  
 
%%  Corrente di stimolazione

%si fissa la "carica" dello stimolo, se poi fissiamo la durata otteniamo la
%corrente Is (carica integrale della corrente nel tempo) 

stimolo=10e-9; %C
durata=0.2; %msec
IS=(stimolo)/(durata*1e-3); %corrente in A
ncampioni=round(durata/dt); %numero di passi o campioni relativi alla durata dell'impulso
IS=IS*10^6;% [uA/cm2]
I(5*ncampioni:5*ncampioni+ncampioni+ncampioni)=IS;
%figure
%plot(t,I);
  
  
%%  setto le condizioni iniziali


for i=1:steps-1 %tempo
for j=1:stepx-1 %Spazio





    
    %conduttanze
      GK=GKM*(n(i,j)^4);
      GNa=GNaM*(m(i,j)^3)*h(i,j);
         
      %correnti ioniche
      IK=GK*(Vm(i,j)-VK); %K
      INa=GNa*(Vm(i,j)-VNa); %Na
      IL=GL*(Vm(i,j)-VL); %L
      Ii=IK+INa+IL; % corrente ionica 
      
      if j==1
          II=1/dx*(I(i)/(2*pi*a)-(a/(R*dx))*(Vm(i,j)-Vm(i,j+1)));
      else
          II=(a/(2*R*dx^2))  * (Vm(i,j-1)-2*Vm(i,j)+Vm(i,j+1));     
      end
      IM(i,j)=II; 
        
          
      Vm(i+1,j)=Vm(i,j)+dt*( II- Ii)/CM;
      
      
      dv=Vm(i,j)-Vr; 
      
      %calcolo alpha beta 
      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); 

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

end




s1=1000; %fisso un certo istante temporale s1*dt
figure
plot(x,Vm(s1,:));
title(['t=' num2str(t(s1)) ' msec']);
xlabel('x');
ylabel('Vm');

figure
GNa=GNaM*(m.^3.*h);
GK=GKM*(n.^4);
Ii=GNa.*(Vm-VNa)+GK.*(Vm-VK)+GL*(Vm-VL);
Ic=IM-Ii;
subplot(3,1,1)
plot(x,Vm(s1,:));
title(['t=' num2str(t(s1)) ' msec']);
hold on 
plot(x,Vr*ones(size(Vm(s1,:))),'-.');
xlabel('x');
ylabel('Vm');
subplot(3,1,2)
plot(x,GNa(s1,:));
hold on
plot(x,GK(s1,:));
xlabel('x');
ylabel('Conduttanze');
legend('GNa','Gk');
subplot(3,1,3)
plot(x,Ii(s1,:));
hold on
plot(x,Ic(s1,:));
plot(x,IM(s1,:));
xlabel('x');
ylabel('Correnti');
legend('Ii','Ic','Im');


sx1=51; %fisso x=sx1*dx
figure
plot(t,Vm(:,sx1));
xlabel('t');
ylabel('Vm');
hold on
plot(t,Vm(:,2*sx1));
plot(t,Vr*ones(size(Vm(:,100))),'-.')
legend(    ['x=' num2str(x(sx1))]   ,  ['x=' num2str(x(2*sx1))]  );

figure 
subplot(3,1,1)
plot(t,Vm(:,sx1));
title(['x=' num2str(x(sx1)) ' cm']);
xlabel('t');
ylabel('Vm');
hold on
plot(t,Vr*ones(size(Vm(:,100))),'-.')
subplot(3,1,2)
plot(t,GNa(:,sx1));
hold on
plot(t,GK(:,sx1));
xlabel('t');
ylabel('Conduttanze');
legend('GNa','Gk');
subplot(3,1,3)
plot(t,Ii(:,sx1));
hold on
plot(t,Ic(:,sx1));
plot(t,IM(:,sx1));
ylabel('Correnti');
legend('Ii','Ic','Im');
xlabel('t');
ylabel('Im');


figure
s1=501;
plot(x,Vm(s1,:));
hold on
s2=1001;
plot(x,Vm(s2,:));
s3=1501;
plot(x,Vm(s3,:));
legend(    ['t=' num2str(t(s1))]   ,  ['t=' num2str(t(s2))] , ['t=' num2str(t(s3))]   );
xlabel('x');
ylabel('Vm');





