%%% THIS CODE CALLS MSO_dae.m %%%

% Computes Vm of MSO neuron model and Ve in extracellular "virtual cylinder"
% Accompanies the manuscript [submitted to J. Neuroscience]:
% "A model of the medial superior olive explains spatiotemporal features of local field potentials"
% JH Goldwyn, M Mc Laughlin, E Verschooten, PX Joris, J Rinzel
%
% MSO model neuron was introduced in:
% "Control of submillisecond synaptic timing in binaural coincidence detectors by Kv1 channels"
% Paul J Mathews, Pablo E Jercog,	John Rinzel, Luisa L Scott, Nace L Golding
% Nature Neuroscience 13 601-609 (2010)

% Simulation code by Joshua H Goldwyn
% Submitted to ModelDB 1/14/13 by Joshua H Goldwyn [jgoldwyn@nyu.edu]

close all
clear all

% Reproduce figures from manuscript [1=Yes, 0=No]
MakeFig4  = 1; 
MakeFig8  = 1;
MakeFig11 = 1;


%%%     MAKE FIGURE 8 -- MONOLATERAL 1kHZ EXCITATION  %%%
if MakeFig4
    
    figure(4), clf
    set(gcf,'position',[1           1        1436         805])
    
    %%% MONOLATERAL 1kHz EXCITATION %%%
    tEnd = 7.;             % simulation duration [ms]
    stimType = 'left';     % monolateral excitation
    gE = 10;               % excitatory conductance [mS / cm2]
    tauE = 0.2;            % excitatory time constant (alpha function) [ms]
    csynE = [2 22];        % location of excitation (compartment number)
    gI = 0;                % inhibitory conductance [mS / cm2]
    tauI = [.4 2];         % inhibitory time constants (double exponential function) [ms]
    csynI = [12];          % location of inhibition (compartment number)
    synFreq = [1000 1001]; % EPSP frequency (Hz) for each dendrite. inhibition freq is first entry
    synDelay = [.0 .0];    % Delay of EPSP onset times in each dendrite [ms]
    inhibDelay = 0;        % Delay of inhibition relative to excitation in first entry of synDelay
    FreezeKLT = 0;         % Whether to Freeze KLT conductance at rest (0=No)
    rB = 11;               % radius of extracellular virtual cylinder (must be larger than soma radius = 10) [micro m]

    % Run model
    out = MSO_dae(tEnd, stimType, gE, tauE, csynE, gI, tauI, csynI, synFreq, synDelay, inhibDelay, FreezeKLT, rB);
    [nt,nx] = size(out.Ve);
    Ve = out.Ve; % Extracellular potential [mV]
    Vm = out.Vm; % Membrane potential [mV]
    Isyn = 1e3*(out.Isyn.*repmat(out.ParamStruct.Surface',nt,1)); % synaptic current [nA]
    Im =1e3* (out.Im.*repmat(out.ParamStruct.Surface',nt,1));  % net membrane current [nA]
    x = out.x; % spatial location of compartments [micro m, 0 is soma center]
    t = out.t; % time (ms)
    ParamStruct = out.ParamStruct; % Input parameters
    dx = ParamStruct.dx*1e4; % distance between compartments [micro m]
    dt = t(2)-t(1); % time step [ms]
    dxG = out.ParamStruct.dxG*1e4; % Distance to ground [microm m]

    % Time and space grid for surface plots of Vm
    [tSurf, xSurf] = meshgrid([t t(end)+dt], [-160 ;-160+cumsum(dx)]); tSurf = tSurf'; xSurf = xSurf';

    % Time and space grid for surface plots of extracellular domain [extended to 0mV at ground via linear decay]
    [tSurfE, xSurfE] = meshgrid([t t(end)+dt], [[-160-dxG:10:-170]' ; -160 ;-160+cumsum(dx) ; [170:10:160+dxG]']); tSurfE = tSurfE'; xSurfE = xSurfE';

    % Extracellular voltage extended to 0mV at ground via linear decay
    for i=1:nt
        VeE(i,:) = [linspace(0,Ve(i,1),100) Ve(i,:) linspace(Ve(i,end),0,100)] ;
    end

    m = mean(VeE(4001:end,:),1); % Mean of "ongoing" response at each spatial location
    VeE0 = VeE -repmat(m,nt,1); % "DC removed" responses

    FontSize = 24;
    
    % Panel A: Vm
    subplot(2,3,1)
    surf(tSurf,xSurf,zeros(size(tSurf)),Vm); shading flat; view([0, 90])
    axis([4 7 min(min(xSurf)) max(max(xSurf))])
    caxis([-80,-40]);
    set(gca,'XTick',0:tEnd,'YTick',[-160 0 160],'FontSize',FontSize)
    cb = colorbar('location','northoutside');
    set(cb,'xlim',[-80 -40],'XTick',[-80:20:-4],'FontSize',FontSize);
    xlabel('Time (ms)', 'FontSize',FontSize)
    ylabel('Distance (\mum)', 'FontSize',FontSize)
    ti = title('Vm (mV)','FontSize',FontSize);
    set(gca,'position',[.13 .6 .2 .2])
    set(ti,'position',get(ti,'position')+[0 160 0])



    % Panel B: Ve 
    subplot(2,3,2)
    surf(tSurfE,xSurfE,zeros(size(tSurfE)),VeE); shading flat; view([0, 90])
    caxis(.325*[-1 1])
    set(gca,'XTick',0:tEnd,'YTick',[-1000 0 1000],'FontSize',FontSize)
    cb = colorbar('location','northoutside');
    set(cb,'xlim',[-.325 .325],'XTick',[-.3 0 .3],'FontSize',FontSize);
    axis([4 7 min(min(xSurfE)) max(max(xSurfE))])
    xlabel('Time (ms)', 'FontSize',FontSize)
    ti = title('Ve (mV)','FontSize',FontSize);
    set(gca,'position',[.43 .6 .2 .2])
    set(ti,'position',get(ti,'position')+[0 1200 0])


    % Panel C -- Ve [DC removed]
    subplot(2,3,3)
    surf(tSurfE,xSurfE,zeros(size(tSurfE)),VeE0); shading flat; view([0, 90])
    caxis(.2*[-1 1])
    set(gca,'XTick',0:tEnd,'YTick',[-1000 0 1000],'FontSize',FontSize)
    cb = colorbar('location','northoutside');
    set(cb,'xlim',[-.2 .2],'XTick',[-.2 0 .2],'FontSize',FontSize);
    axis([4 7 min(min(xSurfE)) max(max(xSurfE))])
    xlabel('Time (ms)', 'FontSize',FontSize)
    ti = title('Ve [No DC] (mV)','FontSize',FontSize);
    set(gca,'position',[.76 .6 .2 .2])
    set(ti,'position',get(ti,'position')+[0 1200 0])



    % Panel D : Membrane current 
    subplot(2,3,4)
    surf(tSurf,xSurf,zeros(size(tSurf)),Im); shading flat; view([0, 90])
    axis([4 7 min(min(xSurf)) max(max(xSurf))])
    caxis([-.75 .75]);
    set(gca,'FontSize',FontSize,'Xtick',0:tEnd)
    cb = colorbar('location','northoutside');
    ti = title('I memb (nA)','FontSize',FontSize);
    set(gca,'YTick',[-160 0 160])
    set(cb,'xlim',[ -.75 .15],'XTick',[ -.6:.3:0],'YTick',[-160 0 160],'FontSize',FontSize);
    set(ti,'position',get(ti,'position')+[0 160 0])
    xlabel('Time (ms)', 'FontSize',FontSize)
    ylabel('Distance (\mum)', 'FontSize',FontSize)
    set(gca,'position',[.13 .1 .2 .2])


    % Panel E: Membrane current without synapse current
    subplot(2,3,5)
    surf(tSurf,xSurf,zeros(size(tSurf)),Im-Isyn); shading flat; view([0, 90])
    axis([4 7 min(min(xSurf)) max(max(xSurf))])
    caxis([-.121,.121]);
    set(gca,'FontSize',FontSize,'Xtick',4:tEnd,'YTick',[-160 0 160])
    xlabel('Time (ms)', 'FontSize',FontSize)
    cb = colorbar('location','northoutside');
    ti = title('I memb w/out  I syn (nA)','FontSize',FontSize);
    set(cb,'xlim',[ -.01 .121],'XTick',[ 0 .05 .1],'FontSize',FontSize);
    set(ti,'position',get(ti,'position')+[0 160 0])
    set(gca,'position',[.43 .1 .2 .2])



    % Panel F: Membrane current by region
    % compartments of cell regions
    iDend1 = 1:10; iSoma = 11:13; iDend2 = 14:23;

    subplot(2,3,6), hold all
    plot(t,sum(Isyn'),'k','linewidth',2)
    plot(t,[sum(Im(:,iDend1)') ;sum(Im(:,iSoma)') ; sum(Im(:,iDend2)')],'linewidth',2)
    leg = legend({'Synapse','Near Dend.','Soma','Far Dend.'},'location','northoutside','fontsize',16);
    set(gca,'XTick',0:tEnd,'YTick',-.8:.4:.4,'FontSize',FontSize)
    xlabel('Time (ms)','FontSize',FontSize)
    ylabel('Current (nA)','FontSize',FontSize)
    axis([4 7 -.85 .4])
    legend('boxoff')
    set(gca,'position',[.76 .1 .2 .2])


end



%%%     MAKE FIGURE 8 -- BILATERAL 1kHz EXCITATION  %%%
if MakeFig8
    
    figure(8)
    set(gcf,'position',[1           1        1436         805])
    
    %%% Out of phase (.5ms delay) %%%
    tEnd = 7.;             % simulation duration [ms]
    stimType = 'both';     % bilateral excitation
    gE = 10;               % excitatory conductance [mS / cm2]
    tauE = 0.2;            % excitatory time constant (alpha function) [ms]
    csynE = [2 22];        % location of excitation (compartment number)
    gI = 0;                % inhibitory conductance [mS / cm2]
    tauI = [.4 2];         % inhibitory time constants (double exponential function) [ms]
    csynI = [12];          % location of inhibition (compartment number)
    synFreq = [1000 1000]; % EPSP frequency (Hz) for each dendrite. inhibition freq is first entry
    synDelay = [0  .5];    % Delay of EPSP onset times in each dendrite [ms]
    inhibDelay = 0;        % Delay of inhibition relative to excitation in first entry of synDelay
    FreezeKLT = 0;         % Whether to Freeze KLT conductance at rest (0=No)
    rB = 11;               % radius of extracellular virtual cylinder (must be larger than soma radius = 10) [micro m]

    % Run model
    out = MSO_dae(tEnd, stimType, gE, tauE, csynE, gI, tauI, csynI, synFreq, synDelay, inhibDelay, FreezeKLT, rB);
    [nt,nx] = size(out.Ve);
    Ve = out.Ve; % Extracellular potential [mV]
    Vm = out.Vm; % Membrane potential [mV]
    Isyn = 1e3*(out.Isyn.*repmat(out.ParamStruct.Surface',nt,1)); % synaptic current [nA]
    Im =1e3* (out.Im.*repmat(out.ParamStruct.Surface',nt,1));  % net membrane current [nA]
    x = out.x; % spatial location of compartments [micro m, 0 is soma center]
    t = out.t; % time (ms)
    ParamStruct = out.ParamStruct; % Input parameters
    dx = ParamStruct.dx*1e4; % distance between compartments [micro m]
    dt = t(2)-t(1); % time step [ms]
    dxG = out.ParamStruct.dxG*1e4; % Distance to ground [microm m]

    % Time and space grid for surface plots of Vm
    [tSurf, xSurf] = meshgrid([t t(end)+dt], [-160 ;-160+cumsum(dx)]); tSurf = tSurf'; xSurf = xSurf';

    % Time and space grid for surface plots of extracellular domain [extended to 0mV at ground via linear decay]
    [tSurfE, xSurfE] = meshgrid([t t(end)+dt], [[-160-dxG:10:-170]' ; -160 ;-160+cumsum(dx) ; [170:10:160+dxG]']); tSurfE = tSurfE'; xSurfE = xSurfE';

    % Extracellular voltage extended to 0mV at ground via linear decay
    for i=1:nt
        VeE(i,:) = [linspace(0,Ve(i,1),100) Ve(i,:) linspace(Ve(i,end),0,100)] ;
    end

    m = mean(VeE(4001:end,:),1); % Mean of "ongoing" response at each spatial location
    VeE0 = VeE -repmat(m,nt,1); % "DC removed" responses

    FontSize = 24;
    
    % Panel A: Vm
    subplot(2,3,1)
    surf(tSurf,xSurf,zeros(size(tSurf)),Vm); shading flat; view([0, 90])
    axis([4 7 min(min(xSurf)) max(max(xSurf))])
    caxis([-80,-40]);
    set(gca,'XTick',0:tEnd,'YTick',[-160 0 160],'FontSize',FontSize)
    cb = colorbar('location','northoutside');
    set(cb,'xlim',[-80 -40],'XTick',[-80:20:-4],'FontSize',FontSize);
    xlabel('Time (ms)', 'FontSize',FontSize)
    ylabel('Distance (\mum)', 'FontSize',FontSize)
    ti = title('Vm (mV)','FontSize',FontSize);
    set(gca,'position',[.13 .6 .2 .2])
    set(ti,'position',get(ti,'position')+[0 160 0])



    % Panel B: Ve 
    subplot(2,3,2)
    surf(tSurfE,xSurfE,zeros(size(tSurfE)),VeE); shading flat; view([0, 90])
    caxis(.35*[-1 1])
    set(gca,'XTick',0:tEnd,'YTick',[-1000 0 1000],'FontSize',FontSize)
    cb = colorbar('location','northoutside');
    set(cb,'xlim',[-.35 .35],'XTick',[-.3 0 .3],'FontSize',FontSize);
    axis([4 7 min(min(xSurfE)) max(max(xSurfE))])
    xlabel('Time (ms)', 'FontSize',FontSize)
    ti = title('Ve (mV)','FontSize',FontSize);
    set(gca,'position',[.43 .6 .2 .2])
    set(ti,'position',get(ti,'position')+[0 1200 0])


    % Panel C -- Ve [DC removed]
    subplot(2,3,3)
    surf(tSurfE,xSurfE,zeros(size(tSurfE)),VeE0); shading flat; view([0, 90])
    caxis([-.21 .21])
    set(gca,'XTick',0:tEnd,'YTick',[-1000 0 1000],'FontSize',FontSize)
    cb = colorbar('location','northoutside');
    set(cb,'xlim',[-.21 .21],'XTick',[-.2 0 .2],'FontSize',FontSize);
    axis([4 7 min(min(xSurfE)) max(max(xSurfE))])
    xlabel('Time (ms)', 'FontSize',FontSize)
    ti = title('Ve [No DC] (mV)','FontSize',FontSize);
    set(gca,'position',[.76 .6 .2 .2])
    set(ti,'position',get(ti,'position')+[0 1200 0])



    %%% In phase (0ms delay) %%%
    tEnd = 7.;             % simulation duration [ms]
    stimType = 'both';     % bilateral excitation
    gE = 10;               % excitatory conductance [mS / cm2]
    tauE = 0.2;            % excitatory time constant (alpha function) [ms]
    csynE = [2 22];        % location of excitation (compartment number)
    gI = 0;                % inhibitory conductance [mS / cm2]
    tauI = [.4 2];         % inhibitory time constants (double exponential function) [ms]
    csynI = [12];          % location of inhibition (compartment number)
    synFreq = [1000 1000]; % EPSP frequency (Hz) for each dendrite. inhibition freq is first entry
    synDelay = [0 .0];    % Delay of EPSP onset times in each dendrite [ms]
    inhibDelay = 0;        % Delay of inhibition relative to excitation in first entry of synDelay
    FreezeKLT = 0;         % Whether to Freeze KLT conductance at rest (0=No)
    rB = 11;               % radius of extracellular virtual cylinder (must be larger than soma radius = 10) [micro m]

    % Run model
    out = MSO_dae(tEnd, stimType, gE, tauE, csynE, gI, tauI, csynI, synFreq, synDelay, inhibDelay, FreezeKLT, rB);
    [nt,nx] = size(out.Ve);
    Ve = out.Ve; % Extracellular potential [mV]
    Vm = out.Vm; % Membrane potential [mV]
    Isyn = 1e3*(out.Isyn.*repmat(out.ParamStruct.Surface',nt,1)); % synaptic current [nA]
    Im =1e3* (out.Im.*repmat(out.ParamStruct.Surface',nt,1));  % net membrane current [nA]
    x = out.x; % spatial location of compartments [micro m, 0 is soma center]
    t = out.t; % time (ms)
    ParamStruct = out.ParamStruct; % Input parameters
    dx = ParamStruct.dx*1e4; % distance between compartments [micro m]
    dt = t(2)-t(1); % time step [ms]
    dxG = out.ParamStruct.dxG*1e4; % Distance to ground [microm m]

    % Time and space grid for surface plots of Vm
    [tSurf, xSurf] = meshgrid([t t(end)+dt], [-160 ;-160+cumsum(dx)]); tSurf = tSurf'; xSurf = xSurf';

    % Time and space grid for surface plots of extracellular domain [extended to 0mV at ground via linear decay]
    [tSurfE, xSurfE] = meshgrid([t t(end)+dt], [[-160-dxG:10:-170]' ; -160 ;-160+cumsum(dx) ; [170:10:160+dxG]']); tSurfE = tSurfE'; xSurfE = xSurfE';

    % Extracellular voltage extended to 0mV at ground via linear decay
    for i=1:nt
        VeE(i,:) = [linspace(0,Ve(i,1),100) Ve(i,:) linspace(Ve(i,end),0,100)] ;
    end

    m = mean(VeE(4001:end,:),1); % Mean of "ongoing" response at each spatial location
    VeE0 = VeE -repmat(m,nt,1); % "DC removed" responses

    %Panel D: Vm
    subplot(2,3,4)
    surf(tSurf,xSurf,zeros(size(tSurf)),Vm); shading flat; view([0, 90])
    axis([4 7 min(min(xSurf)) max(max(xSurf))])
    caxis([-80,-40]);
    set(gca,'XTick',0:tEnd,'YTick',[-160 0 160],'FontSize',FontSize)
    cb = colorbar('location','northoutside');
    set(cb,'xlim',[-80 -40],'XTick',[-80:20:-4],'FontSize',FontSize);
    xlabel('Time (ms)', 'FontSize',FontSize)
    ylabel('Distance (\mum)', 'FontSize',FontSize)
    ti = title('Vm (mV)','FontSize',FontSize);
    set(gca,'position',[.13 .1 .2 .2])
    set(ti,'position',get(ti,'position')+[0 160 0])



    % Panel E: Ve 
    subplot(2,3,5)
    surf(tSurfE,xSurfE,zeros(size(tSurfE)),VeE); shading flat; view([0, 90])
    caxis([-.5 .5])
    set(gca,'XTick',0:tEnd,'YTick',[-1000 0 1000],'FontSize',FontSize)
    cb = colorbar('location','northoutside');
    set(cb,'xlim',[-.5 .5],'XTick',[-.5 0 .5],'FontSize',FontSize);
    axis([4 7 min(min(xSurfE)) max(max(xSurfE))])
    xlabel('Time (ms)', 'FontSize',FontSize)
    ti = title('Ve (mV)','FontSize',FontSize);
    set(gca,'position',[.43 .1 .2 .2])
    set(ti,'position',get(ti,'position')+[0 1200 0])


    % Panel F -- Ve [DC removed]
    subplot(2,3,6)
    surf(tSurfE,xSurfE,zeros(size(tSurfE)),VeE0); shading flat; view([0, 90])
    caxis([-.21 .21])
    set(gca,'XTick',0:tEnd,'YTick',[-1000 0 1000],'FontSize',FontSize)
    cb = colorbar('location','northoutside');
    set(cb,'xlim',[-.21 .21],'XTick',[-.2 0 .2],'FontSize',FontSize);
    axis([4 7 min(min(xSurfE)) max(max(xSurfE))])
    xlabel('Time (ms)', 'FontSize',FontSize)
    ti = title('Ve [No DC] (mV)','FontSize',FontSize);
    set(gca,'position',[.76 .1 .2 .2])
    set(ti,'position',get(ti,'position')+[0 1200 0])


end


%%%     MAKE FIGURE 11 -- MONOLATERAL 1kHz EXCITATION + SOMA INHIBIITON %%%
if MakeFig11
    
    figure(11), clf
    set(gcf,'position',[1           300        1436         400])
    
    %%% MONOLATERAL 1kHz EXCITATION with soma inhibition %%%
    tEnd = 7.;             % simulation duration [ms]
    stimType = 'left';     % monolateral excitation
    gE = 10;               % excitatory conductance [mS / cm2]
    tauE = 0.2;            % excitatory time constant (alpha function) [ms]
    csynE = [2 22];        % location of excitation (compartment number)
    gI = 4;                % inhibitory conductance [mS / cm2]
    tauI = [.4 2];         % inhibitory time constants (double exponential function) [ms]
    csynI = [12];          % location of inhibition (compartment number)
    synFreq = [1000 1001]; % EPSP frequency (Hz) for each dendrite. inhibition freq is first entry
    synDelay = [.0 .0];    % Delay of EPSP onset times in each dendrite [ms]
    inhibDelay = -.35;        % Delay of inhibition relative to excitation in first entry of synDelay
    FreezeKLT = 0;         % Whether to Freeze KLT conductance at rest (0=No)
    rB = 11;               % radius of extracellular virtual cylinder (must be larger than soma radius = 10) [micro m]

    % Run model
    out = MSO_dae(tEnd, stimType, gE, tauE, csynE, gI, tauI, csynI, synFreq, synDelay, inhibDelay, FreezeKLT, rB);
    [nt,nx] = size(out.Ve);
    Ve = out.Ve; % Extracellular potential [mV]
    Vm = out.Vm; % Membrane potential [mV]
    Isyn = 1e3*(out.Isyn.*repmat(out.ParamStruct.Surface',nt,1)); % synaptic current [nA]
    Im =1e3* (out.Im.*repmat(out.ParamStruct.Surface',nt,1));  % net membrane current [nA]
    x = out.x; % spatial location of compartments [micro m, 0 is soma center]
    t = out.t; % time (ms)
    ParamStruct = out.ParamStruct; % Input parameters
    dx = ParamStruct.dx*1e4; % distance between compartments [micro m]
    dt = t(2)-t(1); % time step [ms]
    dxG = out.ParamStruct.dxG*1e4; % Distance to ground [microm m]

    % Time and space grid for surface plots of Vm
    [tSurf, xSurf] = meshgrid([t t(end)+dt], [-160 ;-160+cumsum(dx)]); tSurf = tSurf'; xSurf = xSurf';

    % Time and space grid for surface plots of extracellular domain [extended to 0mV at ground via linear decay]
    [tSurfE, xSurfE] = meshgrid([t t(end)+dt], [[-160-dxG:10:-170]' ; -160 ;-160+cumsum(dx) ; [170:10:160+dxG]']); tSurfE = tSurfE'; xSurfE = xSurfE';

    % Extracellular voltage extended to 0mV at ground via linear decay
    for i=1:nt
        VeE(i,:) = [linspace(0,Ve(i,1),100) Ve(i,:) linspace(Ve(i,end),0,100)] ;
    end
    VeE_withI = VeE;
    
    m = mean(VeE(4001:end,:),1); % Mean of "ongoing" response at each spatial location
    VeE0 = VeE -repmat(m,nt,1); % "DC removed" responses

    FontSize = 24;
    
    % Panel A: Vm
    subplot(1,4,1)
    surf(tSurf,xSurf,zeros(size(tSurf)),Vm); shading flat; view([0, 90])
    axis([4 7 min(min(xSurf)) max(max(xSurf))])
    caxis([-80,-40]);
    set(gca,'XTick',0:tEnd,'YTick',[-160 0 160],'FontSize',FontSize)
    cb = colorbar('location','northoutside');
    set(cb,'xlim',[-80 -40],'XTick',[-80:20:-4],'FontSize',FontSize);
    xlabel('Time (ms)', 'FontSize',FontSize)
    ylabel('Distance (\mum)', 'FontSize',FontSize)
    ti = title('Vm (mV)','FontSize',FontSize);
    set(gca,'position',[.08 .2 .16 .5])
    set(ti,'position',get(ti,'position')+[0 100 0])
    
    % Panel B: Ve 
    subplot(1,4,2)
    surf(tSurfE,xSurfE,zeros(size(tSurfE)),VeE); shading flat; view([0, 90])
    caxis([-.61 .61])
    set(gca,'XTick',0:tEnd,'YTick',[-1000 0 1000],'FontSize',FontSize)
    cb = colorbar('location','northoutside');
    set(cb,'xlim',[-.61 .61],'XTick',[-.6 0 .6],'FontSize',FontSize);
    axis([4 7 min(min(xSurfE)) max(max(xSurfE))])
    xlabel('Time (ms)', 'FontSize',FontSize)
    ti = title('Ve (mV)','FontSize',FontSize);
    set(gca,'position',[.32 .2 .16 .5])
    set(ti,'position',get(ti,'position')+[0 700 0])


    % Panel C -- Ve [DC removed]
    subplot(1,4,3)
    surf(tSurfE,xSurfE,zeros(size(tSurfE)),VeE0); shading flat; view([0, 90])
    caxis([-.2 .2])
    set(gca,'XTick',0:tEnd,'YTick',[-1000 0 1000],'FontSize',FontSize)
    cb = colorbar('location','northoutside');
    set(cb,'xlim',[-.2 .2],'XTick',[-.2 0 .2],'FontSize',FontSize);
    axis([4 7 min(min(xSurfE)) max(max(xSurfE))])
    xlabel('Time (ms)', 'FontSize',FontSize)
    ti = title('Ve [No DC] (mV)','FontSize',FontSize);
    set(gca,'position',[.56 .2 .16 .5])
    set(ti,'position',get(ti,'position')+[0 700 0])
    
    
    %%% MONOLATERAL 1kHz EXCITATION without soma inhibition %%%
    tEnd = 7.;             % simulation duration [ms]
    stimType = 'left';     % monolateral excitation
    gE = 10;               % excitatory conductance [mS / cm2]
    tauE = 0.2;            % excitatory time constant (alpha function) [ms]
    csynE = [2 22];        % location of excitation (compartment number)
    gI = 0;                % inhibitory conductance [mS / cm2]
    tauI = [.4 2];         % inhibitory time constants (double exponential function) [ms]
    csynI = [12];          % location of inhibition (compartment number)
    synFreq = [1000 1001]; % EPSP frequency (Hz) for each dendrite. inhibition freq is first entry
    synDelay = [.0 .0];    % Delay of EPSP onset times in each dendrite [ms]
    inhibDelay = -.35;        % Delay of inhibition relative to excitation in first entry of synDelay
    FreezeKLT = 0;         % Whether to Freeze KLT conductance at rest (0=No)
    rB = 11;               % radius of extracellular virtual cylinder (must be larger than soma radius = 10) [micro m]

    % Run model
    out = MSO_dae(tEnd, stimType, gE, tauE, csynE, gI, tauI, csynI, synFreq, synDelay, inhibDelay, FreezeKLT, rB);
    [nt,nx] = size(out.Ve);
    Ve = out.Ve; % Extracellular potential [mV]
    Vm = out.Vm; % Membrane potential [mV]
    Isyn = 1e3*(out.Isyn.*repmat(out.ParamStruct.Surface',nt,1)); % synaptic current [nA]
    Im =1e3* (out.Im.*repmat(out.ParamStruct.Surface',nt,1));  % net membrane current [nA]
    x = out.x; % spatial location of compartments [micro m, 0 is soma center]
    t = out.t; % time (ms)
    ParamStruct = out.ParamStruct; % Input parameters
    dx = ParamStruct.dx*1e4; % distance between compartments [micro m]
    dt = t(2)-t(1); % time step [ms]
    dxG = out.ParamStruct.dxG*1e4; % Distance to ground [microm m]

    % Time and space grid for surface plots of Vm
    [tSurf, xSurf] = meshgrid([t t(end)+dt], [-160 ;-160+cumsum(dx)]); tSurf = tSurf'; xSurf = xSurf';

    % Time and space grid for surface plots of extracellular domain [extended to 0mV at ground via linear decay]
    [tSurfE, xSurfE] = meshgrid([t t(end)+dt], [[-160-dxG:10:-170]' ; -160 ;-160+cumsum(dx) ; [170:10:160+dxG]']); tSurfE = tSurfE'; xSurfE = xSurfE';

    % Extracellular voltage extended to 0mV at ground via linear decay
    for i=1:nt
        VeE(i,:) = [linspace(0,Ve(i,1),100) Ve(i,:) linspace(Ve(i,end),0,100)] ;
    end
    

    % Panel D -- [Ve w/inhibition - Ve w/out inhibition]
    subplot(1,4,4)
    surf(tSurfE,xSurfE,zeros(size(tSurfE)),VeE_withI - VeE); shading flat; view([0, 90])
    caxis([-.34 .34])
    set(gca,'XTick',0:tEnd,'YTick',[-1000 0 1000],'FontSize',FontSize)
    cb = colorbar('location','northoutside');
    set(cb,'xlim',[-.34 .34],'XTick',[-.3 0 .3],'FontSize',FontSize);
    axis([4 7 min(min(xSurfE)) max(max(xSurfE))])
    xlabel('Time (ms)', 'FontSize',FontSize)
    ti = title('Ve Difference (mV)','FontSize',FontSize);
    set(gca,'position',[.8 .2 .16 .5])
    set(ti,'position',get(ti,'position')+[0 700 0])

    
end