% CODE TO PRODUCE PANEL B of two compartment LSO model in 
% "Glycinergic axonal inhibition subserves acute spatial sensitivity to sudden increases in sound intensity"
% Tom P. Franken, Brian J. Bondy, David B. Haimes, Joshua H. Goldwyn, Nace L. Golding, Philip H. Smith, Philip X. Joris

close all
clear all

    
% coupling coefficiencts
stoa = .95; atos = .60;

POS = [0.13    0.84    0.83    0.1;
       0.13    0.65    0.83    0.1;
       0.13    0.47    0.83    0.1;
       0.13    0.29    0.83    0.1;
       0.13    0.1     0.83    0.1];
       
figure(1), clf, 
    dataSet = 'dclamp'; % PSG kinetics from in vitro study
    ITDplot = [0:.6:2.4];
    areaRatio = .12; KHTfrac = .1; KLTfrac = [0 0];
    uE = 2.3; nE = 20; % excitatory inputs
    uI = 3.1; nI = 8; % inhibitory inputs
    P = getParam2cptWC([stoa, atos],areaRatio,'yes',[],KHTfrac,KLTfrac,dataSet);
    P.tE0 = 3; P.tI1 = 5; P.tI2 = P.tI1;
    P.gE = 0; P.gI1 = uI*(nI-2); P.gI2 = 2*uI;
    P.I1 = 0; P.I2 = 0;
    [ti,x] = ode15s(@twoCptODEWC ,[0 20], P.x0 ,P.options,P);
    ipsp = x(:,1) - P.vRest;

    for i=1:length(ITDplot)
        subplot('position',POS(i,:))
        hold all
        P.tE0 = 5+ITDplot(i) ; P.tI1 = 5; P.tI2 = P.tI1; 
        P.gE = nE*uE; P.gI1 = uI*(nI-2); P.gI2 = 2*uI;
        P.I1 = 0; P.I2 = 0;
        [t,x] = ode15s(@twoCptODEWC ,[0 20], P.x0 ,P.options,P);
        g2=plot(ti-P.tI1,ipsp,'color',.6*[1 1 1],'linewidth',2); 
        g1=plot(t-P.tI1,x(:,1)-P.vRest,'k','linewidth',2);
        axis([-1.5 5.5 -10 33])
        title(['B',num2str(i)],'position',[-1.5 39 0])
        set(gca,'xtick',0:2:4,'ytick',0:30:40)

    end
    
    subplot('position',POS(1,:))
        g3=plot(.0,21,'>r','markerfacecolor','r');
        [hh,icons,plots,txt] = legend([g1 g2], {'E + I', 'I only'},'box','off','position',[.6 .87 .1 .1]);
        p1 = icons(1).Position;
        p2 = icons(2).Position;
        icons(1).Position = [0.15 .4 0];
        icons(2).Position = [0.15 .15  0];
        icons(3).XData = [-0.15 0.1];
        icons(3).YData = icons(1).Position(2)*[1 1];
        icons(5).XData = [-0.15 0.1];
        icons(5).YData = icons(2).Position(2)*[1 1];
        text(3,35,['ITD=',num2str(ITDplot(1)),'ms']);

    subplot('position',POS(2,:))
        text(2.4,25,['ITD=',num2str(ITDplot(2)),'ms']);
        
    subplot('position',POS(3,:))
        text(2.4,25,['ITD=',num2str(ITDplot(3)),'ms']);
        ylabel('soma voltage (mV)')

    subplot('position',POS(4,:))
        g3=plot(2.15,21,'>r','markerfacecolor','r');
        text(-1.4,25,['ITD=',num2str(ITDplot(4)),'ms']);

    subplot('position',POS(5,:))
        g3=plot(2.55,21,'>r','markerfacecolor','r');
        text(-1.4,25,['ITD=',num2str(ITDplot(5)),'ms']);
        xlabel('time (ms)')

set(gcf,'units','inches','position',[0 0 3.5 5])
set(gcf, 'PaperPositionMode','auto')



