function neuron_current_injection_sim
% neuron_current_injection_sim
% 2 s neuron with 200 ms current injection (t = 1.0–1.2 s).
% UI: I_inj (pA), K Conductance, Na Conductance, Replay.
% Plots:
%  • Top-left: Vm (black), fixed [-90, +55] mV
%  • Bottom-left: I_inj (black)
%  • Top-right: AP overlay centered on NEW peak, window [-0.5 ms, +2.0 ms]; previous=gray, new=red.
% Auto layout: pixel-based, responsive to figure resize. Fonts scale gently.

    % ===== Time / model =====
    EL_mV      = -65;   Vth_mV   = -48;   Vpeak_mV = 40;
    tref       = 0.002; Rm       = 100e6; taum      = 0.020;
    T          = 2.0;   dt       = 1e-4;  t         = 0:dt:T;   % 0.1 ms

    % SI units
    gL = 1/Rm; C = taum/Rm;
    EL = EL_mV*1e-3; Vth = Vth_mV*1e-3; Vpeak = Vpeak_mV*1e-3;

    % ===== Ionic parameters =====
    EK       = -80e-3;
    tauK     = 0.100;                 % s  (AHP decay)
    gK_step0 = 30e-9;                 % S  per-spike AHP step at K-gain=1
    gK_fast  = 200e-9;                % S  phasic K+ when depolarized (scaled by K-gain)
    ENa      =  60e-3;
    gNa_fast = 400e-9;                % S  fast Na near threshold (scaled by Na-gain)
    thrNaWin = 5e-3;                  % V  window below threshold where Na assist engages

    % Threshold adaptation (for single spike at rheobase)
    th_step0_mV = 5; tauTh = 0.100; th_step0 = th_step0_mV * 1e-3;

    % Injection timing
    inj.t0 = 1.0; inj.t1 = 1.2;

    % Precompute decay factors
    alphaK  = exp(-dt/tauK);
    alphaTh = exp(-dt/tauTh);

    % ===== Figure & layout =====
    % Preferred size from your measured figure: [left bottom width height]
    defaultPos = [3969 776 775 558];

    f = figure('Name','Neuron Current Injection','NumberTitle','off',...
        'Color','w','Units','pixels','Position',defaultPos);

    % Reposition to secondary monitor (and cap to screen) if present
    monPos = get(0,'MonitorPositions');
    if size(monPos,1) >= 2
        mon = monPos(2,:);
        figW = min(775, mon(3)-200);
        figH = min(558, mon(4)-200);
        x = mon(1) + 80;
        y = mon(2) + mon(4) - figH - 80;
        set(f,'Position',[x y figW figH]);
    else
        % Ensure the default fits the primary monitor too
        mon = monPos(1,:);
        p = get(f,'Position');
        p(3) = min(p(3), mon(3)-120);
        p(4) = min(p(4), mon(4)-120);
        set(f,'Position',p);
    end

    % 2x2 layout (plots in three tiles, controls in the 4th)
    tl = tiledlayout(f,2,2,'TileSpacing','compact','Padding','compact');

    axVm = nexttile(tl,1); hold(axVm,'on'); box(axVm,'on');
    axI  = nexttile(tl,3); hold(axI,'on'); box(axI,'on');
    axAP = nexttile(tl,2); hold(axAP,'on'); box(axAP,'on');

    % Create a placeholder tile, then replace with a panel
    axCtrl = nexttile(tl,4);
    ctrlPosNorm = get(axCtrl,'Position'); delete(axCtrl);
    ctrlPanel = uipanel('Parent',f,'Title','Controls','Units','normalized',...
                        'Position',ctrlPosNorm,'BackgroundColor','w','FontWeight','bold');

    % Axes formatting
    xlabel(axVm,'Time (s)'); ylabel(axVm,'V_m (mV)'); title(axVm,'Membrane potential');
    ylim(axVm,[-90 55]); xlim(axVm,[0 T]);
    yline(axVm, Vth_mV,'--','Threshold','LabelVerticalAlignment','bottom');

    xlabel(axI,'Time (s)'); ylabel(axI,'Current Injection (pA)'); title(axI,'Injected current');
    ylim(axI,[-5 305]); xlim(axI,[0 T]);

    linkaxes([axVm, axI], 'x');

    title(axAP,'Single action potential (−0.5 to +2.0 ms around peak)');
    xlabel(axAP,'Time (ms)'); ylabel(axAP,'V_m (mV)'); ylim(axAP,[-90 55]);

    % Lines (black)
    vmLine = plot(axVm, t(1), EL_mV, 'k', 'LineWidth', 1.5);
    iLine  = plot(axI,  t(1), 0,     'k', 'LineWidth', 1.5);

    % ===== Controls (created once; positions set by relayout) =====
    hLblI   = uicontrol('Style','text','Parent',ctrlPanel,'String','Current Injection (pA)',...
        'BackgroundColor','w','ForegroundColor','k','HorizontalAlignment','left');
    editBox = uicontrol('Style','edit','Parent',ctrlPanel,'String','0','Callback',@onEditChanged);

    kTitle  = uicontrol('Style','text','Parent',ctrlPanel,'String','K Conductance',...
        'BackgroundColor','w','ForegroundColor','k','HorizontalAlignment','center');

    kSlider = uicontrol('Style','slider','Parent',ctrlPanel,'Min',0,'Max',1,'Value',0.5);
    k0Lbl   = uicontrol('Style','text','Parent',ctrlPanel,'String','0','BackgroundColor','w','HorizontalAlignment','left');
    knLbl   = uicontrol('Style','text','Parent',ctrlPanel,'String','normal','BackgroundColor','w','HorizontalAlignment','center');
    kMaxLbl = uicontrol('Style','text','Parent',ctrlPanel,'String','too much','BackgroundColor','w','HorizontalAlignment','right');

    naTitle  = uicontrol('Style','text','Parent',ctrlPanel,'String','Na Conductance',...
        'BackgroundColor','w','ForegroundColor','k','HorizontalAlignment','center');

    naSlider = uicontrol('Style','slider','Parent',ctrlPanel,'Min',0,'Max',1,'Value',0.5);
    na0Lbl   = uicontrol('Style','text','Parent',ctrlPanel,'String','0','BackgroundColor','w','HorizontalAlignment','left');
    nanLbl   = uicontrol('Style','text','Parent',ctrlPanel,'String','normal','BackgroundColor','w','HorizontalAlignment','center');
    naMaxLbl = uicontrol('Style','text','Parent',ctrlPanel,'String','too much','BackgroundColor','w','HorizontalAlignment','right');

    replayBtn = uicontrol('Style','pushbutton','Parent',ctrlPanel,'String','Replay','Callback',@(~,~) replay());

    % ===== Shared side labels spanning both sliders =====
