%Besic Computational Neuron Model for BioNB330
%See web page for user interface description
%Bruce Land Jan2004

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 Cell 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 plot button
x=10; y=figTop-25
plotbutton=uicontrol('style','pushbutton',...
    'string','Run', ...
    'fontsize',12, ...
    'position',[x,y,50,20], ...
    'callback', 'run=1;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

%in2
rectangle('position',[.01,0.30,0.33,0.32], ...
    'FaceColor',[1 1 1]*0.5)
%in1
rectangle('position',[.01,0.63,0.33,0.32], ...
    '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 plotter
plotter = axes('position',[0.05,0.05,.9,.25]);
set(plotter,'color','g') ;

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

%define the DC level
y=y-25;
uicontrol('style','text',...
    'string','DC level', ...
    'fontsize',12, ...
    'position',[x,y,60,20]);
DClevelUI = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[x+60,y,50,20]);

%define the noise level
xt=x+130;
uicontrol('style','text',...
    'string','Noise lvl', ...
    'fontsize',12, ...
    'position',[xt,y,60,20]);
noiselevelUI = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[xt+60,y,50,20]);

%define pulse train
y=y-25;
uicontrol('style','text',...
    'string','Pulse', ...
    'fontsize',12, ...
    'position',[x,y,40,20]);
