I've been slacking on updates--turns out when you get to work on a project full time, it doesn't make logs magically appear. In any case, this log post will focus on one thing: improving the core NeuroBytes firmware. I'll bring ya'll up to speed on other stuff in the next few weeks.

If you ever took a look at the v04 code, you may have noticed that the neuron membrane potential decay algorithm was quite simplistic:

if (potentialTimerCounter >= potentialTimerOverflow) {
	decayPotential = (decayPotential * 95) / 100;
}

This pretty much takes the current value and multiplies it by 0.95. There were a few other snippets here and there that did stuff too--if the current potential was over 100, for example, the code called that a "fire", immediately dropped the membrane potential by _a lot_, and made the LED blink white. The resultant boards looked reasonably correct, but the model didn't respond quite as expected when stimulated at higher levels. Neurons also fired _immediately_ rather than ramping up over time. Other issues included an unacceptably low refresh rate, occasional atomicity problems, and general inflexibility--it wasn't possible to change individual dendrite weighting, for example. Time for a rewrite.

---

In 2003, Eugene Izhikevich published a paper describing a computationally efficient method for simulating a variety of biological neuron types in software. His results were impressive: using a consumer-grade 1 GHz desktop PC, he was able to simulate a 10,000-element spiking cortical neuron network in real time with 1 ms clock resolution. A few years ago, @Bruce Land successfully implemented Izhikevich's model on an FPGA in his #Spiking neural net simulated on FPGA project.

Electronic version of the figure and reproduction permissions are freely available at www.izhikevich.com

You should probably stop reading this now and go read Izhikevich's paper, and then check out Bruce's project. In a nutshell, the model uses a pair of differential equations to keep track of membrane potential (v) and membrane recovery (u), along with a limit check that resets the model when the neuron fires. a, b, c, and d are parameters related to the specific type of neuron being modeled, while I is the cell's constant current input. As far as I can tell, these parameters were determined empirically based on observations of actual neurons, and the resultant model is quite accurate.

The first step was to recreate Izhikevich's model in a spreadsheet so I could wrap my head around the equations. This is what I eventually came up with, graphed over 0 < t < 500 ms:

The sheet is actually pretty fun to putz around with--you can change the parameters as Izhikevich suggests to emulate different types of neurons, such as fast spiking:

I also broke out I, as well as the original coefficients Izhikevich used in his differential equations; I called them E (the square coefficient), F (the linear coefficient), G (the offset coefficient), and H (the reset voltage threshold). They were also fun to putz with, although pretty much every minor change resulted in something like this:

Disastrous results from a tiny change, as you can see. Maybe I shouldn't mess with the coefficients.

---

So what's the issue? Why am I messing with the coefficients? Why not just copy the equations directly from LibreOffice Calc into VIM and call it a day?

In a word, speed. Remember, I'm trying to reduce flicker while PWMing three LED elements using an 8 MHz 8-bit microcontroller. Based on the LEDs I'm using and the desired dimming range (and a lot of waving test boards around), I need to update the LEDs every 30 us or so--that's 240 clock cycles. Since the ATtiny doesn't have anything resembling a floating point unit, everything also needs to be worked out in integers.

I tried adding int() tags in various spots in the equations, but I didn't get very far. The first coefficient divides the square of v by 25, and rounding the quotient inevitably caused huge errors. Once I started introducing these errors--even when I managed to make them pretty small--I began to realize just how delicate the balance between a, b, c, d, E, F, G, H, and I really was. Time for another computer program! Here's a test version that uses floating point math to prove out the concept:

... results:

Yeah, I totally was going to try (E = 0.031384, F = 4.654088, G = 153.545868, H = 75.376656) next!

The program is actually pretty simple. It selects four random coefficients within a predefined range, runs through the entire 0 < t < 500 sheet using Izhikevich's original equation format (and standard spiking neuron parameters), checks the "fit" of the result, and stores the best "fit" coefficients until it's run through an unreasonable number of iterations. It's a brute force approach that works well enough for floating point math, as shown above.

To produce AVR-ready code, I forced the program to use integers in its calculations. I decided to change Izhikevich's original equation a bit to hopefully speed up the calculation--namely, I changed E, a, and b to bitshift operations, meaning they're restricted to various n values in 1/(2^n). I also set up boundary conditions that prevent 16-bit signed integer overflows. Here's what I (eventually) came up with:

It's far from perfect, but the important dynamics work properly--namely, changing the value of I causes the fire rate to increase or drop to zero, and there is a distinct (although not as clear as in the original model) inflection point before each spike. You can see a few "known good" coefficient/parameter sets I came up with--getting to these required fiddling with fitment weights quite a bit, so I wrote 'em down whenever they worked decently well. If you're curious/bored/etc, the complete program is here.

You'll also notice that I left a bunch of new parameters on the screenshot above, along with a second plotted line--those are part of the fitment parameters. The program selects a series of points on the membrane voltage plot where it thinks the neuron experiences an action potential; this is based on the rate of change of this value (the "slope threshold"), and defines the linear and firing regions. The "firing?" line helped me quickly select good slope threshold values. Once this point is selected, the program calculates four parameters, averaged across 0 < t < 500:

linear region slope avg (0.252)

linear region range pct (0.180)

fire ratio (0.132)

number of firings (6)

Combined, these four parameters roughly describe what a neuron membrane potential waveform should look like; the ideal values from Izhikevich's regular spiking neuron are shown above in parentheses.

Up next: implementation in AVR-C. Which is done and works rather well, but I have already written far too much for one log update.