function simulation(network,para,iinp,trial_num,hd_txt,OUTPUT_OPT)

%DO NOT EDIT THIS FILE, UNLESS YOU ARE ABOSLUTE SURE WHAT YOU ARE DOING.
%This is the template file use to generate the ODE simulation .m file
%"simulation.m"
%Written by Ning Jiang, Institute of Biomedical Engineering, Univesity of New
%Brunswick, NB, Canada, E3B 5A3
%Email: ning.jiang@unb.ca
%Created Date: Nov 19, 2006
%Rev. 1.1, Nov 20, 2006 --adding para to the saved variables

global inp;
duration=para.duration*1000; 
neuron_num=size(iinp,2);
if OUTPUT_OPT
    out_v=zeros(1,neuron_num);
    out_t=0;
end
initial=zeros(1,3*neuron_num);
initial(1:3:3*neuron_num-2)=rand(1,neuron_num)*5-60;
initial(2:3:3*neuron_num-1)=rand(1,neuron_num);
initial(3:3:3*neuron_num)=rand(1,neuron_num)*0.5;
options=odeset('InitialStep',para.stepsize);
firing_t=cell(1,neuron_num);
v=initial(1:3:3*neuron_num-2);
t_course=0;
t0=0; 
for n=1:size(iinp,1)
    inp=iinp(n,:);
    [t,y]=ode45(@hhn,[t0,t0+para.stepsize],initial,options);
    initial=y(size(y,1),:);
    v=[v;y((2:size(y,1)),(1:3:neuron_num*3-2))];
    t_course=[t_course;t(2:length(t))];
    if rem(n,1000/para.stepsize) == 0
        if hd_txt
            set(hd_txt,'string',['trial ',num2str(trial_num),':  ',num2str(duration/1000-n*para.stepsize/1000),' sec left']);
            pause(0.01);
        else
            disp(['trial ',num2str(trial_num),':  ',num2str(duration/1000-n*para.stepsize/1000),' sec left']);
        end
        for nn=1:neuron_num
            [ft,abnormal]=get_interval(v(:,nn),t_course,1);
            if ft ~=0
                firing_t{nn}=[firing_t{nn},ft];
            end
        end
        if OUTPUT_OPT
            out_v=[out_v;v(2:size(v,1),:)];
            out_t=[out_t;t_course(2:length(t_course))];
        end
        v=v(size(v,1),:);
        t_course=t_course(length(t_course));
    end
    t0=t(length(t));
end
if OUTPUT_OPT
    out_t=out_t(2:length(out_t));
    out_v=out_v(2:size(out_v,1),:);
    save([para.resultfile,'_',num2str(trial_num)],'firing_t','out_t','out_v','iinp','network','para');
else
    save([para.resultfile,'_',num2str(trial_num)],'firing_t','iinp','network','para');
end