xt = x +45;
uicontrol('style','text',...
    'string','level', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+40;
pulselevelUI = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
xt = x +45; y=y-25;
uicontrol('style','text',...
    'string','dur', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+40;
pulsedurUI = uicontrol('style','edit', ...
    'string','5', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
xt = x +45; y=y-25;
uicontrol('style','text',...
    'string','space', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+40;
pulsespaceUI = uicontrol('style','edit', ...
    'string','100', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);

%sinewave input
y=y-25;
uicontrol('style','text',...
    'string','Sine', ...
    'fontsize',12, ...
    'position',[x,y,40,20]);
xt = x +45;
uicontrol('style','text',...
    'string','level', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+40;
sinelevelUI = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
xt = xt +55; 
uicontrol('style','text',...
    'string','freq', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+40;
sinefreqUI = uicontrol('style','edit', ...
    'string','0.05', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);

%output input
y=y-25;
uicontrol('style','text',...
    'string','Output', ...
    'fontsize',12, ...
    'position',[x,y,40,20]);
xt = x +45;
uicontrol('style','text',...
    'string','level', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+40;
outputlevelUI = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);

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

%define the DC level
y=y-25;
uicontrol('style','text',...
    'string','DC level', ...
    'fontsize',12, ...
    'position',[x,y,60,20]);
DClevelUI2 = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[x+60,y,50,20]);

%define the noise level
xt=x+130;
uicontrol('style','text',...
    'string','Noise lvl', ...
    'fontsize',12, ...
    'position',[xt,y,60,20]);
noiselevelUI2 = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[xt+60,y,50,20]);

%define pulse train
y=y-25;
uicontrol('style','text',...
    'string','Pulse', ...
    'fontsize',12, ...
    'position',[x,y,40,20]);
xt = x +45;
uicontrol('style','text',...
    'string','level', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+40;
pulselevelUI2 = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
xt = x +45; y=y-25;
uicontrol('style','text',...
    'string','dur', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+40;
pulsedurUI2 = uicontrol('style','edit', ...
    'string','5', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
xt = x +45; y=y-25;
uicontrol('style','text',...
    'string','space', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+40;
pulsespaceUI2 = uicontrol('style','edit', ...
    'string','100', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);

%sinewave input
y=y-25;
uicontrol('style','text',...
    'string','Sine', ...
    'fontsize',12, ...
    'position',[x,y,40,20]);
xt = x +45;
uicontrol('style','text',...
    'string','level', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+40;
sinelevelUI2 = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
xt = xt +55; 
uicontrol('style','text',...
    'string','freq', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+40;
sinefreqUI2 = uicontrol('style','edit', ...
    'string','0.05', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);

%output input
y=y-25;
uicontrol('style','text',...
    'string','Output', ...
    'fontsize',12, ...
    'position',[x,y,40,20]);
xt = x +45;
uicontrol('style','text',...
    'string','level', ...
    'fontsize',12, ...
    'position',[xt,y,40,20]);
xt=xt+40;
outputlevelUI2 = uicontrol('style','edit', ...
    'string','0.0', ...
    'fontsize',12, ...
    'position',[xt,y,50,20]);
%===Synaptic controls=======================

%===integration controls====================
x=290;y=596;
uicontrol('style','text',...
    'string','voltage function', ...
    'fontsize',12, ...
    'backgroundcolor',[1,.8,.8],...
    '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]);
%reset voltage button
y=y-25
uicontrol('style','pushbutton', ...
    'string','Reset Voltage', ...
    'fontsize',12, ...
    'callback','resetV=1',...
    'position',[x+30,y,100,20]);

%show C and R currents
y=y-75;
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-70;
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,.8,.8],...
    '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/linear
x=570; y=560
spikemodeui = uicontrol('style','radiobutton',...
    'string','Spike', ...
    'fontsize',12, ...
    'position',[x,y,70,20], ...
    'value',1,...
    'callback',...
    ['set(linearmodeui,''value'',0);',...
    'set(sigmoidmodeui,''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(linearmodeui,''value'',0);',...
    'set(spikemodeui,''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]);
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);',...
    ]);

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')

%===initial state of buttons
stop= 0; %wait for a quit button push
run = 0; %wait for a draw 
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;
resetV = 0; %start with no reset

%===state variables========================
t=0; dt=1;
vrest = 0.0;
threshold = .8;
v = vrest ;
spike=0;
spiketime = -inf;
lastspiketime = -inf;
isi=[];
%now the two voltage animation objects
set(fig,'currentaxes',Vanimation)
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)

%===output plot variables==================
Vout=vrest*ones(1,200);
Sout=zeros(1,200);
Iout=zeros(1,200);
IRout=zeros(1,200);
ICout=zeros(1,200);
Tout=threshold*ones(1,200);
% and init the plot
set(fig,'currentaxes',plotter)
VP = plot(Vout,'r','linewidth',2);
hold on
SP = plot(Sout+1,'w');
IP = plot(Iout,'b','linewidth',2);
IRP = plot(IRout,'color',[0 1 .5]);
ICP = plot(ICout,'color',[0 .8 1]);
set(IRP,'visible','off')
set(ICP,'visible','off')
TP = plot(Tout,'r:');
set(plotter,'ylim',[-1,2])
set(plotter,'color',0.2*[1 1 1]) ;

%===Now the main simulation loop============
while (stop==0)
    tic
    if (run==1) %running?
        t=t+dt;
        %===draw the voltage animation=========================      
        set(Trect,'position',[-.1,-0.7,0.2,threshold+.2])
        set(Vrect,'position',[-.1,-0.7,0.2,v+.2])
        
        %===update the state variables========================
        Rin = str2num(get(inputRUI,'string'));
        DKrate = str2num(get(DKrateUI,'string'));
        DClevel = str2num(get(DClevelUI,'string'));
        pulselevel  = str2num(get(pulselevelUI,'string'));
        pulsedur = str2num(get(pulsedurUI,'string'));
        pulsespace = str2num(get(pulsespaceUI,'string'));
        noiselevel = str2num(get(noiselevelUI,'string'));
        sinelevel = str2num(get(sinelevelUI,'string'));
        sinefreq = str2num(get(sinefreqUI,'string'));
        outputlevel=str2num(get(outputlevelUI,'string'));
        DClevel2 = str2num(get(DClevelUI2,'string'));
        pulselevel2  = str2num(get(pulselevelUI2,'string'));
        pulsedur2 = str2num(get(pulsedurUI2,'string'));
        pulsespace2 = str2num(get(pulsespaceUI2,'string'));
        noiselevel2 = str2num(get(noiselevelUI2,'string'));
        sinelevel2 = str2num(get(sinelevelUI2,'string'));
        sinefreq2 = str2num(get(sinefreqUI2,'string'));
        outputlevel2=str2num(get(outputlevelUI2,'string'));
        spikemode = get(spikemodeui,'value');
        linearmode = get(linearmodeui,'value');
        sigmoidmode = get(sigmoidmodeui,'value');
        sigmoidN = str2num(get(sigmoidNui,'string'));
        
        %total current
        stimTotal = DClevel + ...
            pulselevel*(mod(t,pulsespace)<=pulsedur)+...
            noiselevel*randn +...
            sinelevel*sin(2*pi*t*sinefreq) + ...
            outputlevel*spike + ...
            DClevel2 + ...
        pulselevel2*(mod(t,pulsespace2)<=pulsedur2)+...
            noiselevel2*randn +...
            sinelevel2*sin(2*pi*t*sinefreq2) + ...
            outputlevel2*spike ;
        
        %update voltage -- DKrate is actually time constant
        v = v + dt/DKrate * (Rin * stimTotal - (v-vrest));
        if resetV==1
            resetV=0;
            v=0;
        end
        
        %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
            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=-.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(1:199) = Vout(2:200);
        Vout(200) = v;
        
        Sout(1:198) = Sout(2:199);
        Sout(199) = spike;
        Sout(200) = 0;
        
        Iout(1:199) = Iout(2:200);
        Iout(200) = stimTotal;
        IRout(1:199) = IRout(2:200);
        IRout(200) = (v-vrest)/Rin;
        ICout(1:199) = ICout(2:200);
        ICout(200) = stimTotal-IRout(200);
        
        Tout(1:199) = Tout(2:200);
        Tout(200) = threshold;
        
        set(VP,'ydata',Vout);
        set(SP,'ydata',Sout+1);
        set(IP,'ydata',Iout);
        set(IRP,'ydata',IRout);
        set(ICP,'ydata',ICout);
        set(TP,'ydata',Tout);
        if showcurrent==1
            set(IRP,'visible','on')
            set(ICP,'visible','on')
        else
            set(IRP,'visible','off')
            set(ICP,'visible','off')
        end
    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




