clc
clear all
close all








%%


Data = csvread('nodexg0changed_ALLExtraVoltages_stimulateonlyAbeta0_edgedist0.1_.csv');  % time , Abeta0_vext0_node0515 ,  Abeta0_vext1_node0515  , Abeta0_vext1_node015  , Abeta0_vext1_node115  , Abeta1_vext1_node0515 , boundary0_vext1_section1663 ,  Abeta0_vext0_node115 ,  Abeta0_vext0_MYSA0530 ,  Abeta0_vext1_MYSA0530 , Abeta0_vext0_node015



time = Data(:,1) ;

[satr, soton] = size(Data) ;

Vext5 = Data(:,2) ;
Vext1 = Data(:,3) ;
Vext4 = Data(:,4) ;
Vext2 = Data(:,5) ;
Vext3 = Data(:,6) ;
boundary = Data(:,7)  ;

Vext6 = Data(:,8) ;
Vext7 = Data(:,9) ;
Vext8 = Data(:,10) ;
Vext9 = Data(:,11) ;



%
XG =  7.08e+04 *1e-8  ;   %S/micron2
xr = 95769.75706051444*1e6*1e-4 ;   %ohm/micron
Rpn0 = 429128.527578172*1e6*1e-4 ;  %ohm/micron
mycm=0.1    ;
mygm=0.001  ;
nl = 104 ;

xg0_MYSA  = (mygm/(nl*2))*1e-8 ;  %S/micron2
C_MYSA = (mycm/(nl*2))*1e-6*1e-8   ;     %Farad/micron2

V78 = (Vext7 - Vext8) ;
dvdt78 = gradient(V78(:)) ./ gradient(time(:));

I1 = (Vext2-Vext1)/xr  ;    %mA.micron
I2 = -(Vext1-Vext3)*XG ;     %mA/micron2
I3 = (Vext4-Vext1)/xr   ;    %mA.micron
I4 = -(Vext5-Vext1)*XG       ;  %mA/micron2

IB = -(Vext1 - boundary)*16.1*1e-8 ;    %mA/micron2

IG = -Vext1*1e-9 *1e-8     ;  %mA/micron2

I6 = ((Vext6 - Vext5)/Rpn0)  ;  %mA.micron
I7 = ((Vext7 - Vext6)/Rpn0)  ;  %mA.micron
I8 = (Vext8 - Vext2)/xr  ;     %mA.micron
I9 = (Vext7 - Vext8)* xg0_MYSA  ;   %mA/micron2
I10 = C_MYSA * dvdt78  ;    %mA/micron2

I11 = (Vext9 - Vext5)/Rpn0  ;  %mA.micron

%

%%%% mA
IG = IG* 2.594156*pi ; 
I2 = I2* 2.594156*pi ;   %mA
I4 = I4* 2.594156*pi ;
IB = IB* 2.594156*pi ;

I9 = I9* 3*pi*8 ;
I10 = I10* 3*pi*8 ;


%%%%% mA
I1 = I1/0.5  ;
I3 = I3/0.5  ;
I6 = I6/0.5  ;
I11 = I11/0.5  ;

I7 = I7/1.5  ;
I8 = I8/1.5  ;




 


