
%Visual input tuneing curve demo BioNB330
%Bruce Land
%1March2004

%Three current inputs from receptors to a neuron
%each current is sinusiodal with angle of a visually presented line.

clear all

%build the GUI
%define the figure

colordef white

fig=figure(1);
clf
figTop = 650;
figRight = 800;
set(fig,'resize','off',...
    'doublebuffer','on',...
    'color',[1 1 1]*0.8,...
    'position',[100 100 figRight figTop],...
    'name', 'BioNB 330 Tuning Curve Simulator',...
    'numbertitle','off')

fig2=figure(2);
clf
set(fig2,'resize','off',...
    'doublebuffer','on',...
    'color',[1 1 1]*0.8,...
    'position',[900 200 200 150],...
    'name', 'BioNB 330 ISI histogram',...
    'visible', 'off', ...
    'CloseRequestFcn','set(fig2,''visible'',''off'')', ...
    'numbertitle','off')

set(0,'currentfigure',fig)

%===Main state controls=================
%define the Run button
x=10; y=figTop-25
plotbutton=uicontrol('style','pushbutton',...
    'string','ReRun', ...
    'fontsize',12, ...
    'position',[x,y,50,20], ...
    'callback', 'run=1;resetdone=0;measure=0;');

%define the stop button
erasebutton=uicontrol('style','pushbutton',...
    'string','Stop', ...
    'fontsize',12, ...
    'position',[x+65,y,50,20], ...
    'callback','freeze=1;');

