Model fitting with refractory period

I’m following the model fitting example of Tutorial: TraceFitter — Brian2modelfitting 0.3 documentation, but I’m using an adaptive exponential model instead of the given HH model. Now I want to introduce also refractory period, but I get an error. Here’s my code:

from brian2modelfitting import *
from brian2 import *
import pandas as pd
inp_trace = pd.read_csv('input_traces_hh.csv', index_col=0).to_numpy()
out_trace = pd.read_csv('output_traces_hh.csv', index_col=0).to_numpy()
area = 20000*umetre**2                                                                                                                                                                                                                    

model = '''dv/dt = (gL*(EL-v)+gL*DeltaT*exp((v-VT)/DeltaT) -w + I)/Cm : volt (unless refractory)                                                                                                                                                                                  
  dw/dt=(a*(v-EL)-w)/tauw : amp                                                                                                                                                                                                                                                   
  gL : siemens (constant)                                                                                                                                                                                                                                                         
  EL : volt (constant)                                                                                                                                                                                                                                                            
  DeltaT :  volt (constant)                                                                                                                                                                                                                                                       
  VT : volt (constant)                                                                                                                                                                                                                                                            
  a : siemens (constant)                                                                                                                                                                                                                                                          
  tauw : second (constant)                                                                                                                                                                                                                                                        
  wadd : ampere (constant)                                                                                                                                                                                                                                                        
  vre : volt (constant)                                                                                                                                                                                                                                                           
  Cm : farad (constant)                                                                                                                                                                                                                                                           
  tref : second (constant)                                                                                                                                                                                                                                                        
'''

opt = NevergradOptimizer()
metric = MSEMetric()
fitter = TraceFitter(model=model,
                     input_var='I',
                     output_var='v',
                     input=inp_trace * amp,
                     output=out_trace*mV,
                     dt=0.01*ms,
                     n_samples=1000,
                     method='exponential_euler',
                     param_init={'v': -65*mV},
                     reset='v=vre; w += wadd',
                     threshold='v > 0*mV',
                     refractory='tref')

res, error = fitter.fit(n_rounds=10,
                        optimizer=opt,
                        metric=metric,
                        Cm=[100*pfarad,300*pfarad],
                        gL=[2*psiemens, 300*nsiemens],
                        EL=[-90*mV,-40*mV],
                        VT=[-90*mV,-40*mV],
                        DeltaT=[0.1*mV,20*mV],
                        a=[-1.5*nS,3.0*nS],
                        tauw=[1.0*ms,1000.0*ms],
                        vre=[-80*mV,-60*mV],
                        wadd=[-0.2*nA,0.5*nA],
                        tref=[0.0*ms,5.0*ms])
traces = fitter.generate_traces()

The code works if I remove “(unless refractory)”, that is, remove the refractory period, but otherwise it gives me a division by zero. Any ideas how to fix this issue? Is refractory period supported at all in brian2modelfitting library?

I also seem to get NaNs in the objective function values when I remove “(unless refractory)”, unless I make the length of the refractory period shorter than the simulation time step.

Hi @Tuoma . Using refractoriness in fitted models should work in principle. The error you are getting is due to a bug, though it seems to be a bug in Brian’s code optimization itself rather than in brian2modelfitting. Brian’s code generation simplifies expression by replacing expressions with fixed values, e.g. 1.0 / 2.0 will be replaced by 0.5. Apparently in this case it tries to replace 1.0 / 0.0 by a value, which leads to the division by zero error. I will have to have a closer look into why this is happening in this specific case, but in the mean time I think you could use the following work-around. Instead of (unless refractory), you use dv/dt = not_refractory*(gL*....), i.e. manually multiply the RHS of the equation by 0 or 1. This is basically what (unless refractory) does, even though in general it does a bit more (e.g. discard changes to the membrane potential by synapses) – but this should not matter here.

Thanks, that seems to work!

1 Like

For future reference, this has now been fixed in Brian’s current development version: Make evaluate_expr more robust by mstimberg · Pull Request #1389 · brian-team/brian2 · GitHub