Skip to content

Model fitting

Jump to bottom
romainbrette edited this page May 19, 2016 · 32 revisions

There was a model fitting toolbox in Brian 1 but it has not been ported. It used vectorized and parallelized simulation on GPU and multiple CPUs, using Playdoh. It was mainly for fitting models to spike trains and used non-gradient optimization algorithms (e.g. CMAES). It might be possible to port it to Brian 2. But it would be nice in any case to have a simple toolbox to optimize models to continuous traces. A typical use case is to fit channel models (Hodgkin-Huxley or Markov models) to voltage-clamp recordings. Or to optimize a entire neuron model to current-clamp recordings.

Current status

In the modelfitting branch, there is a basic implementation of a model fitting function, using a vectorized version of Scipy's differential evolution algorithm. It only does least-square fitting. To deal with initialization of state variables, argument t_start can be provided: error is only counted after that time.

Typical use scenarios

Trace fitting

A typical dataset is a set of voltage-clamp recordings of pharmacologically isolated currents, for example Na currents. The command (input) is typically a sequence of voltage steps, and the recordings are currents. Sometimes the command can be an action potential waveform. Typically we have a parametric model of the current that we want to fit to the data. Everyone uses least-square fitting, either on the current or on the conductance. But we might well want to use other criteria, and in particular we might want to put more weight on some part of the data than on others. For example if we are interested in spike initiation, then we might want to put more weight on low-voltage data. So generally we have:

  • A set of voltage and corresponding current traces.
  • A criterion, defined by: a distance function between currents (which would allow quadratic or Lp criteria, for example), and possibly weights for the different traces (alternatively, a weight function of voltage, but this is less general).
  • A model, defined as an Equations object with parameters.
  • A target variable to fit (e.g. current I).
  • Constraints on the parameters.

Here is an example of a model:

eqs=Equations('''
I = -g*n**8*(EK-v) : amp
dn/dt = alphan*(1-n) - betan*n : 1
alphan = a*(v-vn)/(exp(d*(v-vn))-1) : 1/second
betan = b*exp(c*v) : 1/second

a : 1/volt/second (constant)
b : 1/second (constant)
c : 1/volt (constant)
d : 1/volt (constant)
g : siemens (constant)
vn : volt (constant)
''')

Here EK is an external, fixed variable. The input command is v.

Model optimization

We want to optimize the parameters of a model with respect to an intrinsic parameter. For example, we might want to look for a model that minimizes energy consumption, under some constraints. Constraints could be inserted as part of the optimization criterion. In this case, there could be no input. But typically the optimization criterion could be complex.

What is needed for this kind of cases is a set of convenience functions to run the same model and criterion calculation over many parameter sets, using vectorization and/or parallelization, that is, a map function for models (or a network), which should include operations and monitors.

Optimization and fitting of multicompartmental models

It would be useful to be able to run optimizations for multicompartmental models too. In this case, probably parallel runs rather than vectorization would be used.

Spike train prediction

This is the typical use case for the model fitting toolbox in Brian 1.

Initialization

An important issue is initialization. There are two ways to deal with this:

  1. Initial values are fixed (e.g. to 0). The fitting criterion applies only after a certain duration.
  2. Initial values are calculated before running the model, as the steady-state value for a given value of the input. For example in the case above, we would use the first value v(0*ms,i) and calculate the steady-state value of I. This can be done in several ways: running the model for a long time (in fact for this specific case, a single very large time step with exponential Euler would work), using a symbolic equation solver, or using a nonlinear root solver.

Obviously the first idea is simpler.

Pre-processing

Traces are often quite noisy. We could ignore this, or we might want to provide ways to reduce noise. The spectrum of the noise could be analyzed on the initial part of the recordings, and used to filter the traces.

It would be nice to be able to use a relative criterion, for example ((I-Imeasured)/I_measured)**2, but this won't be feasible if the traces are too noisy. Apart from filtering the traces, we could also preprocess traces by curve fitting, e.g. fitting splines to each individual trace.

Optimization algorithms

Non-gradient population-based algorithms

The simplest way would be to use algorithms that are not based on gradients. Preliminary work shows that Scipy's differential_evolution gives ok results. Another possibility is the DEAP package.

However, Scipy uses a sequential algorithm, and this is costly. We would need to edit the code to allow vectorized simulation of the population, and perhaps even multi-core simulation.

Gradient algorithms

Gradient algorithms might be adapted for this kind of problem, and they are generally faster. But in this case we need to automatically calculate gradients. This is actually possible using sympy, although not trivial. Then we can directly plug into Scipy's optimization algorithms with constraints.

General procedure

  1. First, the user passes the equations object, the input and output traces and variable names. The criterion is inserted in the equations (maybe with special name E or criterion). Bounds for the parameters are also passed.
  2. TimedArrays are built from the traces, and the equations are edited (v is replaced by v(t,i)). All traces are simultaneously simulated. This means they must all have the same duration. If not, they are extended to have all the same duration.
  3. For gradient-based algorithms, gradients are calculated symbolically and inserted into the equations.
  4. For population-based algorithms, neurons are duplicated: v is replaced by v(t,i % N), where N is the number of traces.
  5. The criterion is calculated either by integration using an equation, or using a network operation (just summing).
  6. An evaluation function is created that initializes the model, runs it, retrieves the criterion value.
  7. We call an optimization algorithm, or a modified distributed evolution algorithm for vectorized computation.

Technical issues

  • Ideally we would want to use the standalone mode, but a major limitation currently is that it is not trivial (possible?) to run a model several times with different parameter values. This needs to be solved.

Syntax

The syntax is inspired from Brian 1 model fitting toolbox. A first proposition:

model = Equations('''
I = -g*n**8*(EK-v) : amp
dn/dt = alphan*(1-n) - betan*n : 1
alphan = a*(v-vn)/(exp(c*(v-vn))-1) : 1/second
betan = b*exp(c*v) : 1/second

a : 1/volt/second
b : 1/second
c : 1/volt
vn : volt
g : siemens''')

params = fit_traces(model, input_var = 'I', output_var = 'v',
                    input = input_traces, output = output_traces,
                    g = [1*nS,10*nS], alpha_n = [1*Hz,10*Hz],...)

Existing tools

Neuron has a model fitting tool. There are a few other programs that optimize models to traces; some of them interface with Neuron. Here are a few:

and related efforts:

Clone this wiki locally