% PSG kinetics from Beiderbeck et al 2018
figure(2), clf, 
    dataSet = 'beiderbeck'; % PSG kinetics from Beiderbeck et al 2018
    ITDplot = [-.3:.3:0.9];
    areaRatio = .12; KHTfrac = .1; KLTfrac = [0 0];
    uE = 1.6; nE = 20; % excitatory inputs
    uI = 3.9; nI = 8; % inhibitory inputs
    P = getParam2cptWC([stoa, atos],areaRatio,'yes',[],KHTfrac,KLTfrac,dataSet);
    P.tE0 = 3; P.tI1 = 5; P.tI2 = P.tI1;
    P.gE = 0; P.gI1 = uI*(nI-2); P.gI2 = 2*uI;
    P.I1 = 0; P.I2 = 0;
    [ti,x] = ode15s(@twoCptODEWC ,[0 20], P.x0 ,P.options,P);
    ipsp = x(:,1) - P.vRest;

    for i=1:length(ITDplot)
        subplot('position',POS(i,:))
        hold all
        P.tE0 = 5+ITDplot(i) ; P.tI1 = 5; P.tI2 = P.tI1; 
        P.gE = nE*uE; P.gI1 = uI*(nI-2); P.gI2 = 2*uI;
        P.I1 = 0; P.I2 = 0;
        [t,x] = ode15s(@twoCptODEWC ,[0 20], P.x0 ,P.options,P);
        g2=plot(ti-P.tI1,ipsp,'color',.6*[1 1 1],'linewidth',2); 
        g1=plot(t-P.tI1,x(:,1)-P.vRest,'k','linewidth',2);
        axis([-1.5 5.5 -10 37])
        title(['B',num2str(i)],'position',[-1.5 41 0])
        set(gca,'xtick',0:2:4,'ytick',0:30:40)

    end
    
    subplot('position',POS(1,:))
        g3=plot(-.1,21,'>r','markerfacecolor','r');
        [hh,icons,plots,txt] = legend([g1 g2], {'E + I', 'I only'},'box','off','position',[.6 .87 .1 .1]);
        p1 = icons(1).Position;
        p2 = icons(2).Position;
        icons(1).Position = [0.15 .4 0];
        icons(2).Position = [0.15 .15  0];
        icons(3).XData = [-0.15 0.1];
        icons(3).YData = icons(1).Position(2)*[1 1];
        icons(5).XData = [-0.15 0.1];
        icons(5).YData = icons(2).Position(2)*[1 1];
        text(3,37,['ITD=',num2str(ITDplot(1)),'ms']);

    subplot('position',POS(2,:))
        text(2.4,29,['ITD=',num2str(ITDplot(2)),'ms']);
        
    subplot('position',POS(3,:))
        text(2.4,29,['ITD=',num2str(ITDplot(3)),'ms']);
        ylabel('soma voltage (mV)')

    subplot('position',POS(4,:))
        g3=plot(1.12,21,'>r','markerfacecolor','r');
        text(2.4,29,['ITD=',num2str(ITDplot(4)),'ms']);

    subplot('position',POS(5,:))
        g3=plot(1.3,21,'>r','markerfacecolor','r');
        text(2.4,29,['ITD=',num2str(ITDplot(5)),'ms']);
        xlabel('time (ms)')

set(gcf,'units','inches','position',[3.7 0 3.5 5])
set(gcf, 'PaperPositionMode','auto')




