%two pairs of oscillator coupled via inhibtion

function [exit1, oute1,inhib1,outi1, exit2, oute2, inhib2, outi2] = oscillate6 (input1, input2, taue, taui, threshe, threshi, steps, we, wi, wii)
exit1(1)=0; inhib1(1)=0;old_ex1=0;old_in1=0;
oute1(1)=0.0;outi1(1)=0.0;
exit2(1)=0; inhib2(1)=0;old_ex2=0;old_in2=0;
oute2(1)=0.0;outi2(1)=0.0;

for i=2:steps
exit1(i)=exp(-1/taue)*old_ex1+(1-exp(-1/taue))*(input1(i)-wi*outi1(i-1));
if (exit1(i)>threshe)
oute1(i)=1.0;
old_ex1=0.0;
else
oute1(i)=0.0;
old_ex1=exit1(i);
end;

inhib1(i)=exp(-1/taui)*old_in1+(1-exp(-1/taui))*(input1(i)+we*oute1(i-1)-wii*outi2(i-1));
if (inhib1(i)>threshi)
outi1(i)=1.0;
old_in1=0.0;
else
outi1(i)=0.0;
old_in1=inhib1(i);
end;
exit2(i)=exp(-1/taue)*old_ex2+(1-exp(-1/taue))*(input2(i)-wi*outi2(i-1));
if (exit2(i)>threshe)
oute2(i)=1.0;
old_ex2=0.0;
else
oute2(i)=0.0;
old_ex2=exit2(i);
end;

inhib2(i)=exp(-1/taui)*old_in2+(1-exp(-1/taui))*(input2(i)+we*oute2(i-1)-wii*oute1(i-1));
if (inhib2(i)>threshi)
outi2(i)=1.0;
old_in2=0.0;
else
outi2(i)=0.0;
old_in2=inhib2(i);
end;
end;

t=linspace(1,100);
for (i=1:100)
input(i)=1.0;
end;

[e1,oute1, in1 ,outi1, e2, oute2, in2, outi2]=oscillate6(input, 1.2*input, 5, 5, 0.5, 0.8, 100, 1.0, 1.0, 1.0);

figure
subplot (2,1,1);
plot (t,e1+5*oute1);
hold;

plot(t, in1+5*outi1, 'r');

subplot (2,1,2);

plot (t,e2+5*oute2);
hold;

plot(t, in2+5*outi2,