Cyclical Hidden Markov models (HMMs) can be used to sort spikes. For example, in Herbst et al. (2008), wideband data is assumed to be generated by an HMM, where for most of the time the HMM is in the rest state, and with small probability jumps to the spike initiation state. Once it’s in the spike initiation state, it deterministically goes to subsequent states which generate the spike waveform, and finally returns to the initial rest state.

What does this imply for assumed inter-spike intervals? Spikes can’t overlap by construction, and the transition time between the rest and spiking states follows an exponential distribution. Well, technically it’s a geometric distribution but for simplicity let’s assume time is sampled sufficiently high that this is a good approximation. That means the implied interspike interval is a translated exponential distribution that has its peak at a delay given by the length of a spike waveform. So there’s an absolute but no relative refractory period here, which seems unrealistic.

So how can you recover a more realistic interspike interval? Calabrese & Paninski (2010) use a very simple solution involving adding a dummy refractory state. In fact this trick might be mentioned in Sahani’s thesis, and I’m sure it’s been noted countless times in texts on HMMs, but it didn’t seem obvious to me that this could work on first read so here’s the explanation.

If you add a dummy refractory state between the spike end and rest states, then the interspike interval will no longer be given by a variable following an exponential distribution but by the sum of two variables each following exponential distributions. The distribution of the sum of two exponential distribution is a hypoexponential distro, which looks like a gamma distribution. In fact if the means of the two distros are equal the sum of the two follows an Erlang distribution with k = 2, which in turn is a special type of gamma distribution.

So if you add up exponential random variables, with the first having lambda = 4 and the second lambda = 20, then you get a distribution which looks like this:

So by adding a dummy waiting state in your HMM, you can recover a more realistic distribution for waiting times. Applied to a spike-sorting HMM, this gives you the desirable relative refractory period.