3.1 Neurons

Author: Michael Beyeler; Stanislav Listopad

CARLsim currently supports Izhikevich spiking neurons with either current-based or conductance-based synapses (see 3.2 Synapses), but more neuron types are planned for the future. Different groups of neurons can be created from a one-dimensional array to a three-dimensional grid (see 3.3.2 Topography) via createSpikeGeneratorGroup and createGroup.

Differential equations can be integrated either with the forward Euler or Runge-Kutta method (see 12.2 Numerical Integration Methods).

3.1.1 Izhikevich Neurons (4-Parameter Model)

The Izhikevich neuron is a dynamical systems model that can be described by a two-dimensional system of ordinary differential equations:

$\begin{align*} \frac{dv}{dt} = & ~ 0.04v^2 + 5v + 140 - u + I & \text{(1)}\\ \frac{du}{dt} = & ~ a (bv - u) & \text{(2)} \end{align*}$

Here, (1) describes the membrane potential $v$ for a given current $I$ , where $I$ is the sum of all synaptic and external currents; that is, $I = i_{syn} + I_{ext}$ (see 3.2 Synapses). (2) describes a recovery variable $u$ ; the parameter $a$ is the rate constant of the recovery variable, and the parameter $b$ describes the sensitivity of the recovery variable to the subthreshold fluctuations of the membrane potential. All parameters in (1) and (2) are dimensionless; however, the right-hand side of (1) is in a form such that the membrane potential $v$ has mV scale and the time $t$ has ms scale. $a, b, c, d$ are open parameters that have different values for different neuron types.

In contrast to other simple models such as the leaky integrate-and-fire (LIF) neuron, the Izhikevich neuron is able to generate the upstroke of the spike itself. Thus the voltage reset occurs not at the threshold, but at the peak ( $v_{cutoff}=+30$ ) of the spike. The action potential downstroke is modeled using an instantaneous reset of the membrane potential whenever $v$ reaches the spike cutoff, plus a stepping of the recovery variable:

$\begin{align*} v(v>30) = & ~ c & \text{(3)}\\ u(v>30) = & ~ u + d & \text{(4)}\end{align*}$

The inclusion of $u$ in the model allows for the simulation of typical spike patterns observed in biological neurons. The four parameters $a, b, c, d$ can be set to simulate different types of neurons. For example, regular spiking (RS) neurons (class 1 excitable) have $a=0.02, b=0.2, c=-65, d=8$ . Fast spiking (FS) neurons (class 2 excitable) have $a=0.1, b=0.2, c=-65, d=2)$ . For more information on different neuron types, see Izhikevich (2003, 2004).

See also: 3.1.2 Izhikevich Neurons (9-Parameter Model)

3.1.2 Izhikevich Neurons (9-Parameter Model)

The 9-parameter izhikevich neuron model is equivalent to the 4-parameter one. Larger number of parameters is used for convenience in describing various neuron types.

$\begin{align*} \frac{dv}{dt} = & ~ (k (v - vr) (v - vt) - u + I) / C & \text{(5)}\\ \end{align*}$

Here, (5) describes the membrane potential $v$ for a given current $I$ , where $I$ is the sum of all synaptic and external currents; that is, $I = i_{syn} + I_{ext}$ (see 3.2 Synapses).

The parameter $k$ is a rate constant of the membrane potential, and can be found using the rheobase and the input resistance of a neuron. The parameter $vr$ is the resting membrane potential. The parameter $vt$ is the instantaneous threshold potential. The parameter $C$ is the membrane capacitance.

$\begin{align*} \frac{dv}{dt} = & ~ a (b (v - vr) - u) & \text{(6)}\\ \end{align*}$

Here, (15) describes a recovery variable $u$ ;

The parameter $a$ is the rate constant of the recovery variable. The parameter $b$ describes the sensitivity of the recovery variable to the subthreshold fluctuations of the membrane potential. The parameter $vr$ is the resting membrane potential.

All parameters in (5) and (6) are dimensionless; however, the right-hand side of (5) is in a form such that the membrane potential $v$ has mV scale and the time $t$ has ms scale. is voltage, $u$ is the recovery variable, $I$ is the input current, and $C, k, vr, vt, a, b, vpeak, c, d$ are open parameters that have different values for different neuron types.

In contrast to other simple models such as the leaky integrate-and-fire (LIF) neuron, the Izhikevich neuron is able to generate the upstroke of the spike itself. Thus the voltage reset occurs not at the threshold, but at the peak ( $v_{cutoff}=+vpeak$ ) of the spike. The action potential downstroke is modeled using an instantaneous reset of the membrane potential whenever $v$ reaches the spike cutoff, plus a stepping of the recovery variable:

$\begin{align*} v(v>vpeak) = & ~ c & \text{(7)}\\ u(v>vpeak) = & ~ u + d & \text{(8)}\end{align*}$

Regular spiking (RS) neurons (class 1 excitable) have $C=100, k=0.7, vr=-60, vt=-40, a=0.03, b=-2.0, vpeak=35, c=-50, d=100$ . Fast spiking (FS) neurons (class 2 excitable) have $C=20, k=1, vr=-55, vt=-40, a=0.15, b=8, vpeak=25, c=-55, d=200)$ . For more information on different neuron types, see Izhikevich (2003, 2004). For more information about 9-parameter model see Izhikevich (2007).

We do not recommend to use forward-Euler when working with the 9-parameter Izhikevich model. Use Runge-Kutta instead (see 12.2.2 Runge-Kutta Method).

See also: 3.1.1 Izhikevich Neurons (4-Parameter Model)

Since: v3.1

3.1.3 Multi-Compartment Neurons

Author: Stanislav Listopad; Michael Beyeler

Neurons in a group can be extended to include multiple compartments by calling CARLsim::setCompartmentParameters on a group and connecting them via CARLsim::connectCompartments.

In this case, the total current $I$ of (1) and (5) is extended to include a dendritic current, $i_{dendr}$ in pA; that is, $I = i_{syn} + i_{dendr} + I_{ext}$ . The dendritic current at each compartment consists of the currents coming from the down ("mother") compartment and up ("daughter") compartments:

$\begin{align*} i_{dendr} = & G_{down} (v - v_{down}) + \sum \limits_{up} G_{up} (v - v_{up}) & \text{(9)} \end{align*}$

with the values of the conductances $G_{up}$ and $G_{down}$ specified in CARLsim::setCompartmentParameters.

For more information on multi-compartmental neurons, see Izhikevich and Edelman (2008), especially their supporting information.

We do not recommend to use forward-Euler when working with multi-compartment neurons. Use Runge-Kutta instead (see 12.2.2 Runge-Kutta Method).

Note: The maximum number of allowed compartmental connections per neuron is controlled by the parameter MAX_NUM_COMP_CONN in carlsim_definitions.h.

Since: v3.1

3.1.4 Leaky Integrate-and-Fire (LIF) neurons

Author: Hirak J Kashyap

The support for LIF neurons has been added in CARLSim 5.0. LIF model and its variants are probably the most well known model of spiking point neurons, largely due their to simple membrane voltage dynamics favored for dynamical system analysis and hardware implementation.

The parameter $k$ is a rate constant of the membrane potential, and can be found using the rheobase and the input resistance of a neuron. The parameter $vr$ is the resting membrane potential. The parameter $R$ is the membrane resistance.

$\begin{align*} k \frac{dv}{dt} = & ~ - [v - vr] + R I& \text{(10)}\\ \end{align*}$

3.2 Synapses

Author: Kristofor D. Carlson

Spiking neural networks deliver information using spikes via structures called synapses. In biology, synapses pass information from one neuron to another via a chemical signal (electrical synapses also exist but are not modeled in CARLsim). Communication between synapses is, in general, unidirectional. Synapses are broken into two components: one component located on the neuron sending the information (the presynaptic neuron) and one component located on the neuron receiving the information (the postsynaptic neuron). The two synapse components are separated by a physical gap between the pre and postsynaptic neurons called the synaptic cleft. Information is passed from the presynaptic neuron to the postsynaptic neuron when the presynaptic component of the synapse releases neurotransmitters that cross the synaptic cleft and bind to postsynaptic receptors that result in a current influx into the postsynaptic neuron. The sum of many synaptic current contributions change the postsynaptic neuron voltage, causing it to spike if the voltage threshold is crossed.

CARLsim supports two synapse model descriptions. The current-based (CUBA) description uses a single synaptic current term while the conductance-based (COBA) description calculates the synaptic current using more complex conductance equations for each synaptic receptor-type. Both CUBA and COBA current contributions are influenced by the synaptic weight of the synapse. We next discuss each synapse model specifically.

(see Chapter 4: Connections).

3.2.1 CUBA

When current based (CUBA) mode is used, no conductances are present. The strength of the resulting current is directly proportional to the synaptic weight. Additionally, the total synaptic current at postsynaptic neuron $j$ , $I_{j}^{syn}$ , due to a spike from presynaptic neuron $i$ is given at any point in time by (5):

$\begin{align*} I_{j}^{syn} = & ~ \sum \limits_{i=1}^{N} s_{ij}w_{ij} & \text{(5)} \end{align*}$

where $s_{ij}$ is 1 if the neuron is spiking and 0 otherwise, $w_{ij}$ is the strength of the synaptic weight between postsynaptic neuron $j$ and presynaptic neuron $i$ , and N is the total amount of presynaptic connections onto postsynaptic neuron $j$ .