function P = getParam2cptWC(couple,areaRatio,getC,gNa2,KHTfrac,KLTfrac,dataSet)

    % Outputs %
    % P parameter structure [defined at bottom of function]

    P.options = odeset('RelTol',1e-8,'AbsTol',1e-8); %[];%odeset('maxstep',.01);
    
    % Wang Colburn (Ashida 2017) 
    vShift = +5; % 0; % +5mV for adjusted. 0 for normal
    P.winf = @(v) (1 ./ (1 + exp(-(v+-vShift+48)/6) ) ).^(1/4);
    P.tauw = @(v) 1.5 + 100./(6*exp(v-vShift+60)+16*exp(-(v-vShift+60)/45)) ;
    P.zinf = @(v) 0.5 + 0.5./(1+exp((v-vShift+71)/10) );
    P.tauz = @(v) 50+ 1000./( exp( (v-vShift+60)/20) + exp(-(v-vShift+60)/8) );
    P.ninf = @(v) (1./ (1+exp(-(v-vShift+15)/5)) ).^(1/2);
    P.taun = @(v) 0.7 + 100./(11*exp((v-vShift+60)/24) + 21*exp(-(v-vShift+60)/23) );
    P.pinf = @(v) 1./ (1+exp(-(v-vShift+23)/6) ); 
    P.taup = @(v) 5 + 100./(4*exp((v-vShift+60)/32) + 5*exp(-(v-vShift+60)/22));
    P.minf = @(v) 1./ (1+exp(-(v-vShift+38)/7) ); 
    P.taum = @(v) (0.04 + 10./(5*exp( (v-vShift+60)/18) + 36*exp(-(v-vShift+60)/25)));
    P.hinf = @(v) 1./ (1+exp((v-vShift+65)/6));
    P.tauh = @(v) 0.6 + 100./(7*exp((v-vShift+60)/11) + 10*exp(-(v-vShift+60)/25));

    % Fixed parameters %
    P.R1      = 40 * 1e-3; % Input resistance to CPT1 [10^9 Ohm]
    P.tauEst  = 1;         % "estimated time constant" [ms]
    P.vRest   = -60;       % resting membrane potential [mV]
    P.EsynE   = 0;  
    P.EsynI   = -75;
    P.EK      = -70;
    P.ENa     = 50;
    Q10       = 3;
    T         = 37;
    P.phi     = Q10^((T-22)/10);
   
    P.vThresh  = -20;
    
    % Initial values %
    v0 = P.vRest;  w0 = P.winf(P.vRest); z0 = P.zinf(P.vRest); m0 = P.minf(P.vRest); h0 = P.hinf(P.vRest); n0 = P.ninf(P.vRest); p0 = P.pinf(P.vRest);
    P.x0 = [v0;v0;w0;z0;w0;z0;n0;p0;m0;h0];

    % SYNAPTIC INPTUS
    % gE threshold for finding gNa
