clean up benchmark figure

[sg-s]

• [x] see the template in figure_benchmark_2
• [x] see things to do in figure below

[sg-s]

working on this now...

[sg-s]

multiple problems with DynaSim @alec-hoyland

1. you're using rk2, not exponential euler -- is it not possible to specify how you want the eq.s integrated?
2. for the biggest time step, voltage trace is all NaN (model explodes)

[sg-s]

@alec-hoyland does dynasim run in closed loop mode?

[sg-s]

this pleases me: xolotl is upto 10X faster than dynasim, for the same quality of simulation (similar errors)

[alec-hoyland]

@sg-s Yes, DynaSim supports multiple integration methods. Specifically, they support Runge-Kutta 2nd, and 4th, forward euler, and any built-in MATLAB solver. See here for recommended options.

[sg-s]

what's the key for forward euler? i think that's the fairest comparison

[alec-hoyland]

'solver'      : solver for numerical integration (see dsGetSolveFile)
                     {'euler','rk2','rk4', or any built-in matlab solver} (default: 'rk4')

[alec-hoyland]

I believe that DynaSim operates in open-loop, since it requires initial conditions specified in the equations block and doesn't have a structured way of saving output (except as an output, which is quickly forgotten).

[sg-s]

OK, good to know. maybe we should make a point of how easy it is to run xolotl in closed loop

[alec-hoyland]

I think it's a pretty significant feature. xolotl saves the state of the simulation in the struct, meaning that you can scope in and look at any specific value pretty intuitively. It also means that it is present during the next simulation without you having to think about it. If you want more than temporary storage, you can save the state and reload it later.

In comparison to functional-programming paradigms, where you would have to manually pass a data structure of conditions back to your simulation function.

[alec-hoyland]

DynaSim might have the ability to run in closed loop, but if it does, I can't find it.

[sg-s]

any progress on doing this with NEURON? also, when you do the NEURON benchmarks, you have to do two sets:

1. a "fair" comparison (fixed time step)
2. allow NEURON to do whatever it wants -- variable time step, etc. and see how fast it is

[alec-hoyland]

The "fair" comparison should be all set. That is, if you run the files in xolotl-paper/neuron/HH/ and ../neuron/STG/, you should generate .csv files which contain the output vectors. The raw voltage traces need to be processed in MATLAB before an accuracy measure can be plotted. The speed vectors should be all set. Obviously, these are loaded into MATLAB with csvread.

I haven't tested these scripts because I do not have NEURON on my laptop. I could certainly try to get it on here, but I hate it with my heart and soul. Maybe sshing into chaos is the best course of action...would have to figure out how to set that up from afar.

With great difficulty in actually acquiring this information, it seems that vanilla NEURON uses a static time step with either implicit euler or crank-nicholson. Perhaps it uses an exponentially discretized version of CN called cnexp, but fuck me if I actually know what equations it's solving...

Variable time-steps are implemented with CVODE using SUNDIALS. It seems the way to do this is to integrate using

cvode = h.CVode()
cvode.solve(t_end)
cvode.statistics()

[alec-hoyland]

First though, a custom version of NEURON needs to be compiled

cd ~/code/xolotl-paper/neuron/conductances
nrnivmodl acurrent.mod cas.mod cat.mod hcurrent.mod kca.mod kd.mod na.mod

Then

cd ../HH
nrniv -python benchmark1.py

[alec-hoyland]

Ugh. Biting the bullet. Going to install NEURON and break my Python distribution.

[alec-hoyland]

NEURON works on Arch Linux...at least, it hasn't thoroughly broken yet. Significant changes to .mod files here ( #9 ). I've realized that NEURON tries to couch the fact that C is hard to write well by hijacking MODL and representing models in NMODL, which isn't a real language, since it isn't run as a script with JIT compilation, nor is it compiled and saved as an executable. It's more of a format than anything else...which means that I'm using NMODL to make my C easier to write and then calling my C from Python to make it easier to analyze...all without a list of all the properties I can put in my NMODL mechanisms. For example, cai is a reserved keyword which represents the intracellular calcium concentration, such as in the line READ ca WRITE cai (like this is some kind of SQL database).

[alec-hoyland]

• [x] HH benchmark 1 passing

• [x] HH benchmark 2 passing

• [x] HH benchmark 3 passing

• [x] STG benchmark 1 passing

• [x] STG benchmark 2 passing

• [x] STG benchmark 3 passing

[alec-hoyland]

I don't understand how to get CVODE to work. This code should at the very least integrate an HH model for 10 seconds. Whether it saves the data properly is another question, but there should be no indication that the simulation (cvode.solve()) takes negligible time (viz. it doesn't actually do anything. It returns an error code of zero, meaning that it supposedly ran successfully.

# import the graphical interface
from neuron import h, gui
# import arrays and graphics
import numpy as np
from matplotlib import pyplot
import time

# create the neuron
soma = h.Section(name='soma');

# set the size of the soma
soma.L = 28.209; # microns
soma.diam   = 28.209; # microns

# set up the capacitance
soma.cm     = 1; # μF/cm^2

# add conductances
soma.insert('pas')
soma.insert('hh')

v_vec       = h.Vector()
t_vec       = h.Vector()
v_vec.record(soma(0.5)._ref_v)
t_vec.record(h._ref_t)

cvode = h.CVode()
cvode.solve(10000)

[sg-s]

if it ran successfully, cna you pull out the voltage trace and look at it?

[alec-hoyland]

The HH model works great. There were some significant problems with the intracellular calcium in the STG model. I have fixed most of them...it's producing a tonically spiking model in NEURON for some reason...

Here is the HH model This test was 30 s of simulation at dt = 0.1 ms. It took 3.21 s for a speed factor of 9.32x.

[sg-s]

well, the liu paper you cited is not the neuron in test_bursting_model they have some funky time constants and the area is also different.

[sg-s]

and about the speed -- was this at a fixed time step?

[alec-hoyland]

Yes, it was at a fixed time-step.

[sg-s]

i think the calcium eq you're using is different from the calcium eq in xolotl

the difference is in your equation, tau divides the term that governs how much calcium enters the cell when calcium channels open. in xolotl, it does not. as you can see in many of the test examples, there is actually no difference b/w these equations, and you can rewrite one form in the other, but it does change what your parameters are.

[alec-hoyland]

Current version of the benchmark figure, without cartoons and missing DynaSim STG. I'm still working on debugging it. It's more difficult than NEURON somehow.

Note: just found a bug. Don't use this figure anywhere.

[sg-s]

very nice!

[alec-hoyland]

DynaSim has some really profound issues that are going to take me a long time to deal with. This is going to make finishing the final panel (bottom right) pretty difficult, though everything else seems to be working.

[sg-s]

yeah, the top right is clearly wrong

[alec-hoyland]

[sg-s]

better...what's going on with F?

[alec-hoyland]

I don't know what's going on with F. I am rerunning all of the simulations from scratch to fix a few small issues...one fun fact: DynaSim caused chaos to run out of memory during the compartment simulations.

DynaSim uses ~30 GB of RAM for 1000 Hodgkin-Huxley compartments. NEURON and xolotl use ~8 GB.

[alec-hoyland]

Solved issue with DynaSim by optimizing code by clearing useless variables. Cut DynaSim down to 19 simulations (the last step in nComps is missing). Benchmark simulations are all cached on chaos.