function dydt=hhn(t,y)
global inp;dydt=[-(2.25*(y(1)--60)+37.5/(1+exp(-(y(1)+30)/15))^3*(0.5-y(2))*(y(1)-80)+45*y(2)^4*(y(1)--60))+inp(1)-0.05*y(63)*(y(1)-0);
0.1*(1/(1+exp((y(1)--32)/-4))-y(2))/(2+4/(1+exp((y(1)--40)/2)));
2*(1-y(3))*(1/(1+exp(-(y(1)--37)/2)))-0.03*y(3);


-(2.25*(y(4)--60)+37.5/(1+exp(-(y(4)+30)/15))^3*(0.5-y(5))*(y(4)-73.3333)+45*y(5)^4*(y(4)--60))+inp(2)-0.05*y(63)*(y(4)-0);
0.1*(1/(1+exp((y(4)--32)/-4))-y(5))/(2+4/(1+exp((y(4)--40)/2)));
2*(1-y(6))*(1/(1+exp(-(y(4)--37)/2)))-0.03*y(6);


-(2.25*(y(7)--60)+37.5/(1+exp(-(y(7)+30)/15))^3*(0.5-y(8))*(y(7)-66.6667)+45*y(8)^4*(y(7)--60))+inp(3)-0.05*y(63)*(y(7)-0);
0.1*(1/(1+exp((y(7)--32)/-4))-y(8))/(2+4/(1+exp((y(7)--40)/2)));
2*(1-y(9))*(1/(1+exp(-(y(7)--37)/2)))-0.03*y(9);


-(2.25*(y(10)--60)+37.5/(1+exp(-(y(10)+30)/15))^3*(0.5-y(11))*(y(10)-60)+45*y(11)^4*(y(10)--60))+inp(4)-0.05*y(63)*(y(10)-0);
0.1*(1/(1+exp((y(10)--32)/-4))-y(11))/(2+4/(1+exp((y(10)--40)/2)));
2*(1-y(12))*(1/(1+exp(-(y(10)--37)/2)))-0.03*y(12);


-(2.25*(y(13)--60)+37.5/(1+exp(-(y(13)+30)/15))^3*(0.5-y(14))*(y(13)-53.3333)+45*y(14)^4*(y(13)--60))+inp(5)-0.05*y(63)*(y(13)-0);
0.1*(1/(1+exp((y(13)--32)/-4))-y(14))/(2+4/(1+exp((y(13)--40)/2)));
2*(1-y(15))*(1/(1+exp(-(y(13)--37)/2)))-0.03*y(15);


-(2.25*(y(16)--60)+37.5/(1+exp(-(y(16)+30)/15))^3*(0.5-y(17))*(y(16)-46.6667)+45*y(17)^4*(y(16)--60))+inp(6)-0.05*y(63)*(y(16)-0);
0.1*(1/(1+exp((y(16)--32)/-4))-y(17))/(2+4/(1+exp((y(16)--40)/2)));
2*(1-y(18))*(1/(1+exp(-(y(16)--37)/2)))-0.03*y(18);


-(2.25*(y(19)--60)+37.5/(1+exp(-(y(19)+30)/15))^3*(0.5-y(20))*(y(19)-40)+45*y(20)^4*(y(19)--60))+inp(7)-0.05*y(63)*(y(19)-0);
0.1*(1/(1+exp((y(19)--32)/-4))-y(20))/(2+4/(1+exp((y(19)--40)/2)));
2*(1-y(21))*(1/(1+exp(-(y(19)--37)/2)))-0.03*y(21);


-(2.25*(y(22)--60)+37.5/(1+exp(-(y(22)+30)/15))^3*(0.5-y(23))*(y(22)-33.3333)+45*y(23)^4*(y(22)--60))+inp(8)-0.05*y(63)*(y(22)-0);
0.1*(1/(1+exp((y(22)--32)/-4))-y(23))/(2+4/(1+exp((y(22)--40)/2)));
2*(1-y(24))*(1/(1+exp(-(y(22)--37)/2)))-0.03*y(24);


-(2.25*(y(25)--60)+37.5/(1+exp(-(y(25)+30)/15))^3*(0.5-y(26))*(y(25)-26.6667)+45*y(26)^4*(y(25)--60))+inp(9)-0.05*y(63)*(y(25)-0);
0.1*(1/(1+exp((y(25)--32)/-4))-y(26))/(2+4/(1+exp((y(25)--40)/2)));
2*(1-y(27))*(1/(1+exp(-(y(25)--37)/2)))-0.03*y(27);


-(2.25*(y(28)--60)+37.5/(1+exp(-(y(28)+30)/15))^3*(0.5-y(29))*(y(28)-20)+45*y(29)^4*(y(28)--60))+inp(10)-0.05*y(63)*(y(28)-0);
0.1*(1/(1+exp((y(28)--32)/-4))-y(29))/(2+4/(1+exp((y(28)--40)/2)));
2*(1-y(30))*(1/(1+exp(-(y(28)--37)/2)))-0.03*y(30);


-(2.25*(y(31)--60)+37.5/(1+exp(-(y(31)+30)/15))^3*(0.5-y(32))*(y(31)-80)+45*y(32)^4*(y(31)--60))+inp(11)-0.05*y(63)*(y(31)-0);
0.1*(1/(1+exp((y(31)--32)/-4))-y(32))/(2+4/(1+exp((y(31)--40)/2)));
2*(1-y(33))*(1/(1+exp(-(y(31)--37)/2)))-0.03*y(33);


-(2.25*(y(34)--60)+37.5/(1+exp(-(y(34)+30)/15))^3*(0.5-y(35))*(y(34)-73.3333)+45*y(35)^4*(y(34)--60))+inp(12)-0.05*y(63)*(y(34)-0);
0.1*(1/(1+exp((y(34)--32)/-4))-y(35))/(2+4/(1+exp((y(34)--40)/2)));
2*(1-y(36))*(1/(1+exp(-(y(34)--37)/2)))-0.03*y(36);


-(2.25*(y(37)--60)+37.5/(1+exp(-(y(37)+30)/15))^3*(0.5-y(38))*(y(37)-66.6667)+45*y(38)^4*(y(37)--60))+inp(13)-0.05*y(63)*(y(37)-0);
0.1*(1/(1+exp((y(37)--32)/-4))-y(38))/(2+4/(1+exp((y(37)--40)/2)));
2*(1-y(39))*(1/(1+exp(-(y(37)--37)/2)))-0.03*y(39);


-(2.25*(y(40)--60)+37.5/(1+exp(-(y(40)+30)/15))^3*(0.5-y(41))*(y(40)-60)+45*y(41)^4*(y(40)--60))+inp(14)-0.05*y(63)*(y(40)-0);
0.1*(1/(1+exp((y(40)--32)/-4))-y(41))/(2+4/(1+exp((y(40)--40)/2)));
2*(1-y(42))*(1/(1+exp(-(y(40)--37)/2)))-0.03*y(42);


-(2.25*(y(43)--60)+37.5/(1+exp(-(y(43)+30)/15))^3*(0.5-y(44))*(y(43)-53.3333)+45*y(44)^4*(y(43)--60))+inp(15)-0.05*y(63)*(y(43)-0);
0.1*(1/(1+exp((y(43)--32)/-4))-y(44))/(2+4/(1+exp((y(43)--40)/2)));
2*(1-y(45))*(1/(1+exp(-(y(43)--37)/2)))-0.03*y(45);


-(2.25*(y(46)--60)+37.5/(1+exp(-(y(46)+30)/15))^3*(0.5-y(47))*(y(46)-46.6667)+45*y(47)^4*(y(46)--60))+inp(16)-0.05*y(63)*(y(46)-0);
0.1*(1/(1+exp((y(46)--32)/-4))-y(47))/(2+4/(1+exp((y(46)--40)/2)));
2*(1-y(48))*(1/(1+exp(-(y(46)--37)/2)))-0.03*y(48);


-(2.25*(y(49)--60)+37.5/(1+exp(-(y(49)+30)/15))^3*(0.5-y(50))*(y(49)-40)+45*y(50)^4*(y(49)--60))+inp(17)-0.05*y(63)*(y(49)-0);
0.1*(1/(1+exp((y(49)--32)/-4))-y(50))/(2+4/(1+exp((y(49)--40)/2)));
2*(1-y(51))*(1/(1+exp(-(y(49)--37)/2)))-0.03*y(51);


-(2.25*(y(52)--60)+37.5/(1+exp(-(y(52)+30)/15))^3*(0.5-y(53))*(y(52)-33.3333)+45*y(53)^4*(y(52)--60))+inp(18)-0.05*y(63)*(y(52)-0);
0.1*(1/(1+exp((y(52)--32)/-4))-y(53))/(2+4/(1+exp((y(52)--40)/2)));
2*(1-y(54))*(1/(1+exp(-(y(52)--37)/2)))-0.03*y(54);


-(2.25*(y(55)--60)+37.5/(1+exp(-(y(55)+30)/15))^3*(0.5-y(56))*(y(55)-26.6667)+45*y(56)^4*(y(55)--60))+inp(19)-0.05*y(63)*(y(55)-0);
0.1*(1/(1+exp((y(55)--32)/-4))-y(56))/(2+4/(1+exp((y(55)--40)/2)));
2*(1-y(57))*(1/(1+exp(-(y(55)--37)/2)))-0.03*y(57);


-(2.25*(y(58)--60)+37.5/(1+exp(-(y(58)+30)/15))^3*(0.5-y(59))*(y(58)-20)+45*y(59)^4*(y(58)--60))+inp(20)-0.05*y(63)*(y(58)-0);
0.1*(1/(1+exp((y(58)--32)/-4))-y(59))/(2+4/(1+exp((y(58)--40)/2)));
2*(1-y(60))*(1/(1+exp(-(y(58)--37)/2)))-0.03*y(60);


-(2.25*(y(61)--60)+37.5/(1+exp(-(y(61)+30)/15))^3*(0.5-y(62))*(y(61)-55)+45*y(62)^4*(y(61)--60))+inp(21);
0.1*(1/(1+exp((y(61)--32)/-4))-y(62))/(2+4/(1+exp((y(61)--40)/2)));
2*(1-y(63))*(1/(1+exp(-(y(61)--37)/2)))-0.03*y(63);


];