Effect of ionic diffusion on extracellular potentials (Halnes et al 2016)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:225311
"Recorded potentials in the extracellular space (ECS) of the brain is a standard measure of population activity in neural tissue. Computational models that simulate the relationship between the ECS potential and its underlying neurophysiological processes are commonly used in the interpretation of such measurements. Standard methods, such as volume-conductor theory and current-source density theory, assume that diffusion has a negligible effect on the ECS potential, at least in the range of frequencies picked up by most recording systems. This assumption remains to be verified. We here present a hybrid simulation framework that accounts for diffusive effects on the ECS potential. ..."
Reference:
1 . Halnes G, Mäki-Marttunen T, Keller D, Pettersen KH, Andreassen OA, Einevoll GT (2016) Effect of Ionic Diffusion on Extracellular Potentials in Neural Tissue. PLoS Comput Biol 12:e1005193 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Extracellular; Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Neocortex U1 L6 pyramidal corticalthalamic GLU cell;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: MATLAB; NEURON;
Model Concept(s): Extracellular Fields;
Implementer(s): Halnes, Geir [geir.halnes at nmbu.no]; Maki-Marttunen, Tuomo [tuomomm at uio.no];
Search NeuronDB for information about:  Neocortex U1 L6 pyramidal corticalthalamic GLU cell; 
function PLoSdynplot(Sdiff,Snodiff)
%close all; 
c1 = 'r';
c2 = 'b';

h = figure;
set(h, 'Tag', 'vprofiles');
set(h,'Color', [1 1 1]);
set(h, 'Position', [40 370 1200  430]);

h2 = figure;
set(h2, 'Tag', 'cprofiles2');
set(h2,'Color', [1 1 1]);
set(h2, 'Position', [40 370 1200  430]);

hv = figure;
set(hv, 'Tag', 'vvprofiles');
set(hv,'Color', [1 1 1]);
set(hv, 'Position', [40 370 1200  430]);

hh = figure; 
set(hh, 'Tag', 'jprofiles');
set(hh,'Color', [1 1 1]);
set(hh, 'Position', [57    49   914   767]);

hhhh = uicontrol('Style','text');
set(hhhh, 'Units', 'Normalized', 'BackgroundColor', [1 1 1]);
set(hhhh, 'String', 'A', 'Position', [0.06 0.9 0.05 0.05], 'Fontsize', 14);

hhhh = uicontrol('Style','text');
set(hhhh, 'Units', 'Normalized', 'BackgroundColor', [1 1 1]);
set(hhhh, 'String', 'B1', 'Position', [0.06 0.67 0.05 0.05], 'Fontsize', 14);

hhhh = uicontrol('Style','text');
set(hhhh, 'Units', 'Normalized', 'BackgroundColor', [1 1 1]);
set(hhhh, 'String', 'B2', 'Position', [0.06 0.45 0.05 0.05], 'Fontsize', 14);

hhhh = uicontrol('Style','text');
set(hhhh, 'Units', 'Normalized', 'BackgroundColor', [1 1 1]);
set(hhhh, 'String', 'B3', 'Position', [0.06 0.22 0.05 0.05], 'Fontsize', 14);


plotcprof(Snodiff, 1, c1);
plotcprof(Sdiff, 2, c2);
figure(hh);

% Make some text on top
subplot(4,6,2:6);
separator0 = '          ';
separator = '               ';
separator2 = '             ';
separator3 = '                        ';
fmytext = [separator0, 'j^{Na}', separator, 'j^{K}', separator, 'j^{Ca}', separator2, ' j^{X}', separator3, 'i/F'];
haxis = axis;
text(haxis(1),haxis(4)*1.2, fmytext, 'Fontsize', 16)


figure(h);
subplot(1,2,1);
title(['A'], 'Fontsize', 14);
subplot(1,2,2);
xlabel('i^M (nA)', 'Fontsize', 12);
title(['B'], 'Fontsize', 14);

