%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A Spatially Extended Model for Macroscopic Spike-Wave Discharges %
% Peter N. Taylor & Gerold Baier                                   %
% peter.taylor@postgrad.manchester.ac.uk                           %
% J. Comp. Neurosci. 2011                                          %
% ---------------------------------------------------------------- %
% This code reproduces the SWD in figure 3b of the manuscript.     %
% Usage: At the MATLAB command window type >> taylor_baier_one_dim %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function taylor_baier_one_dim
h1=5;
h2=-0.9;
h3=-1.6;

boundary=25;
interval=0.5;
rangelength=((2*boundary)/interval)+1; % This is equivlent to 'n' in the manuscript (101)
w1matrix=[];w2matrix=[];w3matrix=[];
xloc=[-boundary:interval:boundary];

% Create the mexican hat connectivities (requires the file 'mexhat.m')...
for i=1:rangelength
       [ty,tmp]=mexhat(xloc(i), boundary, interval,0.5,6);      %1.8   0.36
       w1matrix=horzcat(w1matrix,tmp);
       [ty,tmp]=mexhat(xloc(i), boundary, interval,0.9,1.25);      %1     0.4
       w2matrix=horzcat(w2matrix,tmp);
       [ty,tmp]=mexhat(xloc(i), boundary, interval,2,3);      %3.3   1.32
       w3matrix=horzcat(w3matrix,tmp);
end
    
w1matrix = triu(w1matrix);
w1matrix=w1matrix+w1matrix';

w2matrix = triu(w2matrix);
w2matrix=w2matrix+w2matrix';

w3matrix = triu(w3matrix);
w3matrix=w3matrix+w3matrix';
for i=1:rangelength
    w1matrix(i,i)=w1matrix(i,i)/2;
    w2matrix(i,i)=w2matrix(i,i)/2;
    w3matrix(i,i)=w3matrix(i,i)/2;
end

% Some initial conditions...
init=[-0.974581960533541,-0.974581956287244,-0.974581954315587,-0.974581955298441,-0.974581958890720,-0.974581963838189,-0.974581968451107,-0.974581971242544,-0.974581971467674,-0.974581969346310,-0.974581965892685,-0.974581962447881,-0.974581960132444,-0.974581959455298,-0.974581960226297,-0.974581961770732,-0.974581963309053,-0.974581964306854,-0.974581964640516,-0.974581964534348,-0.974581964343607,-0.974581964323743,-0.974581964510558,-0.974581964755552,-0.974581964866178,-0.974581964746343,-0.974581964446155,-0.974581964098398,-0.974581963797328,-0.974581963514208,-0.974581963119321,-0.974581962507439,-0.974581961748035,-0.974581961151168,-0.974581961177393,-0.974581962207563,-0.974581964278226,-0.974581966928299,-0.974581969266354,-0.974581970268874,-0.974581969205452,-0.974581966013329,-0.974581961447928,-0.974581956918248,-0.974581954041228,-0.974581954064620,-0.974581957367235,-0.974581963226810,-0.974581969957882,-0.974581975397315,-0.974581977595596,-0.974581975495763,-0.974581969372290,-0.974581960864219,-0.974581952555643,-0.974581947200746,-0.974581946818103,-0.974581951946683,-0.974581961333275,-0.974581972202400,-0.974581981072921,-0.974581984888201,-0.974581982090067,-0.974581973250670,-0.974581961001767,-0.974581949236569,-0.974581941825615,-0.974581941284294,-0.974581947872719,-0.974581959470150,-0.974581972289204,-0.974581982182189,-0.974581986065999,-0.974581982945674,-0.974581974171606,-0.974581962860878,-0.974581952729998,-0.974581946797138,-0.974581946435165,-0.974581951090783,-0.974581958706904,-0.974581966614759,-0.974581972508709,-0.974581975132758,-0.974581974470614,-0.974581971457699,-0.974581967421480,-0.974581963530685,-0.974581960476564,-0.974581958462528,-0.974581957419842,-0.974581957271395,-0.974581958071321,-0.974581959942566,-0.974581962863435,-0.974581966449821,-0.974581969893702,-0.974581972143369,-0.974581972283610,-0.974581969958055,-0.974581965630378,-0.299671199193836,-0.299671219169071,-0.299671255895606,-0.299671303327744,-0.299671352064499,-0.299671391277294,-0.299671411346523,-0.299671406475997,-0.299671376474743,-0.299671327105425,-0.299671268833705,-0.299671214316225,-0.299671175350141,-0.299671160138802,-0.299671171574174,-0.299671206877660,-0.299671258527585,-0.299671316083593,-0.299671368383257,-0.299671405631286,-0.299671421052083,-0.299671411936740,-0.299671380021327,-0.299671331180059,-0.299671274447674,-0.299671220452509,-0.299671179462907,-0.299671159390238,-0.299671164178577,-0.299671192977573,-0.299671240321914,-0.299671297277300,-0.299671353255219,-0.299671398046026,-0.299671423622188,-0.299671425400712,-0.299671402847803,-0.299671359467806,-0.299671302285920,-0.299671240914587,-0.299671186244698,-0.299671148796738,-0.299671136845943,-0.299671154582849,-0.299671200716750,-0.299671267983634,-0.299671343916534,-0.299671412967061,-0.299671459694302,-0.299671472372237,-0.299671446133680,-0.299671384762100,-0.299671300495323,-0.299671211674157,-0.299671138642481,-0.299671098828875,-0.299671102255635,-0.299671148713637,-0.299671227477076,-0.299671319780569,-0.299671403509171,-0.299671458883831,-0.299671473577324,-0.299671445803727,-0.299671384486090,-0.299671306468056,-0.299671231629577,-0.299671177400009,-0.299671154316064,-0.299671163889052,-0.299671199256434,-0.299671248174666,-0.299671297187440,-0.299671335523874,-0.299671357522752,-0.299671363021703,-0.299671355929000,-0.299671341808370,-0.299671325537767,-0.299671309901439,-0.299671295467012,-0.299671281522239,-0.299671267441262,-0.299671253762240,-0.299671242484871,-0.299671236503376,-0.299671238483181,-0.299671249702684,-0.299671269347069,-0.299671294506827,-0.299671320827579,-0.299671343522315,-0.299671358382828,-0.299671362517856,-0.299671354727306,-0.299671335587158,-0.299671307384020,-0.299671273981909,-0.299671240578605,-0.299671213208409,-0.299671197856504,0.114337528798895,0.114350962578483,0.114344169711404,0.114348879146138,0.114335784610123,0.114333920231971,0.114356273551021,0.114354875519427,0.114361439926507,0.114334706481384,0.114338352748343,0.114336540324029,0.114350441817637,0.114340874497322,0.114364620562490,0.114337511002667,0.114335148080898,0.114351810183701,0.114358430642211,0.114342361028852,0.114337532405083,0.114342766098931,0.114338472988715,0.114348767588854,0.114362300432610,0.114353049123442,0.114335155637631,0.114344458573234,0.114332331403704,0.114346299596592,0.114342710565265,0.114360408473891,0.114352835693209,0.114340744042751,0.114353670029776,0.114327765659931,0.114336187854793,0.114354528273218,0.114354670079459,0.114348432270077,0.114345744839908,0.114336898347058,0.114357097290267,0.114329725722601,0.114349829905968,0.114346579616739,0.114352377979337,0.114337426629117,0.114349116680488,0.114352173087403,0.114335336546106,0.114343850769992,0.114359070791844,0.114345346535524,0.114338576255530,0.114348528181529,0.114339630314795,0.114352370793400,0.114339957999708,0.114341609644916,0.114347561709873,0.114356341987156,0.114334539185591,0.114352069395830,0.114349578134634,0.114357926473885,0.114328214056881,0.114340441929071,0.114351840257831,0.114339819241775,0.114338592437362,0.114355022088257,0.114352029068425,0.114350660099676,0.114352748150017,0.114331755196924,0.114341583460491,0.114344758034206,0.114337443349991,0.114353275432103,0.114357983717876,0.114338567182118,0.114353736185331,0.114333716107056,0.114357170159537,0.114333429504266,0.114348268114906,0.114340383700828,0.114354703628718,0.114343970287390,0.114343099149274,0.114355815685634,0.114337994299381,0.114341313527925,0.114358798442439,0.114328146354512,0.114359940513503,0.114333983614518,0.114350478195323,0.114348050255283,0.114351209285196;];