sideLeftLbl  = uicontrol('Style','text','Parent',ctrlPanel,'String','0', ...
    'BackgroundColor','w','HorizontalAlignment','left');
sideRightLbl = uicontrol('Style','text','Parent',ctrlPanel,'String','too much', ...
    'BackgroundColor','w','HorizontalAlignment','right');

% Hide the old per-slider left/right labels
set(k0Lbl,'Visible','off');   set(kMaxLbl,'Visible','off');
set(na0Lbl,'Visible','off');  set(naMaxLbl,'Visible','off');


    % Hook resize → auto layout
    f.SizeChangedFcn = @relayout;
    relayout();  % do once now

    % ===== Timer + run storage for AP overlay =====
    animTimer = [];
    prevVm_mV = []; prevSpikes = [];
    lastVm_mV = []; lastSpikes = [];
    f.CloseRequestFcn = @onClose;

    % Initial run
    doComputeAndPlay(str2double(editBox.String), kSlider.Value, naSlider.Value);

    % ===== Callbacks =====
    kSlider.Callback  = @(~,~) onAnyChanged();
    naSlider.Callback = @(~,~) onAnyChanged();
    function onAnyChanged()
        IpA = str2double(editBox.String); if isnan(IpA), IpA = 0; end
        doComputeAndPlay(IpA, kSlider.Value, naSlider.Value);
    end

    function onClose(~,~), stopAnim(); if isvalid(f), delete(f); end, end

    function onEditChanged(src,~)
        IpA = str2double(src.String); if isnan(IpA), IpA = 0; end
        IpA = max(0, min(300, IpA)); src.String = num2str(IpA);
        doComputeAndPlay(IpA, kSlider.Value, naSlider.Value);
    end

    function replay
        IpA = str2double(editBox.String); if isnan(IpA), IpA = 0; end
        doComputeAndPlay(IpA, kSlider.Value, naSlider.Value);
    end

    % ======== Auto layout (pixel-based, resize-aware) ========
