The Gnuron Manual

Introduction

Gnuron is a first attempt at an extensible, general purpose neuron simulator. The program is extensible in that new neuron models, current inputs, statistics and the like can be added without changing the structure of the program. It is general purpose in that the models, currents and parameters are defined by the user in a text file which is parsed at the start of each run.

Gnuron was written by John D. Hunter; jdhunter@ace.bsd.uchicago.edu.

Running the Program

The program is executed from the command line with 

gnuron [options]

gnuron will parse from the standard input; otherwise specify an input file containing the instructions with the `-p' option. See section Input File, for a detailed description of the input format. The instructions from the input file will be used to generate a population of model neurons which are stimulated with the currents defined in the input over a range of parameters. The program can be used to compute statistics on the neural population, or simply to generate spike trains, voltage traces and current records for analysis and viewing in another program.

The following command line options are recognized.

-c
Save the current traces
-d dir
Place output in dir; See section Output Files.
-p param-file
Parse input form param-file
-q
Quite mode: no standard output
-s
Save the spike trains
-v
Save the voltage traces

Input File

The input file is used to declare everything you need for a run. This includes the neuron model, the integration parameters, the input currents, the parameter space and the statistics to be computed.

Basic Syntax

All statements must end with a `;'. Both C and C++ style comments are allowed. `//' will comment out anything to the end of the line, and anything enclosed between `/*' and `*/' will be ignored. Caution: the second form of comment cannot be nested.

Here is a brief summary of symbols:

=, +, -, *, /
Their usual meanings.
:=
Function definition; See section User-defined Functions.
->
Used to set properties of current objects; See section Current Inputs.
neuron:
Declare the neuron model; See section Neuron Models.
terminate:
Specify a condition to terminate integration; See section Terminating Integration.
compose:
Specify how to compose multiple terminate: conditions; See section Composing Terminators.
repeat:
Specify how many times to repeat each input; See section Repeating Inputs.
current:
Define a new current input to the model neurons; See section Current Inputs.
stat:
Specify a statistic to compute on the neuron population; See section Statistics.

Neuron Models

The models currently provided are leaky integrate and fire, Morris Lecar, and a first and second order delay differential equation. To specify a model, use the syntax: 

neuron: model(number);

where model is either delay, delay2, leaky or morris or . number should be the number of neurons you want to use, i.e. the number of neurons in a network. Currently I only have uncoupled populations, but, coming soon to theaters everywhere: synapses! To declare a population of 100 leaky integrate and fire neurons use: 

neuron: leaky(100);

Variables are also ok here (and everywhere) so this is fine to: 

N = 100;
neuron: leaky(N);

Note there is a difference between declaring a population of neurons and using the repeat: command; See section Repeating Inputs.

The default variables for each model and their values are given in section leaky variables, section morris variables, section delay variables, section delay2 variables.

Functions and Variables

User-defined Variables

You define variables in the natural way: 

N = 2;
f = 20.0;
omega = 2*pi*f;

Note, you can also define a variable to be the left hand side of an assignment to a linspace, ranspace, etc... parameter generator. You can still use the variable naturally in future expressions, and in each such expression the values of the parameters generated by the parameter generators will be substitued in for the variable; See section Parameter Space.

User-defined Functions

You can only define functions of the special variable `t', which is a place holder for the current time in the integration of the neuron models; See section Built-in Variables. To define a function use `:=' rather than `='; this implcity declares the left hand side of the `:=' operator to be a functions of `t'. 

lambda = 20.0;
f := sin(pi*t)+exp(-lambda*t);

Call a user defined function in the same way you call all functions of one argument `func(arg)'. 

x = f(0.0) // x is now 1.0
g := sin(f(t)+3.0);

Functions can be arbitrarily composed combined as above.

Built-in Variables

The constants `pi' and `e' are provided.

`t' is a special variable which is used to defer evaluation of an expression until the integration. At that time, the current time value will be substituted in for `t' in all the expressions where it is used. When you define a function using the `:=' operator, you are implicity declaring a function of `t'. You can not define functions of any other variables; See section User-defined Functions.

Each neuron model defines a number of variables with default values which are loaded and available at the declaration neuron: (see section Neuron Models). Each model's variables are prefixed with the name model from the `neuron: model' declaration, e.g. leaky.rm for the membrane resistance of the leaky integrate and fire model. For the leaky model see section leaky variables, for the morris model see section morris variables, for the delay model see see section delay variables, and for the second order delay2 model see see section delay2 variables.

Built-in Functions

A growing number of deterministic and stochastic functions are provided.

Regular Functions

The built-in functions are called with 

func(arg1[, arg2, ...])