CUBA mode is the default behavior of all neurons. CARLsim does not currently support mixed modes consisting of some neurons with conductances and others without conductances in the same simulation.

Note: Because CUBA mode has synaptic currents that last for a single time step, and COBA mode has synaptic currents that decay over time, switching from CUBA mode to COBA mode usually requires the user to reduce the synaptic weights as the resulting currents will be much larger in COBA mode.

3.2.2 COBA

When conductance-based (COBA) mode is used, there are conductances with exponential decays present. If the synaptic connection is excitatory then AMPA and NMDA decay values are used. If the synaptic connection is inhibitory then GABA A and GABA B decay values are used. Equation (6) shows the corresponding current, $I_{ij}$ in postsynaptic neuron $j$ , due to inputs from N presynaptic neurons (denoted by $i$ ).

$\begin{align*} I_{j}^{syn} = & ~ \sum \limits_{i=1}^{N} i_{ij}w_{ij} & \text{(6)} \end{align*}$

Here we drop the $i,j$ notation and assume the $g$ values pertain to a specific synapse ( $i,j$ ) and the $v$ values pertain to a specific postsynaptic neuron, $j$ . The total current from a particular synapse can either be excitatory ( $i_{e}$ ) or inhibitory ( $i_{i}$ ). They are described in equations (7) and (8), below:

$\begin{align*} i_{e} = & ~ i_{NMDA} + i_{AMPA} & \text{(7)} \\ i_{i} = & ~ i_{GABA_{A}} + i_{GABA_{B}} & \text{(8)} \end{align*}$

The excitatory current is comprised of an AMPA current ( $i_{AMPA}$ ) and a voltage dependent NMDA current ( $i_{NMDA}$ ). They are described in equations (9) and (10).

$\begin{align*} i_{AMPA} = & ~ g_{AMPA}(v-v^{rev}_{AMPA}) & \text{(9)} \\ i_{NMDA} = & ~ g_{NMDA}\frac{\Big[\frac{v+80}{60} \Big]^{2}}{1+\Big[\frac{v+80}{60} \Big]^{2}} (v-v^{rev}_{NMDA}) & \text{(10)} \end{align*}$

Here, $g$ is the conductance, $v$ is the postsynaptic voltage, and $v^{rev}$ is the reversal potential. Both $g$ and $v^{rev}$ are specific to a particular ion channel/receptor. The inhibitory current is comprised of a GABA A current ( $i_{GABA_A}$ ) and a GABA B current ( $i_{GABA_B}$ ) shown in equations (11) and (12).

$\begin{align*} i_{GABA_{A}} = & ~ g_{GABA_{A}}(v-v^{rev}_{GABA_{A}}) & \text{(11)} \\ i_{GABA_{B}} = & ~ g_{GABA_{B}}(v-v^{rev}_{GABA_{B}}) & \text{(12)} \end{align*}$

The conductance, $g$ , for each ion channel/receptor, $k$ is given in the general form equation (13).

$\begin{align*} g_{k} = & ~ \sum \limits_{f} e^{t-t_{f}/ \tau} \Theta(t-t_f) & \text{(13)} \end{align*}$

Where $f$ indicates a particular spike event, $t_f$ represents the time of this spike event, and $\Theta$ is the Heaviside function.

The code snippet below shows how a user would first set a simulation that uses conductances using the CARLsim::setConductances function. If the CARLsim::setConductances is called with just the true argument:

sim.setConductances(true);

then, the default decay constant values will be used for the conductances (tdAMPA=5ms, tdNMDA=150ms, tdGABAa=6ms, tdGABAb=150ms) (instantaneous rise time).

The user can also specify the decay constant times themselves:

sim.setConductances(true,tdAMPA,tdNMDA,tdGABAa,tdGABAb);

It should be noted that users can specify a rise time constant for the slow the synaptic conductances (NMDA and GABAb) in addition to specifying a decay time constant. The rise time constants are set to zero by default. Below we show a user setting the rise times:

sim.setConductances(true,tdAMPA,tdNMDA,trNMDA,tdGABAa,tdGABAb,trGABAb);

Users can also explicitly set the default conductance time constant parameters to be used in the simulation with the command CARLsim::setDefaultConductanceTimeConstants.

sim.setDefaultConductanceTimeConstants(tdAMPA,tdNMDA,trNMDA,tdGABAa,tdGABAb,trGABAb);

Note: A pre-synaptic spike increases both post-synaptic AMPA and NMDA (GABAa and GABAb) for excitatory (inhibitory) connections.; Receptor-specific gain factors can be speficied to vary the AMPA-NMDA and GABAa-GABA ratio

