Model fitting API

[romainbrette]

We have some model fitting code in the model fitting branch. It uses differential evolution to fit traces. But what we could do is a general API, which would work with any optimization function. The user simply defines the criterion to be calculated, and passes the optimization function to use. This could work for any gradient-free algorithm.

What we need to add is a callback for initializing the models, and maybe a few other things.

Then we might be able to use gradient-based algorithms by using automatic differentiation tools such as Autograd. To calculate dE/dp, where E is an error criterion and p is a parameter, we need to be able to calculate dVm/dp. We get this by differentiating the membrane equation with respect to p, and then integrating. We see that we actually need to calculate dXi/dp for every state variable i. We use autograd to calculate the RHS. The main question is whether this can be easily vectorized across neurons.

[romainbrette]

Some optimization/model fitting tools:

• DPSE
• Nevergrad
• Neurofit
• symfit

Other refs on model fitting:

• Brookings et al. (2014) applies Abarbanel's method.
• A recent paper by Abarbanel: Nogaret et al (2016)
• A recent paper using Approximate Bayesian Computation: Lueckmann et al (2017), comes with a Python package.
• Vavoulis et al (2012).

[romainbrette]

It is probably possible to use also gradient-based optimization algorithms. Perhaps we could use Autograd to compute gradients of the dynamical system (but then without code generation), and then use gradient descent or some other optimization procedure with derivatives. If there are n variables and p parameters, then one must simulate n*p variables.

To calculate dE/dp, where E is an error criterion and p is a parameter, we need to be able to calculate dVm/dp. We get this by differentiating the membrane equation with respect to p, and then integrating. We see that we actually need to calculate dXi/dp for every state variable i. We use autograd to calculate the RHS. The main question is whether this can be easily vectorize.

[mstimberg]

I think the biggest challenge here is to find a good balance between reusing existing libraries and at the same time make good use of Brian's code generation facilities (e.g. in your branch, the error is calculated within Brian code instead of externally).

About gradient-based methods: I don't see a use case for Autograd here, if I am not mistaken its strength lies in being able to do derivatives of (almost) arbitrary Python/numpy code, i.e. code involving loops, if statements, and so on. In our case, we only have equations and can do everything symbolically with sympy I think?

Plenty of things to discuss with @alTeska who will work on this in the coming months :smile:

[romainbrette]

You're right about the differentiation, I think! It would be awesome to be able to used gradient-based optimization methods.

[alTeska]

I have extended the list of possible optimization/model fitting tools listed them on the new Project Board if you'd like to have a look. Will make a proper review soon.

As it comes to gradient-based optimization they have been reported to be used for local optimization, used to refine the results, so if possible we should include them in the module.

[thesamovar]

Should this all perhaps go in a separate package (e.g. https://github.com/brian-team/brian2modelfitting) rather than being added to brian2tools?

[mstimberg]

Should this all perhaps go in a separate package (e.g. https://github.com/brian-team/brian2modelfitting) rather than being added to brian2tools?

This would definitely be an option, especially if it the package gets very complex/powerful. On the other hand, I think it would help for the discoverability if there is a single package brian2tools, instead of all the functionality spread out over several packages. A very practical consideration is also to avoid having yet another test suite and readthedocs setup...

But changing from one to the other should be pretty easy later on, we just have to take a final decision before we make a release.

[romainbrette]

One thing I was just thinking about for the API is the possibility to deal with recordings done in different conditions (e.g. drugs). I see two ways of dealing with this:

• having parameters that can be set by the user for each trace (e.g. channel on/off));
• having parameters fitted by the program that can differ between groups of traces. For example, let's say we have conductance gNa as a parameter. There are (say) three conditions with three different levels of TTX, corresponding to three unknown values of gNa. The program would have to fit 3 parameters gNa1, gNa2, gNa3. Technically, this could be done for example by adding a line gNa = a gNa1+b gNa2+c * gNa3, and (a,b,c) are binary parameters set differently for each trace.

Syntax requires some thought! This is something we can think about later anyway.