figure(hv);
subplot(1,3,1);
title(['A'], 'Fontsize', 14);
subplot(1,3,2);
title(['B (diff. off)'], 'Fontsize', 14);
xlabel('V (mV)', 'Fontsize', 12);
subplot(1,3,3);
title(['C: (diff. on)'], 'Fontsize', 14);
xlabel('V (mV)', 'Fontsize', 12);


figure(h2);
for i = 1:5
subplot(2,5,i);
title(['A',num2str(i)], 'Fontsize', 14);
subplot(2,5,5+i);
title(['B',num2str(i)], 'Fontsize', 14);
end

subplot(2,5,7); xlabel('[Na^+] (mM)', 'Fontsize', 12);
subplot(2,5,8); xlabel('[K^+] (mM)', 'Fontsize', 12);
subplot(2,5,9); xlabel('[Ca^{2+}] (mM)', 'Fontsize', 12);
subplot(2,5,10); xlabel('[X^{-}] (mM)', 'Fontsize', 12);


function plotcprof(S, crow, Dcol);
% Useful parameters
C_m = 1.00e-2; % Membrane capacitance (Farad/m^2);
F = 96485.3365; % C/mol
T = 300; % K
R = 8.3; % J/mol/K
psi = R*T/F; % V


N = S.Nvox;
deltax = S.geometry.deltax;
Avox = S.geometry.Avox;
t = S.Simdata.t;
tv = S.Simdata.t;
cNa = S.Simdata.cNa;
cK = S.Simdata.cK;
cCa = S.Simdata.cCa;
cX = S.Simdata.cX;
V = S.Simdata.V;
diffconsts = S.diffconsts; 
lambda_o = diffconsts(1);
D_K = diffconsts(2);
D_Na = diffconsts(3);
D_Ca = diffconsts(4);
D_X = diffconsts(5);
diffon = S.diffon;

jk = S.Neurondata.jk;
jna = S.Neurondata.jna;
jca = S.Neurondata.jca;
jx = S.Neurondata.jx;
icap = S.Neurondata.icap;
%isyn = S.Neurondata.isyn;
imemb = S.Neurondata.imemb;
times = S.Neurondata.times; 

tv = tv(2:end-1); % just to eliminate some endpoint bugs
V = V(2:end-1,:);
V = [V(1,:); V(:,:); V(end,:)]';
cK = cK'; cNa = cNa'; cCa = cCa'; cX = cX'; 
Npsi = 0.0258;

% Extracellular fluxes: (mol/s)
jKEd = - diffon*Avox*( D_K*(cK(2:end,:) - cK(1:end-1,:))/deltax);
jKEf = - Avox*D_K/Npsi*(cK(2:end,:) + cK(1:end-1,:))/2.*(V(2:end,:) - V(1:end-1,:))/deltax; % N+1 elements
jNaEd =  - diffon*Avox*( D_Na*(cNa(2:end,:) - cNa(1:end-1,:))/deltax);
jNaEf = - Avox*D_Na/Npsi*(cNa(2:end,:) + cNa(1:end-1,:))/2.*(V(2:end,:) - V(1:end-1,:))/deltax; % N+1 elements
jCaEd = - diffon*Avox*( D_Ca*(cCa(2:end,:) - cCa(1:end-1,:))/deltax);
jCaEf = - Avox*2*D_Ca/Npsi*(cCa(2:end,:) + cCa(1:end-1,:))/2.*(V(2:end,:) - V(1:end-1,:))/deltax; % N+1 elements
jXEd =  - diffon*Avox*( D_X*(cX(2:end,:) - cX(1:end-1,:))/deltax);
jXEf = + Avox*D_X/Npsi*(cX(2:end,:) + cX(1:end-1,:))/2.*(V(2:end,:) - V(1:end-1,:))/deltax; % N+1 elements

% Extracellular currents: (A)
iEd = F*(2*jCaEd + jKEd + jNaEd - jXEd);
iEf = F*(2*jCaEf + jKEf + jNaEf - jXEf);
%iEd = iEd/F; iEf = iEf/F;

iEdin = iEd(1,:);
iEdout = iEd(end,:);
iEfin = iEf(1,:);
iEfout = iEf(end,:);
netdin = iEdin-iEdout;
netfin = iEfin-iEfout;
netin = netdin + netfin;

