<< Main Manual

Tutorial #8. Implementing Integrate-and-Fire Models

This tutorial assumes that you have completed Tutorial #7.

MacGregor model

To implement integrate-and-fire models with discontinuous regions, we can embed C++ code in the .oden files. The following code is modified from the FORTRAN code for PTNRN10 in Neural and Brain Modeling by Ronald MacGregor. The nne_build parser looks for a function starting with "inline" at the end of the .oden files. The variable type "fp_type", is a generic floating point number (usually implemented as a double precision floating point number).

MacGregor.oden:

############################################################################
# Copyright 2025 John A. Hayes (hayesjohn@ufl.edu, jahaye1@gmail.com)
# Model of PTNRN10 from MacGregor 1987
# https://www.google.com/books/edition/Neural_and_Brain_Modeling/Nih-AAAAIAAJ
############################################################################

# Cellular Variables
p C = 0.5
p TTH=30.0
p B=5.0
p E_K = -10.0
p TGK=3.0
p TH0=10.0
p T_MEM=5.0

# Pulse Variables
p SC=20.0
p SL=0.0
p SCSTRT=1.0
p SCSTP=100.0
p LTSTOP=100.0

# Initial Conditions
init V_m = -50.0
init G_K = 0.0
init TH = TH0

# Auxiliary Equations
SCN(t) = stepPulse(t)
S(t) = calculateS(V_m, TH)

# Differential Equations
# V_m is MacGregor's "E"
dV_m/dt = V_m*exp(-(1.0+G_K)*dt/T_MEM) + (SCN(t) + G_K*E_K)*(1.0-exp(-(1.0+G_K)*dt/T_MEM))/(1.0+G_K) - V_m
dG_K/dt = G_K*(exp(-dt/TGK))+B*S(t)*(1.0-exp(-dt/TGK)) - G_K
dTH/dt = TH0+(TH-TH0)*exp(-dt/TTH) + C*V_m*(1.0-exp(-dt/TTH)) - TH

# Auxiliary Functions
inline fp_type stepPulse(fp_type time) const {
   if (time >= SCSTRT && time < SCSTP) {
      return (SC + SL*(time-SCSTRT));
   }
   return 0.0;
}

inline fp_type calculateS(float V_m, float TH) const {
   //cerr << "S: " << S << endl;
   double S_tmp = 0.0;
   if (V_m >= TH) {
      S_tmp = 1.0;
   }
   else {
      S_tmp = 0.0;
   }
   return S_tmp;
}

Note that this model uses the first-order Euler integration, so we have to change the .setup file to specify 'euler' integration:

Network MacGregor.net
Dt 1
Print Step Size 1
Output V_m G_K TH S SCN
Method: euler
Setup Seed 1
Simulation Seed 1
Alias none
Experiment control 2

This produces something like the following:
MacGregor.net

If we need to modify multiple state variables in our C++ functions we can use the "fp_type *dVars" argument to access the array of dependent variables where the name of the index for state variable V_m will be "Breen_V_m_INDEX". The following is an integrate-and-fire model of a Na⁺-dependent "pacemaker" neuron published by Breen et al.:

Breen.oden:

############################################################################
# Copyright 2025 John A. Hayes (hayesjohn@ufl.edu, jahaye1@gmail.com)
# Based on the model from https://doi.org/10.1162/089976603322518768
############################################################################
# Single compartment Breen (2003) integrate-and-fire bursting model
#

# Parameter Values
p C_m=21.0,

p E_Na=50.0
p E_L=-65.0
p E_syn=0.0

p theta_mNaP=-40.0
p theta_hNaP=-48.0
p theta_s=-44.0

p sigma_mNaP = -6.0
p sigma_hNaP=6.0
p sigma_s=-2.0

p tau_hNaP_max=10000
p tau_s=5

p g_NaP=2.8,
p g_L=2.8
p g_syn=0.025

p I_stim=20.0

# Differential equations for state variables
V_m'=(-I_NaP-I_L-I_syn+I_stim)/C_m
h'=(hNaP_ss-h)/(tau_hNaP)
s_AMPA'=((1-s_AMPA)*(s_inf(V_m))-s_AMPA)/tau_s

# Initial conditions for state variables
init V_m=-50.0
init h=0.48
init s_AMPA=0.0

# Current equations
I_NaP=g_NaP*mNaP_ss*h*(V_m-E_Na)
I_L=g_L*(V_m-E_L)
I_syn=g_syn*sum(s_AMPA)*V_m

# Auxiliary Equations
S(t) = calculateS(dVars, V_m, h)

mNaP_ss=1.0/(1.0+exp((V_m-theta_mNaP)/sigma_mNaP))
hNaP_ss=1.0/(1.0+exp((V_m-theta_hNaP)/sigma_hNaP))
tau_hNaP=tau_hNaP_max/(cosh(((V_m-theta_hNaP)/sigma_hNaP)/2.0))
s_inf(V_m)=1.0/(1.0+exp((V_m-theta_s)/sigma_s))

inline fp_type calculateS(fp_type *dVars, fp_type V_m, fp_type h) const {
   double V_TH = -18.8970 - 200.41*h + 414.62*h*h - 207.88*h*h*h;
   //cerr << "V_m: " << V_m << " h: " << h << " V_TH: " << V_TH << endl;
   if (V_m >= V_TH) {
      //cerr << "resetting!" << endl;
      dVars[Breen_V_m_INDEX] = -47.359 - 6.5825*h + 11.085*h*h;
      dVars[Breen_h_INDEX] = h -0.00078 - 0.0037*h;
      return 1.0;
   }
   return 0.0;
}