% % all currents
% figure
% 
% 
% 
% plot(time, I1*1e6  , 'LineWidth' , 2)   %nano amper
% hold on
% % 
% plot(time, I3*1e6, 'r' , 'LineWidth' , 2)
% % hold on
% plot(time, I6*1e6, 'c' , 'LineWidth' , 2)
% hold on
% 
% % plot(time, I11*1e6 , 'LineWidth' , 2)
% % hold on
% % 
% plot(time, I7*1e6, 'g' , 'LineWidth' , 2)
% hold on
% plot(time, I8*1e6, 'k' , 'LineWidth' , 2)
% hold on
% 
% 
% %%%%%%%
% 
% plot(time, I2*1e6 ,  'LineWidth' , 2)
% hold on
% % 
% plot(time, I4*1e6 , 'b' , 'LineWidth' , 2)
% hold on
% % 
% % 
% plot(time, IB*1e6  , 'm' , 'LineWidth' , 2)
% hold on
% 
% plot(time, I9*1e6 , 'LineWidth' , 2)
% hold on
% % 
% 
% plot(time, I10*1e6 , 'LineWidth' , 2)
% hold on
% 
% 
% 
% 
% 
% % ylim([-2*1e-7 1*1e-7])
%   %I5 = IB
% xlabel('time (ms)')
% ylabel('Current (nA)')
% xlim([0.95 2])
% 
% 
% legend('I1', 'I3' , 'I6' , 'I11' , 'I7' , 'I8' , 'I2' ,'I4' , 'IB' , 'I9' , 'I10')
% 
% 
% % Ic
% 
% figure
% 
% plot(time, I10*1e6 , 'LineWidth' , 2)
% hold on
% 
% xlabel('time (ms)')
% ylabel('Current (nA)')
% xlim([0.95 2])
% 
% % I1 and I3
% 
% figure
% 
% plot(time, I1*1e6  , 'LineWidth' , 2)   %nano amper
% hold on
% 
% plot(time, I3*1e6, 'r' , 'LineWidth' , 2)
% hold on
% 
% xlabel('time (ms)')
% ylabel('Current (nA)')
% xlim([0.95 2])
% 
% % I6 annd I11
% figure
% 
% plot(time, I6*1e6, 'c' , 'LineWidth' , 2)
% hold on
% 
% plot(time, I11*1e6 , 'LineWidth' , 2)
% hold on
% 
% xlabel('time (ms)')
% ylabel('Current (nA)')
% xlim([0.95 2])
% ylim([-0.1 0.8])

%% check if it is correct

% 
% figure
% III2 = IB + I2 +IG ;
% plot(time , III2  , 'g' ,  'LineWidth' , 2)
% hold on
% plot(time , I1+I3+I4  , 'LineWidth' , 2)
% 


















%% Analysis for neuron connected to Ground, NEURON is used for this simulation, only one fiber

% clc
% clear all
% close all




GData = csvread('ConnectedtoGround_ALLExtraVoltages.csv');


Gtime = GData(:,1) ;

[satr, soton] = size(GData) ;

GVext5 = GData(:,2) ;
GVext1 = GData(:,3) ;
GVext4 = GData(:,4) ;
GVext2 = GData(:,5) ;
GVext3 = GData(:,6) ;
Gboundary = GData(:,7)  ;

GVext6 = GData(:,8) ;
GVext7 = GData(:,9) ;
GVext8 = GData(:,10) ;
GVext9 = GData(:,11) ;



%
GXG =  1e9 *1e-8  ;   %S/micron2
Gxr = 1e9*1e6*1e-4 ;   %ohm/micron
GRpn0 = 429128.527578172*1e6*1e-4 ;  %ohm/micron
Gmycm=0.1    ;
Gmygm=0.001  ;
Gnl = 104 ;

Gxg0_MYSA  = (Gmygm/(Gnl*2))*1e-8 ;  %S/micron2
GC_MYSA = (Gmycm/(Gnl*2))*1e-6*1e-8   ;     %Farad/micron2

GV78 = (GVext7 - GVext8) ;
Gdvdt78 = gradient(GV78(:)) ./ gradient(Gtime(:));

GI1 = (GVext2-GVext1)/Gxr  ;    %mA.micron
GI2 = -(GVext1-GVext3)*GXG ;     %mA/micron2
GI3 = (GVext4-GVext1)/Gxr   ;    %mA.micron
GI4 = -(GVext5-GVext1)*GXG       ;  %mA/micron2

GIB = -(GVext1 - Gboundary)*16.1*1e-8 ;    %mA/micron2

GIG = -GVext1*1e9 *1e-8     ;  %mA/micron2