function relayout(~,~)
    % Panel-local size
    p = getpixelposition(ctrlPanel);  % [x y w h]
    w = p(3); h = p(4);

    % Base metrics
    padB  = 14; gapB = 10; rowB = 30; labHB = 18; btnHB = 34; editWB = 110;

    % Compute height budget and scale
    Hreq = padB + ...                    % top pad
           rowB + gapB + ...             % Current injection
           rowB + gapB + labHB + ...     % K (title+slider) + mid label
           rowB + gapB + labHB + ...     % Na (title+slider) + mid label
           btnHB + padB;                 % button + bottom pad
    s = min(1, max(0.6, h / Hreq));

    pad  = round(padB  * s);
    gap  = round(gapB  * s);
    row  = max(20, round(rowB  * s));
    labH = max(14, round(labHB * s));
    btnH = max(24, round(btnHB * s));
    editW = editWB;

    % Fonts
    fsTitle = max(10, min(14, round(11*s + 3*(w/700))));
    fsText  = max(9,  min(12, round(10*s + 2*(w/800))));
    fsBtn   = max(10, min(14, round(11*s + 2*(w/900))));

    % Panel-local columns
    x0   = pad;
    y    = h - pad - row;
    colW = w - 2*pad;

    % ---- Current injection (label + edit) ----
    set(hLblI,  'Units','pixels','Position',[x0 y colW-editW-8 row], ...
        'FontSize',fsTitle,'FontWeight','bold');
    set(editBox,'Units','pixels','Position',[x0+colW-editW y editW row], ...
        'FontSize',fsTitle,'FontWeight','bold','HorizontalAlignment','center');
    y = y - (row + gap);

    % ---- Shared side labels geometry ----
    sideLblW = max(36, round(40 * s));           % width for the left "0"
    sideRbw  = max(78, round(84 * s));           % width for the right "too much"
    sidePad  = 6;                                % spacing between label and slider

    % Slider track geometry (same for both K and Na)
    sliderX = x0 + sideLblW + sidePad;           % left edge of slider track
    sliderW = colW - sideLblW - sideRbw - 2*sidePad; % narrowed slider width

    % ---- K block ----
    set(kTitle,'Units','pixels','Position',[x0 y colW row], ...
        'FontSize',fsTitle,'FontWeight','bold');
    y = y - row;

    set(kSlider,'Units','pixels','Position',[sliderX y sliderW row]);
    % center "normal" directly under K slider
    y = y - labH;
    set(knLbl,'Units','pixels','Position',[sliderX + (sliderW-60)/2 y 60 labH], ...
        'FontSize',fsText,'FontWeight','bold','HorizontalAlignment','center','Visible','on');
    y = y - gap;

    % ---- Na block ----
    set(naTitle,'Units','pixels','Position',[x0 y colW row], ...
        'FontSize',fsTitle,'FontWeight','bold');
    y = y - row;

    set(naSlider,'Units','pixels','Position',[sliderX y sliderW row]);
    % center "normal" directly under Na slider
    y = y - labH;
    set(nanLbl,'Units','pixels','Position',[sliderX + (sliderW-60)/2 y 60 labH], ...
        'FontSize',fsText,'FontWeight','bold','HorizontalAlignment','center','Visible','on');

    % ---- Place the shared side labels spanning both sliders ----
    % vertical band covering both slider rows (midpoint between them)
    % We’ll align them vertically with the midpoint between the two sliders.
    % Compute Y of K slider and Na slider to average:
    kPos  = get(kSlider,'Position');   % [x y w h] in panel pixels
    naPos = get(naSlider,'Position');
    midY  = (kPos(2) + naPos(2) + kPos(4)/2 + naPos(4)/2)/2 - labH/2;

    set(sideLeftLbl, 'Units','pixels','Position',[x0, midY, sideLblW, labH], ...
        'FontSize',fsText,'FontWeight','bold','HorizontalAlignment','left','Visible','on');
    set(sideRightLbl,'Units','pixels','Position',[x0 + colW - sideRbw, midY, sideRbw, labH], ...
        'FontSize',fsText,'FontWeight','bold','HorizontalAlignment','right','Visible','on');

    % ---- Replay button ----
    set(replayBtn,'Units','pixels','Position',[w - pad - 120, pad, 120, btnH], ...
        'FontSize',fsBtn,'FontWeight','bold');