%     P.gEthresh = 17.5; % nS.   OLD. gives 15mV EPSP (from Ashida)
    
	switch(dataSet)
        case('beiderbeck') % Beiderbeck.. .Grothe 
                P.gEthresh = 11; % nS.  gives 10mV EPSP 
                P.epsg = @(t) (1/0.17133) .* ones(size(t)) .* (t>0) .* (exp(-t/.4) - exp(-t/.25) ); 
                P.ipsg = @(t) (1/0.25) .* ones(size(t)) .* (t>0) .* (exp(-t/.7) - exp(-t/.35) ); 

        case('dclamp') % original [same as in vitro]
                P.gEthresh = 19; % nS.  gives 10mV EPSP 
                P.epsg = @(t) (1/0.21317) .* ones(size(t)) .* (t>0) .* (exp(-t/.18) - exp(-t/.1) ); 
                P.ipsg = @(t) (1/0.5025) .* ones(size(t)) .* (t>0) .* (exp(-t/2) - exp(-t/.45) ); 

    end

    % Passive parameters %    
    P.couple12 = couple(1); P.couple21 = couple(2); 
    
    P.gC    = (1/P.R1) * P.couple21 / (1-P.couple12*P.couple21); % coupling conductance [nS]
    P.gTot1  = P.gC * (1/P.couple21 - 1); % CPT1 total conductance [nS]
    P.gTot2  = P.gC * (1/P.couple12 - 1); % CPT2 total conductance [nS]

    if strcmp(getC,'yes')
        
        % minimize L2 error with exponential
        
        %%%% NUMERICAL EIGEN-CALCULATIONS FOR GENERAL CASE %%%%
        M = @(c1) [ [P.gTot1+P.gC ,  -P.gC]/c1 ; [-P.gC , P.gTot2 + P.gC]/(areaRatio*c1)]; % steady state matrix;
        
        % eigen-calculations. eigenvectors both have 1 in first entry
        eval1 = @(c1) [1,0]*eig(M(c1),'vector'); 
        eval2 = @(c1) [0,1]*eig(M(c1),'vector'); 
        evec1 = @(c1) [1 ; -([1 0] *([M(c1) - eval1(c1)*eye(2)])*[1;0]) / (([1 0] *[M(c1) - eval1(c1)*eye(2)])*[0;1])];
        evec2 = @(c1) [1 ; -([1 0] *([M(c1) - eval2(c1)*eye(2)])*[1;0]) / (([1 0] *[M(c1) - eval2(c1)*eye(2)])*[0;1])];
        
        % coordinates of initial values [1;stoa] in eigenbasis 
        w = @(c1) [evec1(c1) , evec2(c1)] \ [1;couple(1)];
        
        % minimize L2 error with numerical integration
        fobj = @(c1) integral( @(t) (exp(-t/P.tauEst) -  exp(-eval1(c1)*t)*([1,0]*w(c1)) - exp(-eval2(c1)*t)*([0,1]*w(c1)) ).^2, 0,5 ) ;
        P.C1 = fminbnd( @(c1) fobj(c1), 0,100);



    else
        P.C1 = 25;
    end
    
    P.C2  = areaRatio*P.C1; %P.tau2 * (P.gTot2 + P.gC); % CPT2 capacitance [pF]

    P.gKLT1 = P.gTot1*KLTfrac(1)/ (P.winf(P.vRest)^4*P.zinf(P.vRest) );
    P.gL1 = P.gTot1- P.gKLT1*(P.winf(P.vRest)^4*P.zinf(P.vRest) );
    P.EL1 = P.vRest-(P.gKLT1*(P.winf(P.vRest)^4*P.zinf(P.vRest) ).*(P.EK-P.vRest) )/P.gL1; %P.EL   = -56;   % mV

    P.gKLT2 = P.gTot2*KLTfrac(2)/ (P.winf(P.vRest)^4*P.zinf(P.vRest) );
    
    if isempty(gNa2)

        % keep original values; 
        Q = P;
        % find gNa with bisection
        gNaLow = 0; gNaHigh = 20000; P.gNa2 = (gNaHigh + gNaLow) / 2;
        for ib=1:20
           P.gKHT2 = KHTfrac*P.gNa2;  
           P.gL2 = P.gTot2 - P.gKLT2*(P.winf(P.vRest)^4*P.zinf(P.vRest) )  - P.gKHT2*(.85*P.ninf(P.vRest)^2 + .15*P.pinf(P.vRest)) - P.gNa2*P.minf(P.vRest)^3*P.hinf(P.vRest); 
           P.EL2 = P.vRest-(P.gKLT2*P.winf(P.vRest)^4*P.zinf(P.vRest).*(P.EK-P.vRest)+ P.gKHT2*(.85*P.ninf(P.vRest)^2 + .15*P.pinf(P.vRest)).*(P.EK-P.vRest) + P.gNa2*P.minf(P.vRest)^3*P.hinf(P.vRest).*(P.ENa-P.vRest) )/P.gL2;
           P.gE = P.gEthresh; P.gI1 = 0; P.gI2 =0; P.I1 = 0; P.I2=0; P.tE0=2; P.tI1=0; P.tI2=0;
           [t,x] = ode15s(@twoCptODEWC ,[0 20], P.x0 ,P.options,P); 
           if max(x(:,2))> P.vThresh; gNaHigh = P.gNa2; % spike, decrease gNa
           elseif max(x(:,2))<P.vThresh; gNaLow = P.gNa2; % no spike, increase gNa
           elseif max(x(:,2))==P.vThresh; break; % exactly at spike threshold
           end
           P.gNa2 = (gNaHigh + gNaLow) / 2;

        end
        % restore values
        if isfield(Q,'gE');  P.gE = Q.gE; end
        if isfield(Q,'gI1'); P.gI1 = Q.gI1; end
        if isfield(Q,'gI2'); P.gI2 =Q.gI2; end
        if isfield(Q,'I1');  P.I1 = Q.I1; end
        if isfield(Q,'I2');  P.I2 = Q.I2; end
        if isfield(Q,'tE0');  P.I2 = Q.tE0; end
        if isfield(Q,'tI1');  P.I2 = Q.tI1; end
        if isfield(Q,'tI2');  P.I2 = Q.tI2; end
    else
        % use gNa given as input
        P.gNa2 = gNa2;
        P.gKHT2 = KHTfrac*P.gNa2; 
        P.gL2 = P.gTot2 - P.gKLT2*(P.winf(P.vRest)^4*P.zinf(P.vRest) )  - P.gKHT2*(.85*P.ninf(P.vRest)^2 + .15*P.pinf(P.vRest)) - P.gNa2*P.minf(P.vRest)^3*P.hinf(P.vRest); 
        P.EL2 = P.vRest-(P.gKLT2*P.winf(P.vRest)^4*P.zinf(P.vRest).*(P.EK-P.vRest)+ P.gKHT2*(.85*P.ninf(P.vRest)^2 + .15*P.pinf(P.vRest)).*(P.EK-P.vRest) + P.gNa2*P.minf(P.vRest)^3*P.hinf(P.vRest).*(P.ENa-P.vRest) )/P.gL2;

    end