GI6 = ((GVext6 - GVext5)/GRpn0)  ;  %mA.micron
GI7 = ((GVext7 - GVext6)/GRpn0)  ;  %mA.micron
GI8 = (GVext8 - GVext2)/Gxr  ;     %mA.micron
GI9 = (GVext7 - GVext8)* Gxg0_MYSA  ;   %mA/micron2
GI10 = GC_MYSA * Gdvdt78  ;    %mA/micron2

GI11 = (GVext9 - GVext5)/GRpn0  ;  %mA.micron

%

%%%% mA
GIG = GIG* 2.594156*pi ; 
GI2 = GI2* 2.594156*pi ;   %mA
GI4 = GI4* 2.594156*pi ;
GIB = GIB* 2.594156*pi ;

GI9 = GI9* 3*pi*8 ;
GI10 = GI10* 3*pi*8 ;


%%%%% mA
GI1 = GI1/0.5  ;
GI3 = GI3/0.5  ;
GI6 = GI6/0.5  ;
GI11 = GI11/0.5  ;

GI7 = GI7/1.5  ;
GI8 = GI8/1.5  ;




 

% 
% % all currents
%  figure
% 
% 
% 
% plot(Gtime, GI1*1e6  , 'LineWidth' , 2)   %nano amper
% hold on
% % 
% plot(Gtime, GI3*1e6, 'r' , 'LineWidth' , 2)
% hold on
% plot(Gtime, GI6*1e6, 'c' , 'LineWidth' , 2)
% hold on
% % 
% plot(Gtime, GI11*1e6 , 'LineWidth' , 2)
% hold on
% 
% plot(Gtime, GI7*1e6, 'g' , 'LineWidth' , 2)
% hold on
% plot(Gtime, GI8*1e6, 'k' , 'LineWidth' , 2)
% hold on
% % 
% 
% %%%%%%%
% 
% % plot(time, I2*1e6 ,  'LineWidth' , 2)
% % hold on
% 
% plot(Gtime, GI4*1e6 , 'b' , 'LineWidth' , 2)
% hold on

% 
% % plot(time, IB*1e6  , 'm' , 'LineWidth' , 2)
% % hold on
% 
% plot(Gtime, GI9*1e6 , 'LineWidth' , 2)
% hold on
% % 
% % 
% plot(Gtime, GI10*1e6 , 'LineWidth' , 2)
% hold on
% 
% plot(Gtime, GIG*1e6 , 'LineWidth' , 2)
% hold on

% 
% 
%  ylim([-2*1e-7 1*1e-7])
% %   %I5 = IB
% xlabel('time (ms)')
% ylabel('Current (nA)')
% xlim([0.95 2])
% 
% 
% legend('I1', 'I3' , 'I6' , 'I11' , 'I7' , 'I8' ,'I4' , 'I9' , 'I10' , 'IG')
% 
% 
% 
% %%%%%%%% Ic
% 
% figure
% 
% plot(Gtime, GI10*1e6 , 'LineWidth' , 2)
% hold on
% 
% xlabel('time (ms)')
% ylabel('Current (nA)')
% xlim([0.95 2])
% 
% %%%%%%%%% I1 and I3
% 
% figure
% 
% plot(Gtime, GI1*1e6  , 'LineWidth' , 2)   %nano amper
% hold on
% 
% plot(Gtime, GI3*1e6, 'r' , 'LineWidth' , 2)
% hold on
% 
% xlabel('time (ms)')
% ylabel('Current (nA)')
% xlim([0.95 2])
% 
% %%%%%%%% I6 annd I11
% figure
% 
% plot(Gtime, GI6*1e6, 'c' , 'LineWidth' , 2)
% hold on
% 
% plot(Gtime, GI11*1e6 , 'LineWidth' , 2)
% hold on
% 
% xlabel('time (ms)')
% ylabel('Current (nA)')
% xlim([0.95 2])
% 



%% Compare ground condition with finite case


