Model fitting

The modelfitting library is used for fitting a neuron model to data.

The library provides a single function modelfitting(), which accepts the model and the data as arguments and returns the model parameters that fit best the data. The model is a spiking neuron model, whereas the data consists of both an input (time-varying signal, for example an injected current) and a set of spike trains. Only spikes are considered for the fitness. Several target spike trains can be specified in order to fit independently several data sets. In this case, the modelfitting() function returns as many parameters sets as there are target spike trains.

The model is defined as any spiking neuron model in Brian, by giving the equations as mathematical equations, and the reset and threshold values. The free parameters of the model that shall be fitted by the library are also specified. The data is specified by the input (a vector containing the time-varying injected current), the timestep of the input, and the data as a list of spike times.

How it works

Fitting a spiking neuron model to electrophysiological data is performed by maximizing a fitness function measuring the adequacy of the model to the data. This function is defined as the gamma factor, which is based on the number of coincidences between the model spikes and the experimentally-recorded spikes, defined as the number of spikes in the experimental train such that there is at least one spike in the model train within plus or minus delta, where delta is the size of the temporal window (typically a few milliseconds). For more details on the gamma factor, see Jolivet et al. 2008, “A benchmark test for a quantitative assessment of simple neuron models”, J. Neurosci. Methods (available in PDF here).

The optimization procedure is performed by an optimization algorithm. The optimization toolbox used by modelfitting is implemented in the external Python package Playdoh. It also supports distributed and parallel optimization across CPUs and machines. Different optimization algorithms are supported, the default one is CMAES. All those algorithms require the evaluation of the fitness function for a large number of parameter sets. Each iteration of the algorithm involves the simulation of a large number of neurons (one neuron corresponding to one parameter set) as well as the computation of the gamma factor for each neuron. The quality of the result depends on the number of neurons used, which is specified in the modelfitting() function.

Playdoh supports the use of graphical processing units (GPUs) in order to accelerate the speed of convergence of the algorithm. If multiple cores are detected, the library will use all of them by default. Also, if a CUDA-enabled GPU is present on the system, and if PyCUDA is installed, the library will automatically use the GPU by default. In addition, several computers can be networked over IP, see Clusters.

Usage example

To import the library, use

from brian.library.modelfitting import *

To fit the parameters of a neuron model with respect to some data, use the modelfitting() function

results = modelfitting(model = equations, reset = 0, threshold = 1,
                     data = spikes,
                     input = input, dt = .1*ms,
                     popsize = 1000, maxiter = 10,
                     R = [1.0e9, 1.0e10], tau = [1*ms, 50*ms])

print_table(results)

Warning

Windows users should read the section Important note for Windows users.

The model is defined by equations (an Equations object), reset (a scalar value or a set of equations as a string) and threshold (a scalar value or a set of equations as a string).

The target spike trains are defined by data (a list of pairs (neuron index, spike time) or a list of spike times if there is only one target spike train).

The input is specified with input (a vector containing the time-varying signal) and dt (the time step of the signal). The input variable should be I in the equations, although the input variable name can be specified with input_var.

The number of particles per target train used in the optimization algorithm is specified with popsize. The total number of neurons is popsize multiplied by the number of target spike trains. The number of iterations in the algorithm is specified with maxiter.

Each free parameter of the model that shall be fitted is defined by two values

param_name = [min, max]

param_name should correspond to the parameter name in the model equations. min and max specify the initial interval from which the parameter values will be uniformly sampled at the beginning of the optimization algorithm. A boundary interval can also be specified by giving four values

param_name = [bound_min, min, max, bound_max]

The parameter values will be forced to stay inside the interval [bound_min, bound_max] during the optimization.

The complete list of arguments can be found in the reference section of the modelfitting() function.

The best parameters and the corresponding best fitness values found by the optimization procedure are returned in the OptimizationResult object result.

Important note for Windows users

The model fitting library uses the Python multiprocessing package to distribute fitting across processors in a single computer or across multiple computers. However, there is a limitation of the Windows version of multiprocessing which you can read about here. The end result is that a script like this:

from brian.library.modelfitting import *
...
results = modelfitting(...)

will crash, going into an endless loop and creating hundreds of Python processes that have to be shut down by hand. Instead, you have to do this:

from brian.library.modelfitting import *
...
if __name__=='__main__':
        results = modelfitting(...)

Clusters

The model fitting package can be used with a cluster of computers connected over IP. Every computer must have Brian and Playdoh installed, and they must run the Playdoh server: see the Playdoh documentation. Then, you can launch the `modelfitting function with the machines keyword, which is the list of the IP addresses of the machines to use in parallel for the fitting procedure. You must also specify the unit_type keyword, which is CPU or GPU, to indicate whether you want to use CPUs or GPUs on these computers. You can’t mix CPUs and GPUs for the same optimization.

IP

To connect several machines via IP, pass a list of host names or IP addresses as strings to the machines keyword of the modelfitting() function. To specify a specific port, use a tuple (IP, port) instead of a string. You can also specify a default port in the Playdoh user preferences, see the Playdoh documentation.

Authentication

You can specify an authentication string on all the computers running the Playdoh server to secure communications. See the Playdoh documentation.

Example

The following script launches a fitting procedure in parallel on two machines:

from brian import loadtxt, ms, Equations
from brian.library.modelfitting import *

if __name__ == '__main__':
    # List of machines IP addresses
    machines = ['bobs-machine.university.com',
                'jims-machine.university.com']

    equations = Equations('''
        dV/dt=(R*I-V)/tau : 1
        I : 1
        R : 1
        tau : second
    ''')
    input = loadtxt('current.txt')
    spikes = loadtxt('spikes.txt')
    results = modelfitting( model = equations,
                            reset = 0,
                            threshold = 1,
                            data = spikes,
                            input = input,
                            dt = .1*ms,
                            popsize = 1000,
                            maxiter = 3,
                            delta = 4*ms,
                            unit_type = 'CPU',
                            machines = machines,
                            R = [1.0e9, 9.0e9],
                            tau = [10*ms, 40*ms],
                            refractory = [0*ms, 10*ms])
    print_table(results)

The two remote machines would run the Playdoh server.