cK = cK'; cNa = cNa'; cCa = cCa'; cX = cX';

comps = 1:size(cNa,2); % number of compartments
timepts = round(linspace(2,length(t),6)); % Selected time pts.

%avgisyn = mean(isyn, 2);
%avgisyn = [0; avgisyn; 0];
ints = 5;
intt = 7/5;
imtpts = 0:intt:7;

inti = round(length(imemb)/ints);
avgimemb = zeros(ints, size(imemb,1)+2);

for i = 1:ints;
    myint = ((i-1)*inti+1):(i*inti-1);
    avgimemb(i,:) = [0;mean(imemb(:,myint),2);0];

end
zeroline = zeros(1,size(imemb,1)+2);


Mycolors = [0    0.4470    0.7410;
    0.8500    0.3250    0.0980;
    0.9290    0.6940    0.1250;
    0.4940    0.1840    0.5560;
    0.4660    0.6740    0.1880];

%%%%%%%%% PLOT SPATIAL PROFILES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% FIGURE %%%%%%%
%%%% Spatial profiles %%%%%%
h = findobj('Tag', 'vprofiles');
figure(h);

subplot(1,2,1);
hold on;
% make a drawing of the neuron
plot(3,3,'^', 'markersize', 20, 'Linewidth', 0.5, 'MarkerFaceColor', Dcol ,'Color', Dcol);
plot([3 1.7], [3 1.7], '-', 'Linewidth', 2, 'Color', Dcol);
plot([3 4.3], [3 1.7], '-', 'Linewidth', 2, 'Color', Dcol)
plot([3 3], [3 13], '-', 'Linewidth', 4, 'Color', Dcol)
plot([3 1.5], [13 14], '-', 'Linewidth', 2, 'Color', Dcol)
plot([3 4.5], [13 14], '-', 'Linewidth', 2, 'Color', Dcol)
%plot(3,3,'^', 'markersize', 8, 'Linewidth', 8, 'Color', Dcol);

axis([0 6 1 15]);
%set(gca, 'Yticklabel', '')
set(gca, 'Xticklabel', '')