end





    % ======== Sim + animate ========
    function doComputeAndPlay(IpA, kslider, naslider)
        stopAnim();

        % AP panel "Updating..." while animating
        cla(axAP); box(axAP,'on'); hold(axAP,'on');
        text(axAP,0.5,0.5,'Updating…','Units','normalized',...
            'HorizontalAlignment','center','FontWeight','bold');
        xlabel(axAP,'Time (ms)'); ylabel(axAP,'V_m (mV)'); ylim(axAP,[-90 55]);

        % Map sliders: 0.5=1×, clamp to 0..2×
        k_gain  = max(0, min(2, kslider  / 0.5));
        na_gain = max(0, min(2, naslider / 0.5));

        % Injection waveform (pA) & simulate
        Iinj_pA = zeros(size(t));
        Iinj_pA(t >= inj.t0 & t < inj.t1) = IpA;
        [Vm_mV, spikeTimes] = simulate_with_ions(Iinj_pA*1e-12, k_gain, na_gain, naslider);

        % Store for overlay
        lastVm_mV = Vm_mV; lastSpikes = spikeTimes;

        % Prepare animation
        set(vmLine,'XData',t(1),'YData',Vm_mV(1));
        set(iLine, 'XData',t(1),'YData',Iinj_pA(1));
        title(axVm, sprintf('Vm (I=%d pA, K=%.2fx, Na=%.2fx) — %d spike(s)',...
                            IpA, k_gain, na_gain, numel(spikeTimes)));

        animTimer = timer('ExecutionMode','fixedSpacing','Period',0.01,...
                          'TimerFcn',@tick,'StartDelay',0.02);
        idx = 1; N = numel(t); batch = max(1, floor(0.01/dt));
        start(animTimer);

        function tick(~,~)
            if ~isvalid(f), stopAnim(); return; end
            i2 = min(N, idx+batch);
            set(vmLine,'XData',t(1:i2),'YData',Vm_mV(1:i2));
            set(iLine, 'XData',t(1:i2),'YData',Iinj_pA(1:i2));
            drawnow limitrate; idx = i2 + 1;

            if idx > N
                stopAnim();
                updateAPOverlay(prevVm_mV, prevSpikes, lastVm_mV, lastSpikes);
                prevVm_mV = lastVm_mV; prevSpikes = lastSpikes;
            end
        end
    end

    function stopAnim
        if ~isempty(animTimer) && isvalid(animTimer)
            try, stop(animTimer); catch, end
            try, delete(animTimer); catch, end
        end
        animTimer = [];
    end

    % ===== AP overlay: NEW centered, [-0.5 ms, +2.0 ms]; old gray, new red =====
    function updateAPOverlay(prevVm, prevSpk, newVm, newSpk)
        cla(axAP); box(axAP,'on'); hold(axAP,'on');

        if isempty(newSpk)
            text(axAP,0.5,0.5,'No spike','Units','normalized',...
                'HorizontalAlignment','center','FontWeight','bold');
            xlabel(axAP,'Time (ms)'); ylabel(axAP,'V_m (mV)'); ylim(axAP,[-90 55]);
            return;
        end

        % Find NEW peak and define window [-0.5 ms, +2.0 ms]
        tpeak_new = find_peak_time(newVm, newSpk);
        t0 = max(0,  tpeak_new - 0.0005);
        t1 = min(T,  tpeak_new + 0.0020);

        % High-res grid for smooth plot
        dt_zoom = 2e-5; t_zoom = t0:dt_zoom:t1;
        Vm_new  = interp1(t, newVm,  t_zoom, 'pchip');
        if ~isempty(prevVm)
            Vm_prev = interp1(t, prevVm, t_zoom, 'pchip');
        else
            Vm_prev = nan(size(t_zoom));
        end

        % Plot previous (gray) then new (red)
        if any(~isnan(Vm_prev))
            plot(axAP, t_zoom*1e3, Vm_prev, 'Color',[0.7 0.7 0.7], 'LineWidth',1.5);
        end
        plot(axAP, t_zoom*1e3, Vm_new, 'r', 'LineWidth',1.8);

        xlim(axAP,[t0 t1]*1e3); ylim(axAP,[-90 55]);
        title(axAP, sprintf('AP overlay: prev (gray) vs new (red). New peak @ %.3f ms', tpeak_new*1e3));
        xlabel(axAP,'Time (ms)'); ylabel(axAP,'V_m (mV)');
        if any(~isnan(Vm_prev)), legend(axAP,{'Previous','New'},'Location','best'); end
    end

    % Find first-spike peak time (in seconds)
    function tp = find_peak_time(Vm_mV, spkTimes)
        ts = spkTimes(1);
        winSearch = 0.005; % ±5 ms around threshold-crossing time
        idxSearch = (t >= ts - winSearch) & (t <= ts + winSearch);
        if ~any(idxSearch), idxSearch = (t >= ts) & (t <= min(T, ts + 0.010)); end
        [~, iRel] = max(Vm_mV(idxSearch));
        idxVec = find(idxSearch);
        kpeak  = idxVec(1) + iRel - 1;
        tp     = t(kpeak);
    end

    % ===== Model with K & Na modulation + threshold adaptation =====
    function [Vm_mV, spikeTimes] = simulate_with_ions(I_A, k_gain, na_gain, na_slider_raw)
        v = EL; gK = 0; refUntil = -inf; spikeTimes = [];
        Vth_eff = Vth;                      % adaptive threshold starts at baseline
        Vm_mV = zeros(size(t)); Vm_mV(1) = v*1e3;

        for k = 2:numel(t)
            % Decays
            gK      = gK * alphaK;                        % AHP K+
            Vth_eff = Vth + (Vth_eff - Vth) * alphaTh;    % threshold relaxes

            % Phasic K+ when depolarized (scaled by K gain)
            if v > (-20e-3), gK_extra = k_gain * gK_fast; else, gK_extra = 0; end

            % Fast Na near threshold (scaled by Na gain), only outside refractory
            if na_slider_raw > 0 && (t(k) >= refUntil) && (v > (Vth_eff - thrNaWin)) && (v < Vpeak)
                gNa = na_gain * gNa_fast;
            else
                gNa = 0;
            end

            % Integrate membrane (Euler)
            I_ion = -gL*(v - EL) ...
                    - (gK + gK_extra)*(v - EK) ...
                    + gNa*(ENa - v);
            dv = (I_ion + I_A(k-1)) / C * dt;
            v  = v + dv;

            % If Na slider is zero, we disable spikes entirely
            if na_slider_raw <= 0
                Vm_mV(k) = v*1e3; 
                continue;
            end

            % Spike detection (no forced reset; refractory gates re-spiking)
            if (t(k) >= refUntil) && (v >= Vth_eff)
                spikeTimes(end+1) = t(k); %#ok<AGROW>
                v = Vpeak;                 % brief peak
                Vm_mV(k) = v*1e3;

                % AHP step scaled by K gain
                gK = gK + k_gain * gK_step0;

                % Threshold adaptation bump (promotes single spike at rheobase)
                Vth_eff = Vth_eff + th_step0;

                % Enter refractory; keep integrating (K+ shapes downstroke)
                refUntil = t(k) + tref;
                continue;
            end

            Vm_mV(k) = v*1e3;
        end
    end
end