function xdot=mlcdfun(t,x,flag,tdel);

% Carol Lucas
% simple model for coincidence detection.
%
% ML cell receives two just-subthreshold inputs that are identical and
% excitatory -  tdel represents timing difference between them.
%
% this ML model uses the standard ML (Type II) params from Chapt 2 of this
% book (same as Rinzel/Ermentrout in  Koch and Segev's book.
%
% Note: treat tdel as a variable (with ode: tdel'=0) so can do Poincare
% map for response tuning curve
%  init v1=-60.9;w1=.0149;s2=0.;s1=0.;tdel=0;


v1=x(1,1);
w1=x(2,1);
s2=x(3,1);
s1=x(4,1);

t0=100;teon=5;vrest=-60;
gsyne=1.5;vsyne=100;
vk=-84;vl=-60;vca=120;
i=0;gk=8;gl=2;c=20;
va=-1.2;vb=18;vc=-5;vd=30;phi=.04;gca=4.;
thetasyn=20;ksyn=2;alphae=1;betae=0.3;

% fix the start time of excitation from R2 and delay the
% start of excitation from R1  by tdel (if tdel<0; R1 input precedes
% R2 input)

vr2=100*(t>=t0)*(t<=t0+teon)+vrest;
vr1=100*(t>=t0+tdel)*(t<=t0+tdel+teon)+vrest;

isyn2=gsyne*s2*(vsyne-v1);
isyn1=gsyne*s1*(vsyne-v1);

minf=.5*(1+tanh((v1-va)/vb));
winf=.5*(1+tanh((v1-vc)/vd));

tauw=1/cosh((v1-vc)/(2*vd));
sinf1=1/(1+exp(-(vr1-thetasyn)/ksyn));
sinf2=1/(1+exp(-(vr2-thetasyn)/ksyn));
%
xdot(1,1) = (i-gca*minf*(v1-vca)-gk*w1*(v1-vk)-gl*(v1-vl)+isyn2+isyn1)/c;
xdot(2,1) = phi*(winf-w1)/tauw;
xdot(3,1) = alphae*sinf2*(1-s2)-betae*s2;
xdot(4,1) = alphae*sinf1*(1-s1)-betae*s1;