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