% run the ODE...
options = odeset('AbsTol',0.00001,'RelTol',0.00001); %set options (error tolerance)
tic;
[t,u] = ode23(@(t,u)ode(t,u,boundary,interval,w1matrix, w2matrix, w3matrix,rangelength,h1, h2, h3),[0:0.1:100],init',options);
toc;

% Plot results...
figure(1);subplot(2,1,1);plot(mean(u'),'k');axis('tight');xlabel('Time');ylabel('Mean all E, I_1, I_2');
subplot(2,1,2);imagesc(u(:,1:end/3)');ylabel('Spatial locations in E');xlabel('Time');


function dudt=ode(t,u,boundary,interval,w1matrix, w2matrix,w3matrix, rangelength,h1, h2, h3)
    w4=4;
    w5=3;

    tau1=0.66;
    tau2=200;

    e1=u(1:rangelength);
    i1=u(rangelength+1:(2*rangelength));
    i2=u(2*rangelength+1:(3*rangelength));

    e1sigmatrix=sigfunc(e1);
    i1sigmatrix=sigfunc(i1);
    i2sigmatrix=sigfunc(i2);

    excitatory=sum(w1matrix(:,:).*repmat(e1sigmatrix,1,rangelength))';
    inhibitory1=sum(w2matrix(:,:).*repmat(i1sigmatrix,1,rangelength))';
    inhibitory2=sum(w3matrix(:,:).*repmat(i2sigmatrix,1,rangelength))';    
    
    % Eqn (2) from the manuscript
    du1dt =  h1 - e1 + (excitatory) - (inhibitory1) - (inhibitory2)  ;
    du2dt = (h2 - i1 + e1sigmatrix*w4) / tau1;
    du3dt = (h3 - i2 + e1sigmatrix*w5) / tau2;
    
    dudt=vertcat(du1dt,du2dt,du3dt);      
    function u=sigfunc(u)
%         beg=-0.3;
%         en=0.3;
%         u=(u-beg).*(1/(en-beg));
%         u(u>1)=1;
%         u(u<0)=0;

         u=(1./(1+5000.^-(u))); %comment this out and uncomment the 
%          preceding lines to use PWL rather than sigmoid
    end
end
end