Today, Hackaday.io, you will be my rubber ducky. This log will cover implementing a crude approximation of Izhikevich neuron dynamics on a really cheap Atmel processor. In case you missed it, I discussed his model in the previous log, along with all the silly ways I corrupted the math in an attempt to make it run faster on an 8-bit, 8 MHz microcontroller.

First, a quick pseudo-code overview of the firmware:

int [a whole bunch of global variables];
main() {
    SystemInit;
    for(;;) { 
        waitforISR; 
        checkDendrites; 
        calcDendrites; 
        calcMembraneCurrent;
        calcIzhikevichModel; 
        if (firing) {fire;} 
        calcLED; 
        updateLED; 
    } 
}

The actual firmware is a bit more scattered since timing ends up being really important; in particular, the model calculations are broken into several pieces with LED updates in between to prevent flicker. If you want to see actual code... read on.

The complete v07 (again, more on that later!) firmware is here. I'll get it up on GitHub once I consider it "release-worthy".. I still need to debug a few timing issues with the 'scope. But feel free to check it out if you like. Below, I'll go over the non-commented portions, skipping the exhaustive intro that covers hardware and Izhikevich model outputs.

Global Variables

/*  Izhikevich model parameters and coefficients */
const int16_t a = 5;
const int16_t b = 2;
const int16_t c = -65;
const int16_t d = 8;
const int16_t E = 6;
const int16_t F = 2;
const int16_t G = 16;
const int16_t H = 161;

/*  Izhikevich model variables */
volatile int8_t I;
volatile int8_t I_rest = 21;
volatile int16_t v = -65;
volatile int16_t u = 0;

/*  Intra-model variables */
volatile int16_t vSq = 0;
volatile int16_t vF = 0;
volatile uint8_t modelStage = 0;
volatile int16_t v_prev;
volatile int16_t u_prev;

Here's where stuff gets a little bit more interesting. One of the more computationally expensive operations required in implementing the Izhikevich model involves squaring the current value of v. I ended up breaking this, along with a few other parts, out of the main calculation to save time. The modelStage variable keeps track of where the model currently is in this process.

/*  Neuron variables */
volatile uint8_t firing = 0;
volatile uint8_t fireTimer = 0;
const uint8_t axonPulseLength = 5;
volatile unsigned char dendStatus = 0; //current status of four dendrites, including types
volatile unsigned char dendStatusPrev = 0; //previous dendrite status/type
volatile int8_t val_Dend[4]; //keeps track of current dendrite's contribution to I (positive or negative)
volatile uint8_t stg_Dend[4] = {0,0,0,0}; //current stage of each dendrite impulse: 0 is REST, 1 is HOLD, 2 is DECAY
const int8_t hldTime_Dend[4] = {20,20,20,20};
volatile int8_t hldTimeCur_Dend[4];
const int8_t hldVal_Dend[4] = {10,10,10,10};
const int8_t dec_Dend[4] = {2,2,2,2};

/*  LED variables */
volatile uint8_t LED[3] = {0,0,0}; //R,G,B, 0-100
volatile uint8_t LEDtick = 0; //0-100
volatile uint8_t LEDfade[3][61] = {//fader for -102 <= v < -41, where I=21 (v=-83) is rest (pure green). R,G,B
{0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,1 ,2 ,2 ,3 ,3 ,4 ,4 ,5 ,5 ,6 ,6 ,7 ,7 ,8 ,8 ,9 ,9 ,10,10,11,11,12,12,13,13,14,14,15,15,16,16,17,17,18,18,19,19,20},
{1 ,2 ,3 ,4 ,5 ,6 ,7 ,8 ,9 ,10,11,12,13,14,15,16,17,18,19,20,20,19,19,18,18,17,17,16,16,15,15,14,14,13,13,12,12,11,11,10,10,9 ,9 ,8 ,8 ,7 ,7 ,6 ,6 ,5 ,5 ,4 ,4 ,3 ,3 ,2 ,2 ,1 ,1 ,0 ,0 },
{19,18,17,16,15,14,13,12,11,10,9 ,8 ,7 ,6 ,5 ,4 ,3 ,2 ,1 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 }
};

/*  Timing variables */
volatile uint8_t tick = 0; //ISR flips this to non-zero to update loop
const uint8_t modelUpdateFrequency = 255; //how many ticks pass before the model recalculates (higher = slower)
const uint8_t modelUpdateMultiplier = 3; //multiplier for modelUpdateFrequency (since they're 8-bit unsigned integers)

void getDendrites(void) {
/*  updates dendStatus to the current input value:
bit MSB 6   5   4   3   2   1   LSB
DEND    4   3   2   1   4   3   2   1
input   SIG SIG SIG SIG TYPE    TYPE    TYPE    TYPE
pin PB1 PB0 PB7 PD6 PB2 PD7 PB6 PD5
*/
    dendStatusPrev = dendStatus;
    dendStatus = 0;
    dendStatus |= ((PINB & (1<<PB1))<<6);
    dendStatus |= ((PINB & (1<<PB0))<<6);
    dendStatus |= ((PINB & (1<<PB7))>>2);
    dendStatus |= ((PIND & (1<<PD6))>>2);
    dendStatus |= ((PINB & (1<<PB2))<<1);
    dendStatus |= ((PIND & (1<<PD7))>>5);
    dendStatus |= ((PINB & (1<<PB6))>>5);
    dendStatus |= ((PIND & (1<<PD5))>>5);
}