where func is one of the functions listed below and arg is any legal variable or function, either user-defined or built-in (with some exceptions listed below). The special variable `t' may also be used in any argument; in this case evaluation will be deferred until integration.

eta(arg)
Gaussian white noise w/ zero mean and unit variance
exp(arg)
exponential
cos(arg)
cosinusoid
delta_conv(point-process, response-function)
First argument must be a point process (periodic_point, poisson or Delta). Convolution of point process with response function.
delta_conv_efficient(point-process, response-function, ignore-distance)
As above but ignoring points that are farther away than ignore-distance.
Heavi(arg)
Heaviside step function; 0 for arg < 0 and 1 otherwise.
ln(arg)
natural logarithm
periodic_point(freq, dt)
A periodic point process with instantaneous frequency (which can be a function of time) given by freq. When the absolute value of the distance between the current time and the last spike is within dt of the instantaneous period, a spike will occur.
poisson(arg)
A poisson process with instantaneous probability given by arg. Multiplying the desired rate by the step-size will give the required instantanouse probability. arg must be between 0 and 1.
sin(arg)
sinusoid
sqrt(arg)
square root

Examples using some of the less straight-forward functions.

Heavi(arg) 
alpha = 100.0;
resp := 0.5*alpha*t*exp(-alpha*t)*Heavi(t); // impulse response
periodic_point(freq, dt) 
dt = 0.001;
f = 2.0;
train := periodic_point(2*f, dt);
poisson(arg) 
dt = 0.001;
train := poisson( 5.0*dt ); // constant 5 Hz Poisson
poisson(arg) 
train2 := poisson( (1.0+sin(pi*t))*dt ); // modulated Poisson
delta_conv(point-process, response-function) 
syn := delta_conv(train1(t), resp(t)); // synaptic current
delta_conv_efficient(point-process, response-function, ignore-distance) 
syn := delta_conv_efficient(train(t), resp(t), 10.0/alpha);

Constructed Functions

Some functions require initialization before they can be used. This is analogous to constructing a new object with initial parameters using a `new' operator in some object oriented languages. To construct a function, define a new function using the function definition operator `:=' (see section User-defined Functions). 

f := FunctionConstructor(init1[, init2, ...]).

with the initialzation arguments init1, init2 etc .... Then f is called in the usual way: `f(arg)'. Function constructors begin with uppercase letters.

Here are the functions which require construction before use:

Dirac(dt)
The Dirac delta-distribution; returns 1/dt for for |arg| < dt/2.
GaussianColored(lambda, dt)
Exponentially-correlated, Gaussian-distributed, colored noise with zero mean, unit variance and correlation rate lambda; the units of lambda are 1/s.

Example usage

Dirac(dt) 
dt  = 0.001;
delta := Dirac(dt); // delta(arg) = 1/dt for |arg| < dt/2
x = delta(0);        // x = 1000.0
x = delta(dt);       // x = 0.0
GaussianColored(lambda, dt) 
// inject colored and white noise
lambda = 1.0;
cnse := GaussianColored(lambda, dt);
current: s = s1*cnse(t) + s2*eta(t);

If `func' is a constructed function, you can use it in the definition of new functions in the usual way; note that it is illegal to use the assignment operator `=' to define new functions. 

cnse := GaussianColored(lambda, dt);
f := cnse(t) + eta(t);   // right 
f = cnse(t) + eta(t);    // wrong
current: s = f(t);

Current Inputs

You can declare as many currents as you like, each of which can have different properties. These include the reversal potential and the neural targets.

Declaring a Current

To declare a current input to the neuron(s), use this syntax: 

current: handle = function;

function can be an arbitrary function of time, composed of user defined and built-in functions and variables; See section Functions and Variables. handle is any string which can be used to refer to the current in the rest of the input file; See section Basic Syntax, for legal names. It is used to set the current properties, such as which neurons receive the current and whether it has a reversal potential. 

f = 2.0;
current: s = sin(2*pi*f*t);  // a 2 Hz sinusoidal input

With a handle to a current, such as `s' is in the example above, you can set current properties with the `->' operator. For example, to make the sinusoid a coherent input into the entire population (see section Setting the Targets) use: 

s->coherent;

This is more efficient that generating an identical sine for each neural element.

Setting the Targets

I have not yet included support for currents going to individual neurons but this will come soon. For now, any declared current will project to every neuron in the population, but a current can be declared either independent or coherent. If the current is independent, every neuron will receive a different instance. For stochastic currents, this means that each element will get a noise input which has the same statisitcal properties (same mean, variance, etc ...) but different time values; for deterministic inputs the independent and coherent currents will generate the same values, but it is more efficient to use the coherent designation because the value will only be computed once in this case, and will be recomputed for each element in the population for the independent case. The default behavior is independent. 

sigma = 0.1;
current: nse = sigma*eta(t);
nse->independent;              // every neuron's copy is different