end



function dx = twoCptODEWC(t,x,P)

    % read in variables
    v1 = x(1); w1 = x(3); z1 = x(4); 
    v2 = x(2); w2 = x(5); z2 = x(6); n2 = x(7); p2 = x(8); m2 = x(9); h2 = x(10);

    % membrane currents
    IL1     = P.gL1*(P.EL1-v1);
    IL2     = P.gL2*(P.EL2-v2);
    IKLT1   = P.gKLT1*w1^4*z1*(P.EK-v1);
    IKLT2   = P.gKLT2*w2^4*z2*(P.EK-v2);
    IKHT2   = P.gKHT2*(0.85*n2^2 + 0.15*p2)*(P.EK-v2);
    INa2    = P.gNa2*m2^3*h2*(P.ENa-v2);
    
    % synaptic currents
    IsynE1  = P.gE*P.epsg(t-P.tE0)*(P.EsynE-v1);
    IsynE2   = 0;
    IsynI1  = P.gI1*P.ipsg(t-P.tI1)*(P.EsynI-v1);
    IsynI2  = P.gI2*P.ipsg(t-P.tI2)*(P.EsynI-v2);

    % DC input currents
    I1 = ((t>5) && (t<15))*P.I1;
    I2 = ((t>5) && (t<15))*P.I2;

    % dynamics
    dv1 = (IL1 + IKLT1 + IsynE1 + IsynI1 + P.gC*(v2-v1) + I1 ) / P.C1;
    dv2 = (IL2 + IKLT2 + IKHT2 + INa2 + IsynE2 + IsynI2 + P.gC*(v1-v2) + I2) / P.C2;
    dw1 = P.phi * ( P.winf(v1) - w1)/ P.tauw(v1); 
    dz1 = P.phi * ( P.zinf(v1) - z1)/ P.tauz(v1); 
    dw2 = P.phi * ( P.winf(v2) - w2)/ P.tauw(v2); 
    dz2 = P.phi * ( P.zinf(v2) - z2)/ P.tauz(v2); 
    dn2 = P.phi * ( P.ninf(v2) - n2)/ P.taun(v2); 
    dp2 = P.phi * ( P.pinf(v2) - p2)/ P.taup(v2); 
    dm2 = P.phi * ( P.minf(v2) - m2)/ P.taum(v2); 
    dh2 = P.phi * ( P.hinf(v2) - h2)/ P.tauh(v2); 
    
    % output for ODE solver
    dx = [dv1 ; dv2; dw1 ; dz1 ; dw2 ;dz2 ; dn2 ; dp2 ; dm2 ; dh2];


end