figure
plot(time, I1*1e6  , 'k', 'LineWidth' , 2) 
hold on
plot(Gtime, GI1*1e6  , 'g','LineWidth' , 2) 
title('IL (right)')
xlim([0 0.12])
 xlabel('time (ms)')
ylabel('Current (nA)')
%%
figure
plot(time, I3*1e6,'k',  'LineWidth' , 2)
hold on
plot(Gtime, GI3*1e6, 'g', 'LineWidth' , 2)
title('IL (left)')
xlim([0 0.12])
 xlabel('time (ms)')
ylabel('Current (nA)')

figure
plot(time, I6*1e6, 'k', 'LineWidth' , 2)
hold on
plot(Gtime, GI6*1e6, 'g', 'LineWidth' , 2)
title('Ip (right)')
xlim([0 0.12])
 xlabel('time (ms)')
ylabel('Current (nA)')


figure
plot(time, I11*1e6 ,'k', 'LineWidth' , 2)
hold on
plot(Gtime, GI11*1e6 ,'g', 'LineWidth' , 2)
title('Ip (left)')
xlim([0 0.12])
 xlabel('time (ms)')
ylabel('Current (nA)')

% figure
% plot(time, I7*1e6, 'k', 'LineWidth' , 2)
% hold on
% plot(Gtime, GI7*1e6,'g',  'LineWidth' , 2)
% title('I7')
% xlim([0 1])
%  xlabel('time (ms)')
% ylabel('Current (nA)')

% figure
% plot(time, I8*1e6,'k',  'LineWidth' , 2)
% hold on
% plot(Gtime, GI8*1e6, 'g', 'LineWidth' , 2)
% title('I8')
% xlim([0 1])
%  xlabel('time (ms)')
% ylabel('Current (nA)')


%%%%%%%
%%
figure
plot(time, I2*1e6  ,'k' , 'LineWidth' , 2)
hold on
plot(Gtime, GI2*1e6 ,'g' , 'LineWidth' , 2)
title('I nearby')
xlim([0 0.12])
  xlabel('time (ms)')
ylabel('Current (nA)')
 
 %%

figure
plot(time, I4*1e6 ,'k',  'LineWidth' , 2)
hold on
plot(Gtime, GI4*1e6 ,'g', 'LineWidth' , 2)
title('IT')
xlim([0 0.12])
  xlabel('time (ms)')
ylabel('Current (nA)')
 
 
 

 
 %%

figure
plot(time, I9*1e6 ,'k', 'LineWidth' , 2)
hold on
plot(Gtime, GI9*1e6 ,'g', 'LineWidth' , 2)
title('Ileak_myelin')
xlim([0 0.12])
  xlabel('time (ms)')
ylabel('Current (nA)')
 %%



figure
plot(time, I10*1e6 ,'k', 'LineWidth' , 2)
hold on
plot(Gtime, GI10*1e6 ,'g', 'LineWidth' , 2)
title('Icap_myelin')

xlim([0 0.12])
 xlabel('time (ms)')
ylabel('Current (nA)')


%%
figure
plot(time, IB*1e6 ,'k' , 'LineWidth' , 2)
hold on
plot(Gtime, GIG*1e6 ,'g' , 'LineWidth' , 2)
title('IG ')

xlim([0 0.12])
 xlabel('time (ms)')
ylabel('Current (nA)')



%%
% close all
%% 
% figure
% plot(time, Vext1 , 'k' , 'LineWidth' , 2)
% hold on
% plot(time, GVext1 , 'g' , 'LineWidth' , 2)
% 
% xlabel('time (ms)')
% ylabel('Extracellular V (mV)')
% xlim([0 1])
% ylim([-4e-8 1e-8])


%%
% 
% figure
% plot(time, Vext5 , 'k' , 'LineWidth' , 2)
% hold on
% plot(time, GVext5 , 'g' , 'LineWidth' , 2)
% 
% %%
% close all
% figure
% plot(time, Vext5 , 'k' , 'LineWidth' , 2)
% hold on
% plot(time, GVext5 , 'g' , 'LineWidth' , 2)