%define the Quit button
quitbutton=uicontrol('style','pushbutton',...
    'string','Quit', ...
    'fontsize',12, ...
    'position',[figRight-65,y,50,20], ...
    'callback','stop=1;set(fig2,''CloseRequestFcn'','''');close all;');

%===define the regions on the window=======
%for in1, in2, synapse, V, output
wholefigure=axes('position', [0 0 1 1]);
set(wholefigure, 'color', [1 1 1]*0.8,'xlim',[0,1], 'ylim',[0,1])
axis off

%in3
rectangle('position',[.01,0.30,0.33,0.15], ...
    'FaceColor',[1 1 1]*0.5)
%in2
rectangle('position',[.01,0.46,0.33,0.15], ...
    'FaceColor',[1 1 1]*0.5)
%in1
rectangle('position',[.01,0.62,0.33,0.15], ...
    'FaceColor',[1 1 1]*0.5)
%volage
rectangle('position',[.35,0.30,0.33,0.65], ...
    'FaceColor',[1 1 1]*0.5)
%ooutput
rectangle('position',[.690,0.30,0.3,0.65], ...
    'FaceColor',[1 1 1]*0.5)

%define the voltage variable animation axes
Vanimation = axes ;
w=.4;
set(Vanimation,'xlim',[-1,1],...
    'ylim',[-1,1],...
    'position', [0.44,0.3,w,w])
axis off
axis square

%define the input line animation axes
Lanimation = axes ;
w=.15;
set(Lanimation,'xlim',[-1,1],...
    'ylim',[-1,1],...
    'position', [0.1,0.8,w,w])
axis off
axis square

%define the plotters -- total height .30
plotterIN = axes('position',[0.05,0.05,.9,.08]);
set(plotterIN,'color','k') ;
set(plotterIN,'xlim',[0,180]);
plotterOUT = axes('position',[0.05,0.135,.9,.16]);
set(plotterOUT,'color','k') ;
set(plotterOUT,'xlim',[0,180]);


%===Input 1 controls========================
%identifier
x=10;y=480;
psat = 0.6 %saturation of plots
uicontrol('style','text',...
    'string','Input 1', ...
    'fontsize',12, ...
    'backgroundcolor',[1,psat,psat],...
    'position',[x,y,50,20])

%sinewave input
y=y-25;
%offset contrl
uicontrol('style','text',...
    'string','Sine', ...
    'fontsize',12, ...
    'position',[x,y,40,20]);
xt = x +45;
uicontrol('style','text',...
    'string','offset', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt = xt + 45;
sineoffsetUI = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
%frequency control
y=y-25;
xt = x + 45;
uicontrol('style','text',...
    'string','freq', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+45;
sinefreqUI = uicontrol('style','edit', ...
    'string','1.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
%weight
xt=xt+60
uicontrol('style','text',...
    'string','weight', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+45;
weightUI = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
%phase control
y=y-25;
xt = x +45;
uicontrol('style','text',...
    'string','phase', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+45;
sinephaseUI = uicontrol('style','edit', ...
    'string','90', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
% 
%===Input 2 controls========================
%%identifier
x=10;y=375;
uicontrol('style','text',...
    'string','Input 2', ...
    'fontsize',12, ...
    'backgroundcolor',[psat 1 psat],...
    'position',[x,y,50,20])

%sinewave input
y=y-25;
%offset contrl
uicontrol('style','text',...
    'string','Sine', ...
    'fontsize',12, ...
    'position',[x,y,40,20]);
xt = x +45;
uicontrol('style','text',...
    'string','offset', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt = xt + 45;
sineoffsetUI2 = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
%frequency control
y=y-25;
xt = x + 45;
uicontrol('style','text',...
    'string','freq', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+45;
sinefreqUI2 = uicontrol('style','edit', ...
    'string','1.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
%weight
xt=xt+60
uicontrol('style','text',...
    'string','weight', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+45;
weightUI2 = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
%phase control
y=y-25;
xt = x +45;
uicontrol('style','text',...
    'string','phase', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+45;
sinephaseUI2 = uicontrol('style','edit', ...
    'string','0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
% 

%===Input 3 controls========================
%identifier
x=10;y=272;
uicontrol('style','text',...
    'string','Input 3', ...
    'fontsize',12, ...
    'backgroundcolor',[psat psat 1],...
    'position',[x,y,50,20])

%sinewave input
y=y-25;
%offset contrl
uicontrol('style','text',...
    'string','Sine', ...
    'fontsize',12, ...
    'position',[x,y,40,20]);
xt = x +45;
uicontrol('style','text',...
    'string','offset', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt = xt + 45;
sineoffsetUI3 = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
%frequency control
y=y-25;
xt = x + 45;
uicontrol('style','text',...
    'string','freq', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+45;
sinefreqUI3 = uicontrol('style','edit', ...
    'string','1.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
%weight
xt=xt+60
uicontrol('style','text',...
    'string','weight', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+45;
weightUI3 = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
%phase control
y=y-25;
xt = x +45;
uicontrol('style','text',...
    'string','phase', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+45;
sinephaseUI3 = uicontrol('style','edit', ...
    'string','-90', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);

%===integration controls====================
x=290;y=596;
uicontrol('style','text',...
    'string','voltage function', ...
    'fontsize',12, ...
    'backgroundcolor',[1,1,1],...
    'position',[x,y,100,20])
x=300;y=550
uicontrol('style','text',...
    'string','Rin', ...
    'fontsize',12, ...
    'position',[x,y,50,20]);
%define Rin
inputRUI = uicontrol('style','edit', ...
    'string','2.0', ...
    'fontsize',12, ...
    'position',[x+60,y,50,20]);
y=y-25;
uicontrol('style','text',...
    'string','Time Const', ...
    'fontsize',12, ...
    'position',[x,y,100,20]);
%define dk rate
DKrateUI = uicontrol('style','edit', ...
    'string','20.0', ...
    'fontsize',12, ...
    'position',[x+110,y,50,20]);

%show C and R currents
% y=y-100;
% uicontrol('style','pushbutton', ...
%     'string','Show Currents', ...
%     'fontsize',12, ...
%     'callback','showcurrent=1',...
%     'position',[x+30,y,100,20]);
% y=y-30;
% uicontrol('style','pushbutton', ...
%     'string','Hide Currents', ...
%     'fontsize',12, ...
%     'callback','showcurrent=0',...
%     'position',[x+30,y,100,20]);

%measure controls
y=y-200;
uicontrol('style','pushbutton', ...
    'string','Measure', ...
    'fontsize',12, ...
    'callback','measure=1;freeze=1;',...
    'position',[x+30,y,100,20]);
y=y-30;
uicontrol('style','pushbutton', ...
    'string','stop measure', ...
    'fontsize',12, ...
    'callback','measure=0;',...
    'position',[x+30,y,100,20]);
y=y-30;
tmeasure = uicontrol('style','text', ...
    'string','---', ...
    'fontsize',12, ...
    'position',[x+20,y,80,20]);
uicontrol('style','text', ...
    'string','Delta T', ...
    'fontsize',12, ...
    'position',[x+80,y,70,20]);
y=y-30;
vmeasure = uicontrol('style','text', ...
    'string','---', ...
    'fontsize',12, ...
    'position',[x+20,y,80,20]);
uicontrol('style','text', ...
    'string','Delta V', ...
    'fontsize',12, ...
    'position',[x+80,y,70,20]);

%threshold 
x=600
y=350;
uicontrol('style','text',...
    'string','Threshold', ...
    'fontsize',12, ...
    'position',[x,y,75,20]);
y=y-25;
uicontrol('style','text',...
    'string','level', ...
    'fontsize',12, ...
    'position',[x,y,50,20]);
thresholdlevelUI = uicontrol('style','edit', ...
    'string','0.8', ...
    'fontsize',12, ...
    'position',[x+60,y,50,20]);
%threshold noise
y=y-25;
uicontrol('style','text',...
    'string','noise level', ...
    'fontsize',12, ...
    'position',[x,y,100,20]);
thresholdnoiselevelUI = uicontrol('style','edit', ...
    'string','0', ...
    'fontsize',12, ...
    'position',[x+110,y,50,20]);
%threshold adaption
y=y-25;
uicontrol('style','text',...
    'string','adapt level', ...
    'fontsize',12, ...
    'position',[x,y,100,20]);
thresholdadaptlevelUI = uicontrol('style','edit', ...
    'string','0', ...
    'fontsize',12, ...
    'position',[x+110,y,50,20]);
y=y-25;
uicontrol('style','text',...
    'string','adapt decay', ...
    'fontsize',12, ...
    'position',[x,y,100,20]);
thresholdadaptDKUI = uicontrol('style','edit', ...
    'string','1.0', ...
    'fontsize',12, ...
    'position',[x+110,y,50,20]);


%===output controls========================
%identifier
x=560;y=596;
uicontrol('style','text',...
    'string','Output Function', ...
    'fontsize',12, ...
    'backgroundcolor',[1,1,1],...
    'position',[x,y,100,20])

%ISI plotter
x=710; y=205
uicontrol('style','pushbutton',...
    'string','Reset ISI', ...
    'fontsize',12, ...
    'position',[x,y,70,20], ...
    'callback', 'resetisi=1;');

isienable = uicontrol('style','pushbutton',...
    'string','Show ISI', ...
    'fontsize',12, ...
    'callback','drawisi=1;', ...
    'position',[x-120,y,80,20]);

%radio buttons for spike/sigmoid
x=570; y=560
spikemodeui = uicontrol('style','radiobutton',...
    'string','Spike', ...
    'fontsize',12, ...
    'position',[x,y,70,20], ...
    'value',1,...
    'callback', [...
        'set(sigmoidmodeui,''value'',0);',...
        'set(linearmodeui,''value'',0);']);
set(gcf,'currentaxes',wholefigure)
spikeLight = ...
    rectangle('position',[.83,.86,.02,.03],'linewidth',2,'facecolor','k')
y=y-40;
sigmoidmodeui = uicontrol('style','radiobutton',...
    'string','Sigmoid', ...
    'fontsize',12, ...
    'position',[x,y,70,20], ...
    'callback', [...
        'set(spikemodeui,''value'',0);',...
        'set(linearmodeui,''value'',0);'...
    ]);
uicontrol('style','text',...
    'string','steepness', ...
    'fontsize',12, ...
    'position',[x+80,y,80,20]);
sigmoidNui = uicontrol('style','edit', ...
    'string','1.0', ...
    'fontsize',12, ...
    'position',[x+160,y,50,20]);
sigmoidplot = axes('position',[.8, .6, .1, .1])
set(gca,'xlim',[0 1.5],'ylim',[0 1])
tempx = [0:.05:2];
sigline = plot(tempx, tempx./(tempx+0.8),'k');
hold on
sigdot = plot(0.2,0.2,'ro');
set(sigmoidplot,'visible','off')
set(sigline,'visible','off')
set(sigdot,'visible','off')
y=y-40;
linearmodeui= uicontrol('style','radiobutton',...
    'string','Linear', ...
    'fontsize',12, ...
    'position',[x,y,70,20], ...
    'callback', [...
        'set(spikemodeui,''value'',0);'...
        'set(sigmoidmodeui,''value'',0);',...
    ]);

%===initial state of buttons
stop= 0; %wait for a quit button push
run = 0; %wait for a draw 
resetdone=0; %are we restarting?

%===Now the main simulation loop============
while (stop==0)
    tic %set timer for loop timing
    
    if (resetdone==0) %set up the state variables
        resetdone = 1;
        freeze = 0; %wait for a freeze
        count=0 ; %spike counter
        resetisi=0 ; %isi starts cleared
        spikeMode=1; %start in spikeing mode
        drawisi=0;
        showcurrent=0 ; %dont show R and C currents
        measure=0 ; %enable interaction with plot
        initialpt(1,1)=0; %point for initial measurement
        initialpt(1,2)=0;
        %===state variables========================
        t=0; dt=1;
        vrest = 0.0;
        threshold = .8;
        v = vrest ;
        spike=0;
        spiketime = -inf;
        lastspiketime = -inf;
        isi=[];
        angle = 0;
        dangle = 1;
        stepcounter = 1;
        %===Animations=============================
        %now the two voltage animation objects
        set(fig,'currentaxes',Vanimation)
        cla
        Trect = rectangle('position',[-.1,-0.7,0.2,threshold+.2],'linewidth',2)
        Vrect = rectangle('position',[-.1,-0.7,0.2,v+.2],'linewidth',2,'facecolor','red')
        rectangle('position',[-.1,-0.5,0.2,.01],'linewidth',2)
        %and the input animation object
        set(fig,'currentaxes',Lanimation)
        cla
        line([-1 1], [0 0], 'color','k') %axis
        line([0 0], [-1 1],'color','k') %axis
        stimline = line([0 1], [0 0],'color','r') %stimulus
        
        %===output plot variables==================
        Vout=vrest*ones(1,180);
        Sout=zeros(1,180);
        I1out=zeros(1,180);
        I2out=zeros(1,180);
        I3out=zeros(1,180);
        Itotal=zeros(1,180);
        %ICout=zeros(1,180);
        Tout=threshold*ones(1,180);
        
        % and init the output  plot
        set(fig,'currentaxes',plotterOUT)
        cla
        VP = plot(Vout,'r','linewidth',2);
        hold on
        SP = plot(Sout+1,'w');
        TP = plot(Tout,'r:');
        IP = plot(Itotal,'color',[1 1 1],'linewidth',2); 
        set(plotterOUT,'ylim',[-1,5])
        set(plotterOUT,'color',0.2*[1 1 1]) ;
        set(plotterOUT,'xtick',[]);
        
        % and init the input  plot
        set(fig,'currentaxes',plotterIN)
        cla
        IP1 = plot(I1out,'color',[1,psat,psat]);
        hold on
        IP2 = plot(I2out,'color',[psat 1 psat]);
        IP3 = plot(I3out,'color',[psat psat 1]);
        set(plotterIN,'ylim',[-1.2,1.2])
        set(plotterIN,'color',0.2*[1 1 1]) ;
        set(plotterIN,'xtick',[0:10:180])
    end
    
    if (run==1) %running?
        t=t+dt;
        angle = angle + dangle;
        wangle = angle *(2*pi)/360;
        if angle==180
            run=0;
        end
        %===draw the voltage animation=========================      
        set(Trect,'position',[-.1,-0.7,0.2,max(.001,threshold+.2)])
        set(Vrect,'position',[-.1,-0.7,0.2,max(.001,v+.2)])
        
        %===draw the input bar animation
        set(stimline,...
            'xdata',[0 cos(wangle)],...
            'ydata',[0 sin(wangle)]) ;
        
        %===update the state variables========================
        Rin = str2num(get(inputRUI,'string'));
        DKrate = str2num(get(DKrateUI,'string'));
        
        sineoffset = str2num(get(sineoffsetUI,'string'));
        sinefreq = str2num(get(sinefreqUI,'string'));
        sinephase = str2num(get(sinephaseUI,'string'))/57.29;
        weight = str2num(get(weightUI,'string'));
        
        sineoffset2 = str2num(get(sineoffsetUI2,'string'));
        sinefreq2 = str2num(get(sinefreqUI2,'string'));
        sinephase2 = str2num(get(sinephaseUI2,'string'))/57.29;
        weight2 = str2num(get(weightUI2,'string'));
        
        sineoffset3 = str2num(get(sineoffsetUI3,'string'));
        sinefreq3 = str2num(get(sinefreqUI3,'string'));
        sinephase3 = str2num(get(sinephaseUI3,'string'))/57.29;
        weight3 = str2num(get(weightUI3,'string'));
        
        spikemode = get(spikemodeui,'value');
        linearmode = get(linearmodeui,'value');
        sigmoidmode = get(sigmoidmodeui,'value');
        sigmoidN = str2num(get(sigmoidNui,'string'));
        
        %total current
        I1 = (sineoffset + sin(wangle*sinefreq+sinephase));
        I2 = (sineoffset2 + sin(wangle*sinefreq2+sinephase2));
        I3 = (sineoffset3 + sin(wangle*sinefreq3+sinephase3));
        stimTotal = weight*I1 + weight2*I2 + weight3*I3 ;
        
        %update voltage -- DKrate is actually time constant
        v = v + dt/DKrate * (Rin * stimTotal - (v-vrest));
        
        %make the spike or sigmoid
        if spikemode %do the spike thing
            set(spikeLight,'visible','on')
            set(sigmoidplot,'visible','off')
            set(sigline,'visible','off')
            set(sigdot,'visible','off')
            if (v>=threshold)
                spike=1;
                v=vrest;
                lastspiketime = spiketime;
                spiketime = t;
                count=count+1;
                isi(count)=spiketime-lastspiketime;
                newspike=1;
                set(spikeLight,'facecolor','r')
            else 
                spike=0;
                %set(gcf,'currentaxes',wholefigure)
                set(spikeLight,'facecolor','k')
            end
        elseif sigmoidmode %do the sigmoid thang
            set(spikeLight,'visible','off') 
            set(sigmoidplot,'visible','on')
            set(sigline,'visible','on')
            set(sigdot,'visible','on')
            spike = v^sigmoidN/(v^sigmoidN + threshold^sigmoidN);
            set(sigdot,'xdata',v,'ydata',spike)
            set(sigline,'ydata',tempx.^sigmoidN./(tempx.^sigmoidN+thresholdlevel.^sigmoidN))
        elseif linearmode
            set(spikeLight,'visible','off')
            set(sigmoidplot,'visible','off')
            set(sigline,'visible','off')
            set(sigdot,'visible','off')
            spike = (v-threshold)*(v>threshold);
        end
        
        %update the threshold
        thresholdlevel = str2num(get(thresholdlevelUI,'string'));
        thresholdnoiselevel = str2num(get(thresholdnoiselevelUI,'string'));
        thresholdadaptlevel=str2num(get(thresholdadaptlevelUI,'string'));
        thresholdadaptDK=str2num(get(thresholdadaptDKUI,'string'));
        threshold = threshold - ...
            thresholdadaptDK*(threshold-thresholdlevel) + ...
            thresholdadaptlevel * spike +thresholdnoiselevel*randn ;
        
        %bound the voltage below
%         if (v<=-0.2)
%             v=-0.199;
%         end
%         
%         %bound the threshold below
%         if (threshold<=0.001)
%             threshold=.001;
%         end
    
        %===draw the history plot=======================
        %update output buffers, then draw the data
        %set(fig,'currentaxes',plotter)
        
        Vout(stepcounter) = v;
        Sout(stepcounter) = spike;
        I1out(stepcounter) = I1;
        I2out(stepcounter) = I2;
        I3out(stepcounter) = I3;
        Itotal(stepcounter) = stimTotal;
        %ICout(stepcounter) = stimTotal-IRout(stepcounter);  
        Tout(stepcounter) = threshold;
        stepcounter = stepcounter + 1;
        
        set(VP,'ydata',Vout);
        set(SP,'ydata',Sout+1);
        set(IP1,'ydata',I1out);
        set(IP2,'ydata',I2out);
        set(IP3,'ydata',I3out);
        set(IP,'ydata',Itotal);
        %set(ICP,'ydata',ICout);
        set(TP,'ydata',Tout);
        
    end %running
    
    %===reset the isi history=========================
    if resetisi==1
        isi=[];
        resetisi=0;
        count=0;
    end
    
    %===draw the ISI histogram========================
    if drawisi==1
        drawisi=0;
        set(0,'currentfigure',fig2)
        set(fig2,'visible','on')
        set(fig2,'position',[900,200,200,150])
        nbins = max(10, length(isi)/10);
        hist(isi(2:end),nbins)
        drawnow
        set(0,'currentfigure',fig)
        drawnow
    end
        
    %===stop the plot=================================
    if (freeze==1)
        run = 0;
        freeze = 0;
    end
    
    %===measure stuff from the plot===================
    if measure==1
        set(gcf,'WindowButtonDownFcn',...
            [  ...
                'plot(initialpt(1,1),initialpt(1,2),''d'',''color'',.2*[1 1 1]);' ...
                'initialpt=get(gca,''currentpoint'');'...
                'plot(initialpt(1,1),initialpt(1,2),''yd'');' ...
            ]);
        
        set(gcf,'WindowButtonMotionFcn',...
            [...
                'finalpt=get(gca,''currentpoint'');'...
                'dtm=(finalpt(1,1)-initialpt(1,1));'...
                'dv=(finalpt(1,2)-initialpt(1,2));'...
                'set(tmeasure,''string'',num2str(dtm,''%6.2f''));'...
                'set(vmeasure,''string'',num2str(dv,''%6.2f''));'...
            ]);
    end
    
    if measure==0
        set(gcf,'WindowButtonDownFcn','')
        set(gcf,'WindowButtonMotionFcn','')
        set(tmeasure,'string','---')
        set(vmeasure,'string','---')
    end
    
    %===draw and control display speed================
    drawnow  %need this in the loop for controls to work
    if (toc<0.05) 
        pause(0.05-toc) %20 frames/sec
    end
end %main simulation loop




