Semi-explicit equation

Support
1

Description of problem

I have a stimuli function which cannot be cast as a simple logical expression. I have interpolated it using scipy.interpolate.interp1d to make a dimensionless interpolator object, say I(t). I wished Brian could handle this semi-explicit input, but apparently it doesn’t.

Is there a way to pass it to a neuron model? When I do so, even after multiplying the correct unit, I get the following error:
TypeError: Object of type <class 'scipy.interpolate._interpolate.interp1d'> does not have dimensions

As a follow-up, how can I do the same if I(t) returns an N-dimensional array, where N is the size of the network? (The use case would be applying a non-uniform time-dependent input to a network).

Minimal code to reproduce problem

import numpy as np
import brian2 as b2

b2.start_scope()

model = """
dv/dt = -v + I(t) : 1
I: Hz
"""

N = 10

def stim(t):
    return np.sin(t)*np.arange(N)
    
G = b2.NeuronGroup(N, model,)
G.I = stim*b2.Hz

Which returns
TypeError: Object of type <class 'function'> does not have dimensions

Thanks in advance for your help!

Best, Arash

2

There is a general mechanism to create new functions for use in Brian. It requires you to annotate the definition with a @check_units decorator so that Brian can verify the units for consistency (see Functions — Brian 2 2.5.1 documentation). However, I think what you are looking for is a TimedArray, see Input stimuli — Brian 2 2.5.1 documentation
It can be a one-dimensional array of time, or a two-dimensional array for a function of time and neuron index.

3

Thanks, Marcel!
My input is changing continuously with time. I guess using TimedArray for them would cause a lot of overhead as the dt argument should be equal to the simulation dt.
I checked the Functions documentation and adjusted my code with a decorator, but still get the same error.

import numpy as np
import brian2 as b2
from brian2 import check_units, implementation, Function

print(b2.__version__)  # out: '2.5.1'

b2.start_scope()

model = """
dv/dt = -v/ms + I : 1
I: Hz
"""

N = 10

@implementation('numpy', discard_units=True) # seemingly doesn't matter
@check_units(t=b2.second, result=b2.Hz)
def stim(t):
    return np.sin(t/b2.ms)*np.arange(N)*b2.Hz
    
fstim = Function(stim, arg_units=[b2.second], return_unit=b2.Hz)

G = b2.NeuronGroup(N, model,)

# both of the following produce the same error
G.I = stim
G.I = fstim

# raised error: TypeError: Object of type <class 'function'> does not have dimensions

Can you let me know what am I missing here? Many thanks!
Arash

4

Hi @arashgmn,

My input is changing continuously with time. I guess using TimedArray for them would cause a lot of overhead as the dt argument should be equal to the simulation dt.

Changing every dt is ok, this is not less efficient than using a function that gets called every time step – TimedArray is implementing a function that efficiently looks up values in an array, if you manually implement the same thing it will most likely not be more efficient than that. Also note that for user-defined functions, you will have to define them in the target language, e.g. if you want your code to work with Cython, you need a Cython implementation, if you want your simulation to run in C++ standalone mode, you need a C++ implementation.

I checked the Functions documentation and adjusted my code with a decorator, but still get the same error.

Your code is not using the function in the right way. You are trying to assign the function as an initial value for the variable I, but it expects a value, not a function. The error message is a bit misleading, since it tries to check whether the value that you provide has the right units, but it is not a value and therefore does not have any units.
The correct way to use your function would be directly in the equations, e.g. with:

model = """
dv/dt = -v/ms + I : 1
I = stim(t, i): Hz
"""

Also note that the function needs a second argument for the neuron indices. It cannot simply unconditionally return an array of size N, since it might be called in other contexts, e.g. in a reset statement. See Functions — Brian 2 2.5.4 documentation for more details.

But again, I think TimedArray does exactly what you want, so I wouldn’t bother with all this.

1 Like
5

@mstimberg , thank you for your detailed instructions. As you suggested, I switched to TimedArray and I believe the issue I had earlier is resolved. However, down the line, cython complains about the TimedArray variable I_stim that I pass to the b2.SpatialNeuron class through a namespace. The error reads:

Error compiling Cython file:
------------------------------------------------------------
...
    cdef double I_point
    # namespace for function I_stim
    global _namespace_timedarray_values
    global _namespace_num_timedarray_values
    cdef _numpy.ndarray[double, ndim=1, mode='c'] _buf__timedarray_values = _namespace['_timedarray_values']
    _namespace_timedarray_values = <double *> _buf__timedarray_values.data
                                  ^
------------------------------------------------------------

You can reproduce the error using this minimal example:

import brian2 as b2
from brian2.units import (ms, us, second, 
                          amp, nA,  uA, pA, 
                          meter, mm, um, 
                          volt, mV, ohm, siemens, farad)

b2.defaultclock.dt = 50*b2.us
dt = b2.defaultclock.dt 
b2.start_scope()

def simulate(as_timedarray=False):
    print('Started simulation with I as_timedarray = ' + str(as_timedarray))

    if as_timedarray:
        I_stim = b2.TimedArray([1]* uA/mm**2, dt=dt) 
        model='''Im = -gl * (v-El) + I_stim(t) : amp/meter**2'''
    else:
        I_stim = 1*uA/mm**2 
        model ='''Im = -gl * (v-El) + I_stim : amp/meter**2'''
        
            
    # some morphology (taken from the docs)
    morpho = b2.Section(n=5, 
                        diameter=[15, 5, 10, 5, 10, 5]*um, 
                        length=[10, 20, 5, 5, 10]*um)
    
    namespace = {
        'El': -84. * b2.mvolt,
        'gl': 400. * siemens / (meter ** 2),
        'I_stim': I_stim,
        }
    
    # Setting up a spatial neuron
    nrn = b2.SpatialNeuron(morphology=morpho, model=model,
                           namespace=namespace ,
                           Cm=0.028 * farad / (meter ** 2),
                           Ri=0.4 * ohm * meter,
                           method='exponential_euler', 
                           refractory="v > -50*mV", threshold="v > 0*mV",
                           reset=''
                           )
    
    nrn.v = namespace['El'] # initialization
    
    # Setting up a monitor
    mon = b2.StateMonitor(nrn, ['v'], True) 
    
    # Setting up a network
    net = b2.Network()
    net.add(nrn)
    net.add(mon)
    
    # running simulation
    net.run(100*ms) 
    
    print('Simulation finished with I as_timedarray = ' + str(as_timedarray))

simulate(False) # runs just fine
simulate(True)  # doesn't run :(

I appreciate your help.

Cheers, Arash

6

Hi @arashgmn. This is this bug here: Compilation error for SpatialNeuron, TimedArray and Cython · Issue #1427 · brian-team/brian2 · GitHub

It is already fixed in the latest release. If you use pip, you can upgrade with pip install --upgrade brian2. If you installed brian2 with conda, you’ll have to wait a bit until the conda package has been updated (the missing conda package is the reason I didn’t officially announce the release yet).

In the meantime, you can also work around the issue by switching to Python mode with

b2.prefs.codegen.target = 'numpy'