%%
% clc
% clear all
% close all

v_01 = csvread('v_Abeta0_stimulateonlyAbeta0_edgedist0.1_.csv');
v_01_ground = csvread('Onefiber_connectedGround_node_v.csv');

time = v_01(:,1) ; 


figure
plot(time, v_01(:,17) , 'k' , 'LineWidth' , 2)
hold on
plot(time, v_01_ground(:,17) , 'g' , 'LineWidth' , 2)
xlabel('time (ms)')
ylabel('Transmembrane V (mV)')
xlim([0 0.12])




%% 

% clc
% clear all

Surface_area_node = 2.59415600000000*pi*1e-8 ;  %cm2

imembrane_01 = csvread('imembrane_Abeta0_stimulateonlyAbeta0_edgedist0.1_.csv');
imembrane_01_ground = csvread('Onefiber_connectGround_nodes_imembrane.csv');

time = imembrane_01(:,1) ;

figure
plot(time, -imembrane_01(:,17)*Surface_area_node*1e6 , 'k' , 'LineWidth' , 2)
hold on
plot(time, -imembrane_01_ground(:,17)*Surface_area_node*1e6 , 'g' , 'LineWidth' , 2)
xlabel('time (ms)')
ylabel('Transmembrane Current (mA/cm2)')
xlim([0 0.12])


%% Iin

Iin = (-imembrane_01(:,17)*Surface_area_node*1e6)/2  ;

GIin = (-imembrane_01_ground(:,17)*Surface_area_node*1e6)/2  ;


figure
plot(time, (-imembrane_01(:,17)*Surface_area_node*1e6)/2 , 'k' , 'LineWidth' , 2)
hold on
plot(time, (-imembrane_01_ground(:,17)*Surface_area_node*1e6)/2 , 'g' , 'LineWidth' , 2)
xlabel('time (ms)')
ylabel('Transmembrane Current (mA/cm2)')
xlim([0 0.12])


%%  IT + Ip (right) + Ip(left)


I_sum = (I4 + I6 + I11)*1e6 ;
IG_sum = (GI4 + GI6 + GI11)*1e6 ;


figure
plot(time, I_sum , 'k' , 'LineWidth' , 2)
hold on
plot(time, IG_sum , 'g' , 'LineWidth' , 2)
xlabel('time (ms)')

xlim([0 0.12])
title('Sum')


%%
figure
plot(time, I10 , 'k' , 'LineWidth' , 2)
hold on
plot(time, GI10 , 'g' , 'LineWidth' , 2)
xlabel('time (ms)')
xlim([0 0.12])
title('I10')


%%
MYSA30_v = csvread('xg0changed_MYSA30_05_1_v_stimulateonlyAbeta0_edgedist0.1_.csv');
MYSA30_vext0 = csvread('xg0changed_MYSA30_05_1_vext0_stimulateonlyAbeta0_edgedist0.1_.csv');


M = ( (MYSA30_v(:,2)+MYSA30_vext0(:,2)) - (MYSA30_v(:,3)+MYSA30_vext0(:,3)) ) / (0.19865854  *1e6) ;   %mA

M = M*1e6  ;   %nA



GMYSA30_v = csvread('ConnectGround_MYSA30_v.csv');
GMYSA30_vext0 = csvread('ConnectGround_MYSA30_vext0.csv');


GM =  ( (GMYSA30_v(:,2)+GMYSA30_vext0(:,2)) - (GMYSA30_v(:,3)+GMYSA30_vext0(:,3)) ) / (0.19865854  *1e6) ;   %mA
 
GM =  GM *1e6  ;

figure
plot(time, M , 'k' , 'LineWidth' , 2)
hold on
% plot(time, GM , 'g' , 'LineWidth' , 2)
% xlabel('time (ms)')
xlim([0 0.12])
% title('M')


