# How to fit this neuron firing model with R?

I originally posted this as an answer elsewhere but in retrospect it seems more like a question: What is the sample-size range for which the median should be preferred to the mean as a measure of central tendency? And why?

I have R code that simulated neuron activity based on the model described at these links: http://www.izhikevich.org/publications/spikes.htm

http://www.izhikevich.org/publications/spikes.pdf

I know how to simulate one but not fit the parameters to existing data. Take this pattern which was generated with the code below. Each neuron could have different a, b, c, d, k1, k2, k3 parameters and also the I parameter (Input Current) could vary over time or be constant for each neuron. The model is described in the pdf above. Anyone know how to fit this model using R? What I would like to be able to do is fit a set of parameters to actual firing data (which would be voltage over time, same as what is simulated). R code to simulate Neuron:

a=.1; b=.2; c=-65; d=2
k1=.04; k2=5; k3=140
I=4

Ntime=1000
v<-matrix(nrow=Ntime)
u<-matrix(nrow=Ntime)

v=-65
u=0
for(i in 1:Ntime){
#Fire Neuron if Voltage>30
if(v[i]>=30){
v[i]<-c
u[i]<-u[i]+d
}

#Calculate Voltage Change
dv= k1*v[i]^2 + k2*v[i] + k3 - u[i] + I
du= a*(b*v[i] -u[i])

#Update Voltage
v[i+1]<-v[i]+dv
u[i+1]<-u[i]+du

}

plot(v, type="l", xlab="Time (ms)", ylab="Voltage")

• What mode of neuron firing are you trying to fit? Does your experimental data show discrete firing (as the model) or does it display burst firing, trains, etc?
– user32490
Apr 1, 2014 at 21:05
• @leonardo I do not actually have experimental data, but I am wondering how to accomplish this. Pretend that the simulated data output by the code provided is the real data. Apr 1, 2014 at 21:10
• In that case, I was going to suggest trying to divide the parameter space into sections based on the observed firing pattern. One approach may be to get a sense of periodicity from the spikes (cropped above a threshold, say of 30 mV). A simple first approach walking through the parameters may be to minimize the root-mean square error.
– user32490
Apr 1, 2014 at 21:15
• @leonardo I am not clear on what you mean by "divide the parameter space into sections" or "walking through the parameters". Minimizing the rmse sounds fine as a goal, but how to achieve it? This seems like it may be a problem good for gradient descent, but I am not sure how I would code it. Apr 1, 2014 at 21:22
• Dividing parameter space, refering to the top-right panels in the web link, is a way of only checking subsets of parameter values if you knew the mode(s) of firing a priori. Walking through the parameters would be to define a range of each parameter and the running the simulation with incrementally small steps through the range.
– user32490
Apr 1, 2014 at 22:18