function nextState = lp_step_template(state,input,noise,param)

steady_n  = 1/(1+exp((stateV-hVn)/sn));
relax_n   = cn/(1+exp((stateV-hVkn)/skn));

steady_ao = 1/(1+exp((stateV-hVao1+cf*stateCa)/sao1))*1/(1+exp((stateV-hVao2+ ...
    cf*stateCa)/sao2))*(stateCa/(c1+stateCa));
relax_ao  = koa;
steady_bo = c2/(c3+stateCa);
relax_bo  = kob;


steady_aA = 1/(1+exp((stateV-hVA1)/sA1));
relax_aA  = kA1;
steady_bA1= 1/(1+exp((stateV-hVB)/sB));
relax_bA1 = kB;
steady_bA2= 1/(1+exp((stateV-hVA2)/sA2));
relax_bA2 = cA2/(1+exp((stateV-hVA2)/sA2));
weightX   = 1/(1+exp((stateV-hVX)/sX));

steady_aCa1 = 1/(1+exp((stateV-hVaCa1)/saCa1));
relax_aCa1  = kaCa1;
steady_bCa1 = 1/(1+exp((stateV-hVbCa1)/sbCa1));
relax_bCa1  = kbCa1;
steady_aCa2 = 1/(1+exp((stateV-hVaCa2)/saCa2));
relax_aCa2  = kaCa2;
ECa    = 8.314472*283/(2*96485.3399)*log(13000/stateCa)*1000;

steady_r = 1/(1+exp((stateV-hVr)/sr));
relax_r  = cr*(1+exp((stateV-hVkr)/skr));

alpha_m = (-4 - 0.1*stateV)./(exp(-4 - 0.1*stateV) - 1);
alpha_h = 0.0027*exp(-stateV./20);
beta_m = 0.1081*exp(-stateV./18);
beta_h = 1./(exp(-3.5-0.1*stateV))+1);
alpha_n = (-0.55 - 0.01*stateV)./(exp(-5.5 - 0.1*stateV) - 1);
beta_n = 0.05547*exp(-stateV./80);

%alpha_m  = cam*(stateV-hVam)/(1-exp((stateV-hVam)/sam));
%beta_m  = cbm*exp((stateV-hVbm)/sbm);
steady_m = alpha_m/(alpha_m+beta_m);
relax_m  = alpha_m+beta_m;
%alpha_h  = cah*exp((stateV-hVah)/sah);
%beta_h  = 1/(1+exp((stateV-hVbh)/sbh));
steady_h = alpha_h/(alpha_h+beta_h);
relax_h  = alpha_h+beta_h;


id   = gd*stateN^4*(stateV+65-EK);
ioCa = goCa*stateAo*stateBo*(stateV-EK);
iA   = gA*stateAA^3*(weightX*stateBA1+(1-weightX)*stateBA2)*(stateV-EK);
iCa  = (gCa1*stateACa1*stateBCa1+gCa2*stateACa2)*(stateV-ECa);
ih   = gH*stateR*(stateV-Eh);
iNa  = gNa*stateM^3*stateH*(stateV+65-ENa);
il   = gl*(stateV+65-El);

ionicCurrent = id + ioCa + iA + iCa + ih + iNa + il;
%ionicCurrent = id + iNa + il;

nextState = [stateV + oneDT./cm*(injectedCurrent - ionicCurrent) ...
                    + sqrtDT*PNoiseLevel*PNoise, ...
		 ... % stateN + oneDT*(steady_n-stateN)*relax_n, ...
         stateN + oneDT*(alpha_n*(1-stateN) - beta_n*stateN), ...
		 stateCa + oneDT*(-cICa*iCa-kCa*stateCa+kCa*CaO), ...
		 stateAo + oneDT*(steady_ao-stateAo)*relax_ao, ...
		 stateBo + oneDT*(steady_bo-stateBo)*relax_bo, ...
         stateAA + oneDT*(steady_aA-stateAA)*relax_aA, ...
		 stateBA1 + oneDT*(steady_bA1-stateBA1)*relax_bA1, ...
		 stateBA2 + oneDT*(steady_bA2-stateBA2)*relax_bA2, ...
		 stateACa1 + oneDT*(steady_aCa1-stateACa1)*relax_aCa1, ...
		 stateBCa1 + oneDT*(steady_bCa1-stateBCa1)*relax_bCa1, ...
		 stateACa2 + oneDT*(steady_aCa2-stateACa2)*relax_aCa2, ...
		 stateR + oneDT*(steady_r-stateR)*relax_r, ...
		 ... %stateM + oneDT*(steady_m-stateM)*relax_m, ...
		 ... %stateH + oneDT*(steady_h-stateH)*relax_h];
         stateM + oneDT*(alpha_m*(1-stateM) - beta_m*stateM), ...
         stateH + oneDT*(alpha_h*(1-stateH) - beta_h*stateH)];