getDendrites() works with a global 8-bit char called dendStatus, where the four least significant bits represent the four dendrite types (inhibitory or excitatory) and the four most significant bits represent the four dendrite statuses (high or low). I'll get into dendrite types more when I go over hardware--suffice it to say that each connector has four conductors (ground, power, type, signal) allowing upstream NeuroBytes and input modules to decide whether stimulus causes excitation or inhibition.

After saving the previous dendrite state into dendStatusPrev, the eight similar looking expressions use bitshift masks (such as 1<<PB7) to check various input pin states, then a further bitshift operation to store that state in the appropriate dendStatus bit. For example, the signal state of dendrite 3 is found in PB7 (i.e. Port B's MSB) but needs to be recorded in bit 5 of dendStatus; thus, the masked port read is shifted right two bits before getting incorporated into dendStatus.

void calcDend(uint8_t dend) {
/*  checks the dendrites and updates I based on hold and decay variables. function only runs on the
    dendrite indicated in the argument so it can be staged via updateModel(). */
    uint8_t inh = 1; //toggles to 0 if the selected dendrite is excitatory
    if (dendStatus & (1<<dend)) {
        inh = 0;
    }
    switch (stg_Dend[dend]) {
        case 0: //REST
            val_Dend[dend] = 0;
            if (dendStatus & (1<<(dend + 4))) {
                stg_Dend[dend]++;
            }  
            break;
        case 1: //HOLD
            val_Dend[dend] = hldVal_Dend[dend];
            if (dendStatus & (1<<(dend + 4))) {
                hldTimeCur_Dend[dend] = hldTime_Dend[dend];
            }
            else {
                hldTimeCur_Dend[dend] -= 1;
                if (hldTimeCur_Dend[dend] == 0) {
                    stg_Dend[dend]++;
                }
            }
            break;
        case 2: //DECAY
            if (dendStatus & (1<<(dend + 4))) {
                stg_Dend[dend]--;
            }
            else {
                if ((val_Dend[dend] - dec_Dend[dend]) > 0) {
                    val_Dend[dend] -= dec_Dend[dend];
                }
                else {
                    stg_Dend[dend] = 0;
                }
            }
            break;
    }
    if (inh == 1) {
        val_Dend[dend] = -val_Dend[dend];
    }
}

This function first checks to see if the currently selected dendrite is inhibitory or excitatory; if the latter, it sets the local variable inh to 0. After that, the function checks the stg_Dend array to see whether that dendrite is currently zero (i.e. contributing nothing to I), holding, or decaying. If the dendrite is stimulated, the stage is reset to HOLD and the appropriate timers are reset. Finally, the resultant value stored in the val_Dend array is negated if the input is inhibitory.

void calcI(void) {
/*  updates I based on current dendrite values */
    uint8_t i;
    I = I_rest;
    for (i=0;i<4;i++) {
        I += val_Dend[i];  
    }
}

void updateModel(stage) {
/*  This function updates the membrane potential (v) and recovery potential (u) variables based on the
    current I value. Since the Izhikevich model includes a reset check to determine firing, this function
    also updates the 'firing' variable to 1 when the neuron fires. The model uses 16-bit integer math, so
    interrupts are temporarily disabled to maintain atomicity.
    Execution time: 29.0 us max, depending on current stage. */
    uint8_t temp = SREG;
    cli();
    if (stage == 0) {
        getDendrites();
    }
    else if ((stage >= 1) & (stage < 5)) {
        calcDend(stage - 1);
    }  
    else if (stage == 5) {
        calcI();
    }
    else if (stage == 6) { //17.4 us
        vSq = square(v);
    }
    else if (stage == 7) {
        v_prev = v;
        u_prev = u;
        vF = v * F;
    }
    else if (stage == 8) { //29.0 us    
        if(v_prev > H) {
            v = c;
            u = u_prev + d;
            firing = 1;
        }
        else {
            v = v_prev + (vSq >> E) + vF + G - u_prev + I;
            u = u_prev + (((v_prev >> b) - u_prev) >> a);
            firing = 0;
        }
    }
    else if (stage == 9) { //28.9 us
        translateColor();
    }
    SREG = temp;
}

Stage 0: update the current dendrite status using getDendrites().

Stage 1, 2, 3, 4: calculate the four dendrite contributions to I using calcDendrites(dend).

Stage 5: sums the four dendrite contributions to update I using calcI().

Stage 6: starts the Izhikevich calculations by squaring the membrane potential, v, and storing the result in vSq. The square(int) function also checks for 16-bit signed integer overflow, and if present, instead returns 32767.

Stage 7: continues the Izhikevich calculations by multiplying v by F, and storing the result in vF. Also stores the current values of v and u since both of the next calculations need to run with each variable.

Stage 8: finish the Izhikevich calculations by checking for a fire state (v > H) and either resetting the model or finishing the calculation of v and u. Note, again, that E, a, and b have been replaced with right bit-shifts (i.e. 1/n^2).

Stage 9: translate the current membrane potential value into LED values (shown below).

Finally, the interrupt register is restored to its previous value. Note that when Stage > 9, this function doesn't do anything.

void translateColor() {
/*  translates membrane potential (v) values into the RGB array
    resting membrane potential:
    v < -102        pure blue
    -102 <= v < -90     fade blue to green
    v = -90         pure green
    -90 < v =< -71      fade green to red
    v > -71         firing (pure white)
    Execution time: 28.9 us    
*/
    if(v < -102) { //3.5 us
        LED[0] = 0;
        LED[1] = 1;
        LED[2] = 20;
    }
    else if((v >= -102) & (v < -41)) { //24.1 us
        LED[0] = LEDfade[0][(uint8_t)(v + 102)];
        LED[1] = LEDfade[1][(uint8_t)(v + 102)];
        LED[2] = LEDfade[2][(uint8_t)(v + 102)];
    }
    else {
        LED[0] = 20;
        LED[1] = 1;
        LED[2] = 0;
    }
}

void updateLED(void) {
/*  Updates the RGB LED based on the current membrane potential value. Uses PWM fading and global variables to
    figure out when the various elements should be on. Resolution controlled by LEDtick reset.
    Execution time: 5.6 us */
    if(LED[0] > LEDtick) {
        PORTC &= ~(1<<PC3);
    }
    else {
        PORTC |= (1<<PC3);
    }
    if(LED[1] > LEDtick) {
        PORTC &= ~(1<<PC4);
    }
    else {
        PORTC |= (1<<PC4);
    }
    if(LED[2] > LEDtick) {
        PORTC &= ~(1<<PC2);
    }
    else {
        PORTC |= (1<<PC2);
    }
    LEDtick++;   
    if (LEDtick > 100) {
        LEDtick = 0;
    }  
}

ISR(TIMER0_COMPA_vect) {
/* Makes ticks fire at a regular interval in the main loop. */
    tick = 1;
}

void systemInit(void) {
/* set up D1 */
    DDRC |= ((1<<PC2) | (1<<PC3) | (1<<PC4));
    PORTC |= ((1<<PC2) | (1<<PC3) | (1<<PC4)); //set pins = LED off
 
/* set up dendrites */
    DDRB &= ~((1<<PB7) | (1<<PB0) | (1<<PB1) | (1<<PB6) | (1<<PB2));
    DDRD &= ~((1<<PD5) | (1<<PD6) | (1<<PD7));
/* set up axons */
    DDRD |= (1<<PD0);
    DDRC |= ((1<<PC1) | (1<<PC0) | (1<<PC5));
    PORTC |= ((1<<PC0) | (1<<PC5)); //set type pins
 
/* set up Timer/Counter0 */
    TCCR0A |= ((1<<CTC0) | (1<<CS00) | (1<<CS01)); //CTC, clk/64
    TCNT0 = 0;
    OCR0A = 4; //loop time = 8 * (OCR0A + 1) uS
    TIMSK0 |= (1<<OCIE0A); //enables Output Compare Match A ISR
 
/* misc */
    sei(); //enable global interrupts
}

void fire(void) {
    PORTC &= ~((1<<PC2) | (1<<PC3) | (1<<PC4));
    fireTimer = axonPulseLength;
}

Main Loop

int main(void) {
    systemInit();
    uint8_t i = 0; //counters for update delay
    uint8_t j = 0;
    for(;;) {
        while(tick == 0) { /* idle loop */ }
        /* This stuff happens every 40 uS or so */
        tick = 0;
        if (i < 10) {
            updateModel(i);
        }
        if (firing == 0) {
            updateLED(); //5.6 us
        }
        else {
            fire();
        }
        j++;
        if (j == modelUpdateMultiplier) {
            i++;
            j = 0;
        }
        if (i == modelUpdateFrequency) {
            i = 0;
            if (fireTimer > 0) {
                PORTD |= (1<<PD0);
                PORTC |= (1<<PC1);
                fireTimer--;
            }
            else {
                PORTD &= ~(1<<PD0);
                PORTC &= ~(1<<PC1);
            }
        }
    }
}

This loop ties everything together. The main loop uses modelUpdateMultiplier and modelUpdateFrequency in conjunction with i and j, two local counter variables, to run updateModel() and iterate fireTimer when needed. Note the use of a while loop to take up all the tick slack; during testing, I can easily comment this line out, add in a debug pulse on an extra I/O line, and use my oscilloscope to quickly figure out the execution time for each stage.

And here it is running on a pair of teeny v07 boards:

I'll get back to hardware next. Comments on the code are welcome!