In my research, I often want to stimulate the population with a coherent aperiodic current that is generated by a stochastic process. For example, if a neural population receives a common synaptic input which is poisson distributed, it would be declared like this: 

dt = 0.001;                                     // the stepsize
train := poisson(5.0*dt);                       // 5 Hz Poisson train
alpha = 200.0;
resp := 0.5*alpha*t*exp(-alpha*t)*Heavi(t);     // impulse response
current: epsps = delta_conv(train(t), resp(t)); // convolution 
epsps->coherent;     // every neuron gets an identical copy

Setting the Reversal

Certain current inputs, such as synaptic currents, have reversal potentials in real neurons. Once you have defined a current (see section Current Inputs), as in: 

current: epsps = delta_conv(train(t), resp(t));

you can set the reversal potential with: 

Vr = -0.060;               // the reversal potential in volts
epsps->set_reversal(Vr);

By default, currents do not reverse.

Terminating Integration

Use the terminate: flag to specify a condition for terminating integration with this syntax: 

terminate: criterion(arg1[, arg2]);

where criteria is either time or spikes. The time criterion specifies how long (in seconds) you should integrate for, which is given as a single argument 

terminate: time(5.0);  // integrate for 5 seconds

The spikes criterion specifies how many spikes the neuron must fire before integration is terminated, and is declared as 

terminate: spikes(spike-number[, neuron-number]);

If the optional second argument neuron-number is omitted, it will default to `0', the index of the first neuron in the population. [The indexing of neuron-number is C-style which starts at 0.]

Integration will terminate when the specified neuron has fired spike-number spikes. 

terminate: spikes(1);   // terminate on the 1st spike of neuron 0
terminate: spikes(5,3); // terminate on the 5th spike of neuron 3

Caution: when using the spikes criterion, it is strongly recommended to also use a time criterion to insure that the integration eventually terminates. If you inadvertantly use only spikes with a subthreshold current: input, the program will never terminate; See section Composing Terminators.

Composing Terminators

Since, multiple terminate: conditions may be specified, you need a way of composing them. This can be useful, for example, if you are investigating a parameter space and you are not sure one of you conditions will be satisfied for all parameter values and you want to insure that the program eventually terminates. So it it advisable to always put an upper time limit on how long you want to integrate. The syntax of the compose: declaration is: 

compose: rule

where rule is either and, or and xor. Thse will be satisified if all of the terminate: conditions are true (and), any of them are true or, or exactly one is true xor. It is easy to add more sophisticated composition rules. If multiple terminate: conditions are given, but no compose: condition is specified, the default rule is or. 

terminate: spikes(1);  
terminate: time(100.0);
compose: or;  // terminate with first spike or 100 seconds

Repeating Inputs

The repeat: declaration is used to repeat a give set of inputs to a neural population, and is useful if you are using stochastic inputs and want to measure some statistical property of the response. The syntax of the declaration is: 

repeat: number;

where number is the number of times you want to repeat the set of input currents; See section Current Inputs.

For example, if you want to compute the first passage time of a first order delay differential equation you declare the model as 

neuron: delay(1);  // a 'population' of 1 element
repeat: 500;       // repeat current inputs 500 times

This would repeat whatever inputs you specify 500 times to a single delay model.

Omitting the repeat: declaration is equivalent to `repeat: 1;'.

Parameter Space

It is very easy to sample a parameter space with the following syntax: 

var = grid(min, max, num-points);