see 4.1.4 Synaptic Receptor-Specific Gain Factors

3.3 Groups

Author: Michael Beyeler

3.3.1 Creating Groups

To create a group of Izhikevich neurons, simply specify a name (e.g., "output"), the number of neurons (e.g., 100), and a type:

int gOut = sim.createGroup("output", 100, EXCITATORY_NEURON);

Here, EXCITATORY_NEURON denotes that the neurons in the group are glutamatergic. Neurons with GABAergic synapses are supported with the INHIBITORY_NEURON keyword.

To refer to this group in later method calls, the CARLsim::createGroup method returns a group ID, gOut.

Next, specify the Izhikevich parameters, in this case for class 1 excitability (regular spiking) neurons (see 3.1.1 Izhikevich Neurons (4-Parameter Model)):

sim.setNeuronParameters(gOut, 0.02f, 0.2f, -65.0f, 8.0f);

where 0.02f, 0.2f, -65.0f, and 8.0f correspond respectively to the a, b, c, and d parameters of the Izhikevich neuron. Note we used gOut from the createGroup call above to reference the group ID.

To create a group of spike generators, the user also specifies a name, size, and type:

int gIn = sim.createSpikeGeneratorGroup("input", 10, EXCITATORY_NEURON);

3.3.2 Topography

Neurons in a group can be arranged into a (up to) three-dimensional (primitive cubic) grid using the Grid3D struct, and connections can be specified depending on the relative placement of neurons via CARLsim::connect. This allows for the creation of networks with complex spatial structure.

Each neuron in the group gets assigned a (x,y,z) location on a 3D grid centered around the origin, so that calling Grid3D(Nx,Ny,Nz) creates coordinates that fall in the range [-(Nx-1)/2, (Nx-1)/2], [-(Ny-1)/2, (Ny-1)/2], and [-(Nz-1)/2, (Nz-1)/2]. The resulting grid is a primitive cubic Bravais lattice with cubic side length 1 (arbitrary units). The primitive (or simple) cubic crystal system consists of one lattice point (neuron) on each corner of the cube. Each neuron at a lattice point is then shared equally between eight adjacent cubes. An example is shown in the figure below.

The 3D grid generated by Grid3D(2,2,2). Black circles indicate the 3D coordinates of all the neurons (with corresponding neuron IDs) in the group.

Examples:

Grid3D(1,1,1) will create a single neuron with location (0,0,0).
Grid3D(2,1,1) will create two neurons, where the first neuron (ID 0) has location (-0.5,0,0), and the second neuron (ID 1) has location (0.5,0,0).
Grid3D(1,1,2) will create two neurons, where the first neuron (ID 0) has location (0,0,-0.5), and the second neuron (ID 1) has location (0,0,0.5).
Grid3D(2,2,2) will create eight neurons, where the first neuron (ID 0) has location (-0.5,-0.5,-0.5), the second neuron has location (0.5,-0.5,-0.5), the third has (-0.5,0.5,-0.5), and so forth (see figure below).
Grid3D(3,3,3) will create 3x3x3=27 neurons, where the first neuron (ID 0) has location (-1,-1,-1), the second neuron has location (0,-1,-1), the third has (1,-1,-1), the fourth has (-1,0,-1), ..., and the last one has (1,1,1).
etc.

The 3D location of a neuron can be queried in all CARLsim states using the method CARLsim::getNeuronLocation3D. This allows, for example, for user-defined connections to be built based on 3D location (see 4.3 User-Defined Connections).

The following code snippet creates a group of 500 excitatory neurons arranged on a 10x10x5 three-dimensional (primitive cubic) grid:

int gOut = sim.createGroup("output", Grid3D(10,10,5), EXCITATORY_NEURON);

Note that in this example, the Grid3D parameter replaces the integer values (for number of neurons) used for CARLsim::createGroup in 3.3.1 Creating Groups. In fact, creating a group with N neurons is the same as arranging N neurons on a Nx1x1 grid:

int gOut = sim.createGroup("output", N, EXCITATORY_NEURON); // these two function calls

int gOut = sim.createGroup("output", Grid3D(N,1,1), EXCITATORY_NEURON); // are equivalent

See also: 4.1.3 The RadiusRF Struct; 4.1.8 Gaussian Connectivity

Since: v3.0

3.5 References

Izhikevich E.M. (2003) Simple Model of Spiking Neurons. IEEE Transactions on Neural Networks, 14:1569-1572.

Izhikevich E.M. (2004) Which Model to Use for Cortical Spiking Neurons? IEEE Transactions on Neural Networks, 15:1063-1070.

Izhikevich E.M. (2007) Dynamical Systems in Neuroscience - The Geometry of Excitability and Bursting.

Table of Contents