function [J] = jhujed_obj(theta,uspks,Topt,dt,Nita)
be=theta(1);
bi=theta(2);
we=theta(3);
wi=theta(4);
wee=theta(5);
wei=theta(6);
wie=theta(7);
wii=theta(8);
Fmaxe=theta(9);
Fmaxi=theta(10);
vce=theta(11);
vci=theta(12);
ae=theta(13);
ai=theta(14);

Nopt=Topt/dt+1;
ta=0:dt:Topt;
rea=zeros(1,Nopt);
llkl=zeros(1,Nita);

for oi=1:Nita
    Tspks=uspks(oi).spks;
    AE=uspks(oi).A;
    WE=uspks(oi).W;
    PE=uspks(oi).P;
    ve=0;
    vi=0;
    rei=0;
    for i=1:Nopt
        t=ta(i);
        u=sum(AE.*cos(WE.*t+PE));
        ge=Fmaxe/(exp(ae*(vce - ve)) + 1);
        gev=(Fmaxe*ae*exp(ae*(vce - ve)))/(exp(ae*(vce - ve)) + 1)^2;
        gevv=(2*Fmaxe*ae^2*exp(2*ae*(vce - ve)))/(exp(ae*(vce - ve)) + 1)^3 - (Fmaxe*ae^2*exp(ae*(vce - ve)))/(exp(ae*(vce - ve)) + 1)^2;
        gi=Fmaxi/(exp(ai*(vci - vi)) + 1);
        giv=(Fmaxi*ai*exp(ai*(vci - vi)))/(exp(ai*(vci - vi)) + 1)^2;
        givv=(2*Fmaxi*ai^2*exp(2*ai*(vci - vi)))/(exp(ai*(vci - vi)) + 1)^3 - (Fmaxi*ai^2*exp(ai*(vci - vi)))/(exp(ai*(vci - vi)) + 1)^2;
        vedot=be*(-ve+ge*wee-gi*wei+we*u);
        vidot=bi*(-vi-gi*wii+ge*wie+wi*u);
        re=ge;
        rea(i)=re;
        ve=ve+vedot*dt;
        vi=vi+vidot*dt;
        rei=rei+re*dt;
    end
    reas=interp1(ta,rea,Tspks);
    lreas=log(reas);
    sreas=sum(lreas);
    llkl(oi)=-rei+sreas;
end
J=-sum(llkl);