This will generate num-points between min and max (inclusive) with a spacing determined by the grid function you specify. Currently two grid functions are provided: linspace and ranspace. The former is a grid with equal spacing between the points and the latter samples from a uniform distribution. Coming soon to theaters everywhere: logspace and gaussianspace. var is a user defined or built-in variable (see section Functions and Variables.

The following declaration will cause the program to loop 100 times, generating a new neural population each time; see section Neuron Models. Every variable, function and current containg alpha will be recomputed for all the points in the parameter space: 

N = 100;                                
min = 0.0;
max = 2.0;                           
alpha = linspace(min, max, N);   // linear spacing

Each point in the parameter space will be repeated as specified in the repeat: declaration; see section Repeating Inputs.

You can generate a multi-dimensional parameter space by using multiple parameter generating commands. The following will sample 10,000 points at all amplitude/frequency combinations: 

N = 100;                                
amp = linspace(0.0, 1.0, N);    // amplitude
freq = linspace(0.0, 2.0, N);   // freq

Statistics

Still in its infancy but promising to grow soon. Currently only one statistic method which is useful for first passage time problems. Call stat: first_spiketime; to compute the first passage time.

Output Files

By default, only the results of the requested statistics (see section Statistics) are output. The voltage, current and spike train records are obtained by passing the proper flags at the command line; See section Running the Program.

Standard Output

The standard output will look something like 

pnt     iter    dde.alpha   first_spike
1       1       0.9            -1
1       2       0.9            -1
1       3       0.9            -1
2       1       1              15.288
2       2       1              16.441
2       3       1              17.463
3       1       1.1             5.559
3       2       1.1             3.977
3       3       1.1             7.839

The first column gives the point in param space, the second column the iteration number, the next columns will indicate the values of the parameters specified in the param space (if any; see section Parameter Space) and the remaining columns will give the values of the statistics you specified (if any; see section Statistics).

The header of these columns will be saved to a file called header and the data in the columns will be stored in stats.dat.

Default Output Files

These files will be put in the default output dir; See section Running the Program.

`param.input'
A copy of the input param-file specified at the command line (if any).
`stats.dat'
The data displayed to the standard output, sans header; See section Standard Output.
`header'
The header of the data displayed to the standard output.

Voltage Current and Spikes

The voltage traces, current traces and spike trains are all available optionally by specifying appropriate command line options; See section Running the Program. These will be placed in respective sub-directories of the default output directory: `volts', `currents' and `spikes'. Files in these three sub-directories have the same naming convention. 

pnt#_iter#.dat

where the first # symbol is the identifier for the point in the parameter space begin computed (from the first column of `stats.dat'), and the second # symbol is the iteration number (the second column of `stats.dat'); e.g., the first point computed will use the output file name `pnt1_iter1.dat'.

Each of these files will have N lines, where N is the number of neural elements in the declaration: 

neuron: model(N);

See section Neuron Models. The first line of each file is the voltage/current/spike record of the first model neuron, the second line is the record of the second model, etc ....

Examples

Some sample input files are given here with a description of what they compute.

First Passage Time Density

This file computes the first passage time to threshold for a delay differential equation: 

// --- dvdt = alpha*v(t) + beta*v(t-tau) + input(t) --- //
neuron: delay(1);       // initialize DDE
repeat: 1000;           // number of repetitions of each point
terminate: time(200.0); // maximum time
terminate: spikes(1);   // stop at first spike
compose: or;            

// using default values of model variables except for dde.alpha 
dde.alpha = 1.1;

// current: is the input(t) to the DDE
sigma = 1.0;                  // std of noise 
current: nse = sigma*eta(t);  // eta(t) is Gaussian white  
stat: first_spiketime;

Arnold Tongue for periodic IPSP input

Compute the Arnold Tongue diagram of the Morris Lecar model with periodic IPSP input 

neuron: morris(20);
dt = morris.dt = 0.1;
dc = morris.dc;
fdc = 0.0607091;    // the unperturbed firing rate
Tdc = 1/fdc;        // the period of the oscillator
num_periods = 200;
terminate: time(num_periods*Tdc);  

// generate the parameter space
N = 150;                       
fmin = 0.0;        
fmax  = 2.5;   // scale by fdc below
ampmin = 0.0;
ampmax = 3.0;
f = linspace(fmin*fdc, fmax*fdc, N);   
amp = linspace(ampmin, ampmax, N);   

// coherent IPSP input
alpha = 3;
train := periodic_point(f, dt);    // periodic spike train at f Hz
resp := alpha*t*exp(-alpha*t)*Heavi(t); // impulse response	
current: s = -amp*dc*delta_conv_efficient(train(t), resp(t), 10/alpha);
s->coherent;                       // every neuron gets this

// independent noise current
sigma = 0.05;
current: nse = dc*sigma*eta(t); // independent by default

lambda = 2*5*fdc;  // the time scale of the sync stat
stat: synchrony(-lambda);

Neuron Model Values

Each neural model defines a number of variables with default values which are loaded and available at the declaration neuron: (see section Neuron Models).

leaky variables

The default values for the leaky integrate and fire mode, (the `leaky.' prefix has been stripped for clarity) written as: 

d/dt V(t) = -V(t)/ (rm cm) + (dc + I(t) ) / cm

where I(t) is the sum of the user defined currents. When V(t) >= theta, a spike is triggered and the membrane potential is reset to vreset.

morris variables

delay variables

The default values for the delay differential equation model (written without the `dde.' prefix for clarity), written as: 

V'(t) = alpha*V(t) + beta*V(t-tau) + I(t)

where I(t) is the sum of the user defined currents. When |V(t)| >= theta, a spike is triggered but there is no reset.

delay2 variables

The default values for the delay differential equation model (written without the `dde2.' prefix for clarity), written as: 

$$V''(t)=alpha*V(t)+beta*V(t-\tau)+gamma*V'(t)+delta*V'(t-\tau)+I(t)$$

where I(t) is the sum of the user defined currents. When |V(t)| >= theta, a spike is triggered but there is no reset.

This document was generated on 14 June 2000 using texi2html 1.56k.