subplot(1,2,2); 
set(gca, 'ColorOrder', Mycolors);      
hold on;
%plot(avgisyn*1e9, comps, 'Color', [0.5 0.5 0.5]);
plot(avgimemb*1e9, comps);
plot(zeroline, comps, 'k:');
%haxis = axis; haxis(1) = -0.5; haxis(2) = 1; haxis(3) = 1; haxis(4) = 15; axis(haxis);
Lfroms = num2str(imtpts(1:end-1)');
Ltos = num2str(imtpts(2:end)');
mylegend = [Lfroms, ['     ']', ['-----']', ['     ']', Ltos, ['sssss']'];
legend(mylegend);
set(gca, 'Yticklabel', '')



%%%% FIGURE %%%%%%%
%%%% Spatial profiles ONLY CONCENTRATIONS%%%%%%
h2 = findobj('Tag', 'cprofiles2');
figure(h2);

subplot(2,5,5*(crow-1)+1); 
hold on;
% make a drawing of the neuron
plot(3,3,'^', 'markersize', 10, 'Linewidth', 0.5, 'MarkerFaceColor', Dcol ,'Color', Dcol);
plot([3 1.7], [3 1.7], '-', 'Linewidth', 2, 'Color', Dcol);
plot([3 4.3], [3 1.7], '-', 'Linewidth', 2, 'Color', Dcol)
plot([3 3], [3 13], '-', 'Linewidth', 4, 'Color', Dcol)
plot([3 1.5], [13 14], '-', 'Linewidth', 2, 'Color', Dcol)
plot([3 4.5], [13 14], '-', 'Linewidth', 2, 'Color', Dcol)
axis([0 6 1 15]);
%set(gca, 'Yticklabel', '')
set(gca, 'Xticklabel', '')

subplot(2,5,5*(crow-1)+2); 
plot(cNa(timepts,:)', comps);
set(gca, 'Yticklabel', '')
haxis = axis; haxis(3) = 1; haxis(4) = 15; axis(haxis);

subplot(2,5,5*(crow-1)+3); 
plot(cK(timepts,:)', comps);
haxis = axis; haxis(3) = 1; haxis(4) = 15; axis(haxis);
set(gca, 'Yticklabel', '')

subplot(2,5,5*(crow-1)+4); 
plot(cCa(timepts,:)', comps);
haxis = axis; haxis(3) = 1; haxis(4) = 15; axis(haxis);
set(gca, 'Yticklabel', '')

subplot(2,5,5*(crow-1)+5); 
plot(cX(timepts,:)', comps);
haxis = axis; haxis(3) = 1; haxis(4) = 15; axis(haxis);
set(gca, 'Yticklabel', '')

if crow == 1
    legt2 = round(10*t(timepts(2:end)))/10; 
    legt1 = [0; legt2(1:end-1)];
    legt = [0;legt2];

    subplot(2,5,5); 
    mylegend = [num2str(legt), ['ssssss']'];
    legend(mylegend);
    set(gca, 'Yticklabel', '');
end



%%%% SPATIAL PROFILES ONLY V
hv = findobj('Tag', 'vvprofiles');
figure(hv);
subplot(1,3,1)
hold on;
Dcol = 'k';
% make a drawing of the neuron
plot(3,3,'^', 'markersize', 20, 'Linewidth', 0.5, 'MarkerFaceColor', Dcol ,'Color', Dcol);
plot([3 1.7], [3 1.7], '-', 'Linewidth', 2, 'Color', Dcol);
plot([3 4.3], [3 1.7], '-', 'Linewidth', 2, 'Color', Dcol)
plot([3 3], [3 13], '-', 'Linewidth', 4, 'Color', Dcol)
plot([3 1.5], [13 14], '-', 'Linewidth', 2, 'Color', Dcol)
plot([3 4.5], [13 14], '-', 'Linewidth', 2, 'Color', Dcol)
axis([0 6 1 15]);
%set(gca, 'Yticklabel', '')
set(gca, 'Xticklabel', '')

if crow == 1
    subplot(1,3,2)
end
if crow == 2
    subplot(1,3,3)
end
avgV = zeros(15,length(timepts)-1);
for ii = 2:length(timepts)
avgV(:,ii-1) = mean(V(:,timepts(ii-1):timepts(ii)),2);
end
hold on; 
plot(1000*avgV', comps);
    plot(zeroline, comps, 'k:');

set(gca, 'Yticklabel', '')
haxis = axis; haxis(1) = -1.5; haxis(2) = 0.2; haxis(3) = 1; haxis(4) = 15; axis(haxis);
%haxis = axis; haxis(3) = 1; haxis(4) = 15; axis(haxis);
set(gca, 'Yticklabel', '')

if 1 
if crow == 1
    subplot(1,3,3)
    avgV = mean(V(:,:),2);
    hold on;
    plot(1000*avgV', comps, 'k', 'LineWidth', 1.5);
%    plot(zeroline, comps, 'k:');
    set(gca, 'Yticklabel', '')
end    
end

if crow == 2
    subplot(1,3,3)
    legt2 = round(10*t(timepts(2:end)))/10;
    legt1 = [0; legt2(1:end-1)];
    mylegend = [num2str(legt1), ['     ']', ['-----']', ['     ']', num2str(legt2), ['sssss']'];
    mylegend = ['  diff. off '; mylegend];
    legend(mylegend);
    plot(zeroline, comps, 'k:');
    set(gca, 'Yticklabel', '')
end




%%%% FIGURE %%%%%%%
%%%% SPATIAL PROFILE OF FLUXES
avgjNaEd = zeros(14,length(timepts)-1);
avgjNaEf = zeros(14,length(timepts)-1);
avgjKEd = zeros(14,length(timepts)-1);
avgjKEf = zeros(14,length(timepts)-1);
avgjCaEd = zeros(14,length(timepts)-1);
avgjCaEf = zeros(14,length(timepts)-1);
avgjXEd = zeros(14,length(timepts)-1);
avgjXEf = zeros(14,length(timepts)-1);


for ii = 2:length(timepts)
avgjNaEd(:,ii-1) = mean(jNaEd(:,timepts(ii-1):timepts(ii)),2);
avgjNaEf(:,ii-1) = mean(jNaEf(:,timepts(ii-1):timepts(ii)),2);
avgjKEd(:,ii-1) = mean(jKEd(:,timepts(ii-1):timepts(ii)),2);
avgjKEf(:,ii-1) = mean(jKEf(:,timepts(ii-1):timepts(ii)),2);
avgjCaEd(:,ii-1) = mean(jCaEd(:,timepts(ii-1):timepts(ii)),2);
avgjCaEf(:,ii-1) = mean(jCaEf(:,timepts(ii-1):timepts(ii)),2);
avgjXEd(:,ii-1) = mean(jXEd(:,timepts(ii-1):timepts(ii)),2);
avgjXEf(:,ii-1) = mean(jXEf(:,timepts(ii-1):timepts(ii)),2);
end

% Extracellular currents: (A)
avgiEd = F*(2*avgjCaEd + avgjKEd + avgjNaEd - avgjXEd);
avgiEf = F*(2*avgjCaEf + avgjKEf + avgjNaEf - avgjXEf);
avgiEd = avgiEd/F; avgiEf = avgiEf/F; 
avgiE = avgiEd + avgiEf;

fcomps = comps(1:end-1); % One less

Mycolors = [0    0.4470    0.7410;
    0.8500    0.3250    0.0980;
    0.9290    0.6940    0.1250;
    0.4940    0.1840    0.5560;
    0.4660    0.6740    0.1880];


%%%%% FIELD STUFF
hh = findobj('Tag', 'jprofiles');
figure(hh);

if crow == 2
    subplot(4,6,7);
    hold on;
    % make a drawing of the neuron
    Dcol = 'b';
    plot(3,3,'^', 'markersize', 10, 'Linewidth', 0.5, 'MarkerFaceColor', Dcol ,'Color', Dcol);
    plot([3 1.7], [3 1.7], '-', 'Linewidth', 2, 'Color', Dcol);
    plot([3 4.3], [3 1.7], '-', 'Linewidth', 2, 'Color', Dcol)
    plot([3 3], [3 13], '-', 'Linewidth', 4, 'Color', Dcol)
    plot([3 1.5], [13 14], '-', 'Linewidth', 2, 'Color', Dcol)
    plot([3 4.5], [13 14], '-', 'Linewidth', 2, 'Color', Dcol)
    axis([0 6 1 15]);
    set(gca, 'Xticklabel', '')
    
    
    haxis = [-2.0000e-14 1.5400e-13 0.5 14.5];
    %xsep = 1.5e-14;
%    xsep = 2.8e-14;
    xsep = 3.2e-14;
    subplot(4,6,8:12);
    set(gca, 'ColorOrder', Mycolors);      
    hold on;
    plot(avgjNaEf, fcomps);
    plot([0 0], [1 15], 'k--');
    
    plot(avgjKEf+xsep, fcomps);
    plot([xsep xsep], [1 15], 'k--');
    
    plot(avgjCaEf+2*xsep, fcomps);
    plot([2*xsep 2*xsep], [1 15], 'k--');
    
    plot(avgjXEf+3*xsep, fcomps);
    plot([3*xsep 3*xsep], [1 15], 'k--');
    
    plot(avgiEf+4.5*xsep, fcomps);
    plot([4.5*xsep 4.5*xsep], [1 15], 'k--');
    
    axis(haxis);
    set(gca, 'Yticklabel', '')
    set(gca, 'Xticklabel', '')
    set(gca, 'Xtick',1000)
    set(gca, 'Visible', 'off')
    
    
    plot([-1 1],[0.5 0.5], 'k-');
    plot([-1 1],[14.5 14.5], 'k-');
    text(-2e-14, 11.5, 'field', 'Fontsize', 14)
    
    %set(gcf,'Position',[ 27 82 1001 733]);
    
    
    %%%%% DIFFUSIVE STUFF
    subplot(4,6,13);
    hold on;
    Dcol = 'b';
    % make a drawing of the neuron
    plot(3,3,'^', 'markersize', 10, 'Linewidth', 0.5, 'MarkerFaceColor', Dcol ,'Color', Dcol);
    plot([3 1.7], [3 1.7], '-', 'Linewidth', 2, 'Color', Dcol);
    plot([3 4.3], [3 1.7], '-', 'Linewidth', 2, 'Color', Dcol)
    plot([3 3], [3 13], '-', 'Linewidth', 4, 'Color', Dcol)
    plot([3 1.5], [13 14], '-', 'Linewidth', 2, 'Color', Dcol)
    plot([3 4.5], [13 14], '-', 'Linewidth', 2, 'Color', Dcol)
    axis([0 6 1 15]);
    %set(gca, 'Yticklabel', '')
    set(gca, 'Xticklabel', '')
    
    subplot(4,6,14:18);
    set(gca, 'ColorOrder', Mycolors);      
    hold on;
    
    plot(avgjNaEd, fcomps); %title('jNaEf')
    plot([0 0], [1 15], 'k--');

    plot(avgjKEd+xsep, fcomps); %title('jKEf')
    plot([xsep xsep], [1 15], 'k--');
    
    plot(avgjCaEd+2*xsep, fcomps); %title('jCaEf')
    plot([2*xsep 2*xsep], [1 15], 'k--');
    
    plot(avgjXEd+3*xsep, fcomps); %title('jXEf')
    plot([3*xsep 3*xsep], [1 15], 'k--');
    
    plot(avgiEd+4.5*xsep, fcomps); %title('iEf')
    plot([4.5*xsep 4.5*xsep], [1 15], 'k--');
    
    plot([-1 1],[0.5 0.5], 'k-');
    plot([-1 1],[14.5 14.5], 'k-');
    text(-2e-14, 11.5, 'diff', 'Fontsize', 14)
    
    axis(haxis);
    set(gca, 'Yticklabel', '')
    set(gca, 'Xticklabel', '')
    set(gca, 'Xtick',1000)
    set(gca, 'Visible', 'off')
    
    
    %%%%% TOTAL STUFF
    subplot(4,6,19);
    hold on;
    % make a drawing of the neuron
    Dcol = 'b';
    plot(3,3,'^', 'markersize', 10, 'Linewidth', 0.5, 'MarkerFaceColor', Dcol ,'Color', Dcol);
    plot([3 1.7], [3 1.7], '-', 'Linewidth', 2, 'Color', Dcol);
    plot([3 4.3], [3 1.7], '-', 'Linewidth', 2, 'Color', Dcol)
    plot([3 3], [3 13], '-', 'Linewidth', 4, 'Color', Dcol)
    plot([3 1.5], [13 14], '-', 'Linewidth', 2, 'Color', Dcol)
    plot([3 4.5], [13 14], '-', 'Linewidth', 2, 'Color', Dcol)
    axis([0 6 1 15]);
    %set(gca, 'Yticklabel', '')
    set(gca, 'Xticklabel', '')
    
    subplot(4,6,20:24);
    set(gca, 'ColorOrder', Mycolors);      
    hold on;
    plot(avgjNaEd + avgjNaEf, fcomps); %title('jNaE')
    plot([0 0], [1 15], 'k--');
    plot(avgjKEd  + avgjKEf + xsep, fcomps); %title('jKE')
    plot([xsep xsep], [1 15], 'k--');
    plot(avgjCaEd + avgjCaEf +2*xsep, fcomps); %title('jCaE')    
    plot([2*xsep 2*xsep], [1 15], 'k--');
    plot(avgjXEd + avgjXEf + 3*xsep, fcomps); %title('jXE')
    plot([3*xsep 3*xsep], [1 15], 'k--');
    plot(avgiE + 4.5*xsep, fcomps); %title('iE')
    plot([4.5*xsep 4.5*xsep], [1 15], 'k--');    
    plot([-1 1],[0.5 0.5], 'k-');
    plot([-1 1],[14.5 14.5], 'k-');
    text(-2e-14, 11.5, 'tot', 'Fontsize', 14)
    

    
    axis(haxis);
    set(gca, 'Yticklabel', '')
    set(gca, 'Xticklabel', '')
    set(gca, 'Xtick',1000)
    set(gca, 'Visible', 'off')
    
end
if crow == 1
    %%%%% TOTAL STUFF
    subplot(4,6,1);
    hold on;     
    haxis = [-2.0000e-14 1.5400e-13 0.5 14.5];
    xsep = 3.2e-14;

    % make a drawing of the neuron
    Dcol = 'r';
    plot(3,3,'^', 'markersize', 10, 'Linewidth', 0.5, 'MarkerFaceColor', Dcol ,'Color', Dcol);
    plot([3 1.7], [3 1.7], '-', 'Linewidth', 2, 'Color', Dcol);
    plot([3 4.3], [3 1.7], '-', 'Linewidth', 2, 'Color', Dcol)
    plot([3 3], [3 13], '-', 'Linewidth', 4, 'Color', Dcol)
    plot([3 1.5], [13 14], '-', 'Linewidth', 2, 'Color', Dcol)
    plot([3 4.5], [13 14], '-', 'Linewidth', 2, 'Color', Dcol)
    axis([0 6 1 15]);
    %set(gca, 'Yticklabel', '')
    set(gca, 'Xticklabel', '')
    
    subplot(4,6,2:6);
    set(gca, 'ColorOrder', Mycolors);      
    hold on;
    plot(avgjNaEd + avgjNaEf, fcomps); %title('jNaE')
    plot([0 0], [1 15], 'k--');
    
    plot(avgjKEd  + avgjKEf + xsep, fcomps); %title('jKE')
    plot([xsep xsep], [1 15], 'k--');
    
    plot(avgjCaEd + avgjCaEf +2*xsep, fcomps); %title('jCaE')
    plot([2*xsep 2*xsep], [1 15], 'k--');
    
    plot(avgjXEd + avgjXEf + 3*xsep, fcomps); %title('jXE')
    plot([3*xsep 3*xsep], [1 15], 'k--');
    
    plot(avgiE + 4.5*xsep, fcomps); %title('iE')
    plot([4.5*xsep 4.5*xsep], [1 15], 'k--');
    
    legt2 = round(10*t(timepts(2:end)))/10;
    legt1 = [0; legt2(1:end-1)];
    mylegend = [num2str(legt1), ['     ']', ['-----']', ['     ']', num2str(legt2), ['sssss']'];
    
    hmfr = legend(mylegend, 'orientation', 'horizontal');
    set(hmfr, 'Fontsize', 12);
    plot([-1 1],[0.5 0.5], 'k-');
    plot([-1 1],[14.5 14.5], 'k-');    
    text(-2e-14, 11.5, 'field', 'Fontsize', 14)

%scale bar
% For 10 mu mol/m^2/s
%    barlength = 10e-6*Avox
%    plot([4.4*xsep 4.4*xsep+barlength], [2 2] , '-k', 'LineWidth', 3); %scalebar 10 mu mol/m^2/s
%    text(4.4*xsep, 3.5, '10 {\mu}mol/(m^2s)' );

    barlength = 10e-6*Avox
    plot([2.02*xsep 2.02*xsep+barlength], [10 10] , '-k', 'LineWidth', 2.5); %scalebar 10 mu mol/m^2/s
    plot([1.98*xsep-barlength 1.98*xsep], [10 10] , '-k', 'LineWidth', 2.5); %scalebar 10 mu mol/m^2/s

    text(2.01*xsep, 8, '10 {\mu}mol/(m^2s)', 'Fontsize',12);
    text(2.01*xsep + 0.2*barlength, 12, '+', 'Fontsize', 13 );
    text(1.99*xsep-0.7*barlength, 12, '-', 'Fontsize', 13 );
    
    axis(haxis);
    set(gca, 'Yticklabel', '')
    set(gca, 'Xticklabel', '')
    set(gca, 'Xtick',1000)
    set(gca, 'Visible', 'off')
end