%%




Surface_area_MYSA = 8*3*pi*1e-8 ;  %cm2

GMYSA30_il_icap = csvread('ConnectGround_MYSA30_ipas_icap.csv');

GMYSA30_il = GMYSA30_il_icap(:,2)*Surface_area_MYSA*1e6 ;

GMYSA30_icap = GMYSA30_il_icap(:,3)*Surface_area_MYSA*1e6 ;

figure
plot(time, GMYSA30_icap , 'g' , 'LineWidth' , 2)
xlim([0 0.12])
hold on



%%

Surface_area_MYSA = 8*3*pi*1e-8 ;  %cm2

MYSA30_il_icap = csvread('xg0changed_MYSA30_il_icap_stimulateonlyAbeta0_edgedist0.1_.csv');

MYSA30_il = MYSA30_il_icap(:,2)*Surface_area_MYSA*1e6 ;

MYSA30_icap = MYSA30_il_icap(:,3)*Surface_area_MYSA*1e6 ;


plot(time, MYSA30_icap , 'k' , 'LineWidth' , 2)
xlim([0 0.12])

%%

IS = Iin - MYSA30_icap+MYSA30_il ;

figure
plot(time, IS , 'k' , 'LineWidth' , 2)
xlim([0 0.12])
hold on
%%

GIS = GIin - GMYSA30_icap+GMYSA30_il ;

plot(time, GIS , 'g' , 'LineWidth' , 2)
xlim([0 0.12])












%% node ina , inap , icap , ik , ileak

node15_current = csvread('node15_currents.csv'); 
ConnectGround_node15_current = csvread('ConnectGround_node15_currents.csv'); 


ina = node15_current(:,1)*Surface_area_node*1e6  ;
inap = node15_current(:,2)*Surface_area_node*1e6  ;
ik = node15_current(:,3)*Surface_area_node*1e6  ;
il = node15_current(:,4)*Surface_area_node*1e6  ;
icap = node15_current(:,5)*Surface_area_node*1e6  ;


Gina = ConnectGround_node15_current(:,1)*Surface_area_node*1e6  ;
Ginap = ConnectGround_node15_current(:,2)*Surface_area_node*1e6  ;
Gik = ConnectGround_node15_current(:,3)*Surface_area_node*1e6  ;
Gil = ConnectGround_node15_current(:,4)*Surface_area_node*1e6  ;
Gicap = ConnectGround_node15_current(:,5)*Surface_area_node*1e6  ;

% 
% figure
% plot(time, ina , 'k' , 'LineWidth' , 2)
% hold on
% plot(time, Gina , 'g' , 'LineWidth' , 2)
% 
% title('ina')
% xlim([0 0.12])
% 
% 
% figure
% plot(time, inap , 'k' , 'LineWidth' , 2)
% hold on
% plot(time, Ginap , 'g' , 'LineWidth' , 2)
% 
% title('inap')
% xlim([0 0.12])
% 
% 
% figure
% plot(time, ik , 'k' , 'LineWidth' , 2)
% hold on
% plot(time, Gik , 'g' , 'LineWidth' , 2)
% 
% title('ik')
% xlim([0 0.12])
% 
% 
% figure
% plot(time, il , 'k' , 'LineWidth' , 2)
% hold on
% plot(time, Gil , 'g' , 'LineWidth' , 2)
% 
% title('il')
% xlim([0 0.12])
% 
% 
% 
% figure
% plot(time, icap , 'k' , 'LineWidth' , 2)
% hold on
% plot(time, Gicap , 'g' , 'LineWidth' , 2)
% 
% title('icap')
% xlim([0 0.12])

%%


figure
plot(time, (ina+inap+ik+il+icap) , 'k' , 'LineWidth' , 2)
hold on
plot(time, (Gina+Ginap+Gik+Gil+Gicap), 'g' , 'LineWidth' , 2)

title('sum')
xlim([0 0.12])



