Custom models

With great power comes great responsibility! - Uncle Ben

In the following tutorial, we’ll explore how to modify neuron models created by Dendrify to generate custom models. Specifically:

• First, we’ll create a simple 2-compartment neuron with adaptation using Dendrify, which will serve as our “control” model.

• Next, we’ll attempt to replicate this model using only custom equations and parameters.

Disclaimer

This tutorial is intended for experienced Dendrify/Brian users who are comfortable working with custom model equations. Modifying Dendrify’s equations offers significant flexibility but also increases the risk of introducing errors. While Dendrify and Brian will catch most of these and raise exceptions, it’s crucial to thoroughly test any custom models before using them in experiments to ensure they work as expected and don’t contain hidden bugs. If you need assistance with creating a custom model, please refer to the available support options.

Importing Libraries & Defining Helper Functions

To avoid repetitive code, we will first create some helper functions that handle the generation of all necessary Brian objects, run the simulation, and plot the results. You don’t need to focus too much on how they work. Instead, concentrate on the code blocks that demonstrate how to modify equations in Dendrify. At the end of the tutorial, we’ll bring everything together, and it will be clearer how the entire pipeline works.

import brian2 as b
from brian2.units import Hz, ms, mV, nS, pF, pA
from dendrify import Soma, Dendrite, NeuronModel


b.prefs.codegen.target = 'numpy'  # faster for simple simulations

def create_neuron_model(soma, dend):
    model = NeuronModel([(soma, dend, 15*nS)], v_rest=-60*mV)
    model.add_params({
        'Vth': -40*mV,
        'tauw_soma': 150*ms,
        'a_soma': 0*nS,
        'b_soma': 50*pA,
        'Vr': -50*mV})
    return model

def create_neuron_group(model):
    return model.make_neurongroup(
        1, method='euler',
        threshold='V_soma > Vth',
        reset='V_soma = Vr; w_soma += b_soma',
        refractory=4*ms)

def create_brian_objects(neuron):
    poisson = b.PoissonGroup(10, rates=20*Hz)
    syn = b.Synapses(poisson, neuron, on_pre='s_AMPA_x_dend += 1; s_NMDA_x_dend += 1')
    syn.connect(p=1)
    mon = b.StateMonitor(neuron, ['V_soma', 'V_dend', 'w_soma'], record=True)
    return poisson, syn, mon

def run_simulation(net):
    b.start_scope()
    b.seed(42)
    net.run(500*ms)

def plot_results(mon):
    time = mon.t / ms
    v_soma = mon.V_soma[0] / mV
    v_dend = mon.V_dend[0] / mV
    w = mon.w_soma[0] / pA

    fig, (ax1, ax2) = b.subplots(2, 1, figsize=[6, 6], sharex=True)
    ax1.plot(time, v_soma, label='soma')
    ax1.plot(time, v_dend, label='dendrite', c='C3')
    ax1.set_ylabel('Voltage (mV)')
    ax1.legend()
    ax2.plot(time, w, color='black', label='w')
    ax2.set_ylabel('Adaptation current (pA)')
    ax2.set_xlabel('Time (ms)')
    ax2.legend()
    fig.tight_layout()
    b.show()

def one_function_to_rule_them_all(soma, dend):
    model = create_neuron_model(soma, dend)
    neuron = create_neuron_group(model)
    poisson, syn, mon = create_brian_objects(neuron)
    net = b.Network(neuron, poisson, syn, mon)
    run_simulation(net)
    plot_results(mon)

Control model behavior

soma = Soma('soma', model='adaptiveIF', cm_abs=200*pF, gl_abs=10*nS)
dend = Dendrite('dend', cm_abs=50*pF, gl_abs=2.5*nS)
dend.synapse('AMPA', tag='x', g=3*nS,  t_decay=2*ms)
dend.synapse('NMDA', tag='x', g=3*nS,  t_decay=60*ms)

one_function_to_rule_them_all(soma, dend)

The upper panel shows how the somatic and dendritic voltages evolve over time in response to random synaptic input on the dendrite. The lower panel illustrates the adaptation current dynamics and how it is influenced by somatic spiking. Next, we’ll try to reproduce the same model behavior using custom equations.

Customizing a Soma

Let’s first inspect the somatic equations that were generated by Dendrify:

print(soma.equations)

dV_soma/dt = (gL_soma * (EL_soma-V_soma) + I_soma - w_soma) / C_soma  :volt
dw_soma/dt = (a_soma * (V_soma-EL_soma) - w_soma) / tauw_soma  :amp
I_soma = I_ext_soma  :amp
I_ext_soma  :amp

Now let’s try to create the same equations starting from a simple passive soma:

soma_custom1 = Soma('soma', cm_abs=200*pF, gl_abs=10*nS)
print(soma_custom1.equations)

dV_soma/dt = (gL_soma * (EL_soma-V_soma) + I_soma) / C_soma  :volt
I_soma = I_ext_soma  :amp
I_ext_soma  :amp

Option 1: Replacing existing equations

custom_eqs = '''
dV_soma/dt = (gL_soma * (EL_soma-V_soma) + I_soma - w_soma) / C_soma  :volt
dw_soma/dt = (a_soma * (V_soma-EL_soma) - w_soma) / tauw_soma  :amp'''

soma_custom1.replace_equations(
    'dV_soma/dt = (gL_soma * (EL_soma-V_soma) + I_soma) / C_soma  :volt',
    custom_eqs)

print(soma_custom1.equations)

dV_soma/dt = (gL_soma * (EL_soma-V_soma) + I_soma - w_soma) / C_soma  :volt
dw_soma/dt = (a_soma * (V_soma-EL_soma) - w_soma) / tauw_soma  :amp
I_soma = I_ext_soma  :amp
I_ext_soma  :amp

one_function_to_rule_them_all(soma_custom1, dend)

Notice how both the equations and the model response are identical to the control model.

Option 2: targeted string replacements and additions

soma_custom2 = Soma('soma', cm_abs=200*pF, gl_abs=10*nS)
soma_custom2.replace_equations(
    "I_soma)",
    "I_soma - w_soma)")

soma_custom2.add_equations(
    'dw_soma/dt = (a_soma * (V_soma-EL_soma) - w_soma) / tauw_soma  :amp')

print(soma_custom2.equations)

dV_soma/dt = (gL_soma * (EL_soma-V_soma) + I_soma - w_soma) / C_soma  :volt
I_soma = I_ext_soma  :amp
I_ext_soma  :amp
dw_soma/dt = (a_soma * (V_soma-EL_soma) - w_soma) / tauw_soma  :amp

one_function_to_rule_them_all(soma_custom2, dend)

Once again, we’ve replicated the control model using more targeted string formatting.

Customizing a Dendrite

Let’s first take a look at the dendritic equations generated by Dendrify. Notice that this time, the equations are more complex, as they describe the AMPA and NMDA synaptic dynamics.

print(dend.equations)

dV_dend/dt = (gL_dend * (EL_dend-V_dend) + I_dend) / C_dend  :volt
I_dend = I_ext_dend + I_NMDA_x_dend + I_AMPA_x_dend  :amp
I_ext_dend  :amp
I_AMPA_x_dend = g_AMPA_x_dend * (E_AMPA-V_dend) * s_AMPA_x_dend * w_AMPA_x_dend  :amp
ds_AMPA_x_dend/dt = -s_AMPA_x_dend / t_AMPA_decay_x_dend  :1
I_NMDA_x_dend = g_NMDA_x_dend * (E_NMDA-V_dend) * s_NMDA_x_dend / (1 + Mg_con * exp(-Alpha_NMDA*(V_dend/mV+Gamma_NMDA)) / Beta_NMDA) * w_NMDA_x_dend  :amp
ds_NMDA_x_dend/dt = -s_NMDA_x_dend/t_NMDA_decay_x_dend  :1

Like in the previous example, we will start from a passive dendritic compartment and then modify its equations using the replace_equations and add_equations methods provided by Dendrify.

dend_custom = Dendrite('dend', cm_abs=50*pF, gl_abs=2.5*nS)
print(dend_custom.equations)

dV_dend/dt = (gL_dend * (EL_dend-V_dend) + I_dend) / C_dend  :volt
I_dend = I_ext_dend  :amp
I_ext_dend  :amp

dend_custom.replace_equations(
    'I_dend = I_ext_dend  :amp',
    'I_dend = I_ext_dend + I_NMDA_x_dend + I_AMPA_x_dend  :amp')

custom_synaptic_eqs = """
I_AMPA_x_dend = g_AMPA_x_dend * (E_AMPA-V_dend) * s_AMPA_x_dend * w_AMPA_x_dend  :amp
ds_AMPA_x_dend/dt = -s_AMPA_x_dend / t_AMPA_decay_x_dend  :1
I_NMDA_x_dend = g_NMDA_x_dend * (E_NMDA-V_dend) * s_NMDA_x_dend / (1 + Mg_con * exp(-Alpha_NMDA*(V_dend/mV+Gamma_NMDA)) / Beta_NMDA) * w_NMDA_x_dend  :amp
ds_NMDA_x_dend/dt = -s_NMDA_x_dend/t_NMDA_decay_x_dend  :1"""

dend_custom.add_equations(custom_synaptic_eqs)

print(dend_custom.equations)

dV_dend/dt = (gL_dend * (EL_dend-V_dend) + I_dend) / C_dend  :volt
I_dend = I_ext_dend + I_NMDA_x_dend + I_AMPA_x_dend  :amp
I_ext_dend  :amp

I_AMPA_x_dend = g_AMPA_x_dend * (E_AMPA-V_dend) * s_AMPA_x_dend * w_AMPA_x_dend  :amp
ds_AMPA_x_dend/dt = -s_AMPA_x_dend / t_AMPA_decay_x_dend  :1
I_NMDA_x_dend = g_NMDA_x_dend * (E_NMDA-V_dend) * s_NMDA_x_dend / (1 + Mg_con * exp(-Alpha_NMDA*(V_dend/mV+Gamma_NMDA)) / Beta_NMDA) * w_NMDA_x_dend  :amp
ds_NMDA_x_dend/dt = -s_NMDA_x_dend/t_NMDA_decay_x_dend  :1

SPOILER ALERT!

Although the custom dendrite has identical equations to the control dendrite, running the same simulation will result in an error. Scroll down to find out why.

one_function_to_rule_them_all(soma, dend_custom)

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/.conda/envs/brian/lib/python3.12/site-packages/brian2/core/network.py:1003, in Network.before_run(self, run_namespace)
   1002 try:
-> 1003     obj.before_run(run_namespace)
   1004 except Exception as ex:

File ~/.conda/envs/brian/lib/python3.12/site-packages/brian2/groups/neurongroup.py:989, in NeuronGroup.before_run(self, run_namespace)
    987 def before_run(self, run_namespace=None):
    988     # Check units
--> 989     self.equations.check_units(self, run_namespace=run_namespace)
    990     # Check that subexpressions that refer to stateful functions are labeled
    991     # as "constant over dt"

File ~/.conda/envs/brian/lib/python3.12/site-packages/brian2/equations/equations.py:1165, in Equations.check_units(self, group, run_namespace)
   1163 external -= set(all_variables.keys())
-> 1165 resolved_namespace = group.resolve_all(
   1166     external, run_namespace, user_identifiers=external
   1167 )  # all variables are user defined
   1169 all_variables.update(resolved_namespace)

File ~/.conda/envs/brian/lib/python3.12/site-packages/brian2/groups/group.py:801, in Group.resolve_all(self, identifiers, run_namespace, user_identifiers, additional_variables)
    800 for identifier in identifiers:
--> 801     resolved[identifier] = self._resolve(
    802         identifier,
    803         user_identifier=identifier in user_identifiers,
    804         additional_variables=additional_variables,
    805         run_namespace=run_namespace,
    806     )
    807 return resolved

File ~/.conda/envs/brian/lib/python3.12/site-packages/brian2/groups/group.py:756, in Group._resolve(self, identifier, run_namespace, user_identifier, additional_variables)
    754 # We did not find the name internally, try to resolve it in the external
    755 # namespace
--> 756 return self._resolve_external(identifier, run_namespace=run_namespace)

File ~/.conda/envs/brian/lib/python3.12/site-packages/brian2/groups/group.py:890, in Group._resolve_external(self, identifier, run_namespace, user_identifier, internal_variable)
    889             error_msg = f'The identifier "{identifier}" could not be resolved.'
--> 890         raise KeyError(error_msg)
    892 elif len(matches) > 1:
    893     # Possibly, all matches refer to the same object

KeyError: 'The identifier "w_NMDA_x_dend" could not be resolved.'

The above exception was the direct cause of the following exception:

BrianObjectException                      Traceback (most recent call last)
Cell In[14], line 1
----> 1 one_function_to_rule_them_all(soma, dend_custom)

Cell In[2], line 53, in one_function_to_rule_them_all(soma, dend)
     51 poisson, syn, mon = create_brian_objects(neuron)
     52 net = b.Network(neuron, poisson, syn, mon)
---> 53 run_simulation(net)
     54 plot_results(mon)

Cell In[2], line 28, in run_simulation(net)
     26 b.start_scope()
     27 b.seed(42)
---> 28 net.run(500*ms)

File ~/.conda/envs/brian/lib/python3.12/site-packages/brian2/core/base.py:346, in device_override.<locals>.device_override_decorator.<locals>.device_override_decorated_function(*args, **kwds)
    344     return getattr(curdev, name)(*args, **kwds)
    345 else:
--> 346     return func(*args, **kwds)

File ~/.conda/envs/brian/lib/python3.12/site-packages/brian2/units/fundamentalunits.py:2652, in check_units.<locals>.do_check_units.<locals>.new_f(*args, **kwds)
   2642             error_message = (
   2643                 f"Function '{f.__name__}' "
   2644                 "expected a quantity with unit "
   2645                 f"{unit} for argument '{k}' but got "
   2646                 f"'{value}'"
   2647             )
   2648             raise DimensionMismatchError(
   2649                 error_message, get_dimensions(newkeyset[k])
   2650             )
-> 2652 result = f(*args, **kwds)
   2653 if "result" in au:
   2654     if isinstance(au["result"], Callable) and au["result"] not in (
   2655         bool,
   2656         np.bool_,
   2657     ):

File ~/.conda/envs/brian/lib/python3.12/site-packages/brian2/core/network.py:1141, in Network.run(self, duration, report, report_period, namespace, profile, level)
   1138 if namespace is None:
   1139     namespace = get_local_namespace(level=level + 3)
-> 1141 self.before_run(namespace)
   1143 if len(all_objects) == 0:
   1144     return  # TODO: raise an error? warning?

File ~/.conda/envs/brian/lib/python3.12/site-packages/brian2/core/base.py:346, in device_override.<locals>.device_override_decorator.<locals>.device_override_decorated_function(*args, **kwds)
    344     return getattr(curdev, name)(*args, **kwds)
    345 else:
--> 346     return func(*args, **kwds)

File ~/.conda/envs/brian/lib/python3.12/site-packages/brian2/core/network.py:1005, in Network.before_run(self, run_namespace)
   1003             obj.before_run(run_namespace)
   1004         except Exception as ex:
-> 1005             raise BrianObjectException(
   1006                 "An error occurred when preparing an object.", obj
   1007             ) from ex
   1009 # Check that no object has been run as part of another network before
   1010 for obj in all_objects:

BrianObjectException: Error encountered with object named 'neurongroup_1'.
Object was created here (most recent call only, full details in debug log):
  File '/home/michalis/.conda/envs/brian/lib/python3.12/site-packages/dendrify/neuronmodel.py', line 362, in make_neurongroup
    group = NeuronGroup(N,

An error occurred when preparing an object. (See above for original error message and traceback.)

Although we matched the default dendrite equations, we didn’t provide values for some key parameters. This can be fixed using the add_params method. See below:

dend_custom_params = {
    't_AMPA_decay_x_dend': 2*ms,  # tau AMPA
    't_NMDA_decay_x_dend': 60*ms, # tau NMDA
    'g_AMPA_x_dend': 3*nS,        # AMPA conductance
    'g_NMDA_x_dend': 3*nS,        # NMDA conductance
    'w_NMDA_x_dend': 1,           # NMDA weight
    'w_AMPA_x_dend': 1,           # AMPA weight
    }
dend_custom.add_params(dend_custom_params)

one_function_to_rule_them_all(soma, dend_custom)

Creating and simulating a custom model from scratch

Note

When creating a custom model from scratch, we have the flexibility to simplify certain equations or variable names. This is because we know in advance the exact model we want to build. In contrast, Dendrify is agnostic to the user’s needs and must account for all possible model features to avoid breaking or causing simulation errors. This is why you might encounter some ugly or long variable or parameter names when using stock Dendrify.

# Create and customize soma
soma = Soma('soma', cm_abs=200*pF, gl_abs=10*nS)
soma.replace_equations("I_soma)", "I_soma - w)")
soma.add_equations('dw/dt = (a * (V_soma-EL_soma) - w) / tauw  :amp')

# Create and customize dendrite
dend = Dendrite('dend', cm_abs=50*pF, gl_abs=2.5*nS)
dend.replace_equations(
    'I_dend = I_ext_dend  :amp',
    'I_dend = I_ext_dend + I_NMDA + I_AMPA  :amp')
custom_synaptic_eqs = """
I_AMPA = g_AMPA * (E_AMPA-V_dend) * s_AMPA  :amp
ds_AMPA/dt = -s_AMPA / t_AMPA_decay  :1
I_NMDA = g_NMDA * (E_NMDA-V_dend) * s_NMDA / (1 + Mg_con * exp(-Alpha_NMDA*(V_dend/mV+Gamma_NMDA)) / Beta_NMDA)  :amp
ds_NMDA/dt = -s_NMDA/t_NMDA_decay  :1"""
dend.add_equations(custom_synaptic_eqs)

# Create neuron model and add parameters
model = NeuronModel([(soma, dend, 15*nS)], v_rest=-60*mV)
model.add_params({
        # spiking parameters
        'Vth': -40*mV,
        'tauw': 150*ms,
        'a': 0*nS,
        'b': 50*pA,
        'Vr': -50*mV,
        # synaptic parameters
        't_AMPA_decay': 2*ms,
        't_NMDA_decay': 60*ms,
        'g_AMPA': 3*nS,
        'g_NMDA': 3*nS})

# Create neuron group
neuron = model.make_neurongroup(
        1, method='euler',
        threshold='V_soma > Vth',
        reset='V_soma = Vr; w += b',
        refractory=4*ms)

# Create necessary Brian objects
poisson = b.PoissonGroup(10, rates=20*Hz)
syn = b.Synapses(poisson, neuron, on_pre='s_AMPA += 1; s_NMDA += 1')
syn.connect(p=1)
mon = b.StateMonitor(neuron, ['V_soma', 'V_dend', 'w'], record=True)
net = b.Network(neuron, poisson, syn, mon)

# Run simulation
b.start_scope()
b.seed(42)
net.run(500*ms)

# Plot results
time = mon.t / ms
v_soma = mon.V_soma[0] / mV
v_dend = mon.V_dend[0] / mV
w = mon.w[0] / pA

fig, (ax1, ax2) = b.subplots(2, 1, figsize=[6, 6], sharex=True)
ax1.plot(time, v_soma, label='soma')
ax1.plot(time, v_dend, label='dendrite', c='C3')
ax1.set_ylabel('Voltage (mV)')
ax1.legend()
ax2.plot(time, w, color='black', label='w')
ax2.set_ylabel('Adaptation current (pA)')
ax2.set_xlabel('Time (ms)')
ax2.legend()
fig.tight_layout()
b.show()

Customizing a NeuronModel

Note that all the customization options available for Soma and Dendrite objects are also available for the final NeuronModel. The methods add_equations, replace_equations, and add_params are shared among these three classes.

print(model.equations)

dV_soma/dt = (gL_soma * (EL_soma-V_soma) + I_soma - w) / C_soma  :volt
I_soma = I_ext_soma + I_dend_soma   :amp
I_ext_soma  :amp
dw/dt = (a * (V_soma-EL_soma) - w) / tauw  :amp
I_dend_soma = (V_dend-V_soma) * g_dend_soma  :amp

dV_dend/dt = (gL_dend * (EL_dend-V_dend) + I_dend) / C_dend  :volt
I_dend = I_ext_dend + I_soma_dend  + I_NMDA + I_AMPA  :amp
I_ext_dend  :amp

I_AMPA = g_AMPA * (E_AMPA-V_dend) * s_AMPA  :amp
ds_AMPA/dt = -s_AMPA / t_AMPA_decay  :1
I_NMDA = g_NMDA * (E_NMDA-V_dend) * s_NMDA / (1 + Mg_con * exp(-Alpha_NMDA*(V_dend/mV+Gamma_NMDA)) / Beta_NMDA)  :amp
ds_NMDA/dt = -s_NMDA/t_NMDA_decay  :1
I_soma_dend = (V_soma-V_dend) * g_soma_dend  :amp

model.replace_equations(
    'I_soma_dend = (V_soma-V_dend) * g_soma_dend  :amp',
    '\n\nHello World!'
)
model.add_equations('Have a beautiful and productive day!')
print(model.equations)

dV_soma/dt = (gL_soma * (EL_soma-V_soma) + I_soma - w) / C_soma  :volt
I_soma = I_ext_soma + I_dend_soma   :amp
I_ext_soma  :amp
dw/dt = (a * (V_soma-EL_soma) - w) / tauw  :amp
I_dend_soma = (V_dend-V_soma) * g_dend_soma  :amp

dV_dend/dt = (gL_dend * (EL_dend-V_dend) + I_dend) / C_dend  :volt
I_dend = I_ext_dend + I_soma_dend  + I_NMDA + I_AMPA  :amp
I_ext_dend  :amp

I_AMPA = g_AMPA * (E_AMPA-V_dend) * s_AMPA  :amp
ds_AMPA/dt = -s_AMPA / t_AMPA_decay  :1
I_NMDA = g_NMDA * (E_NMDA-V_dend) * s_NMDA / (1 + Mg_con * exp(-Alpha_NMDA*(V_dend/mV+Gamma_NMDA)) / Beta_NMDA)  :amp
ds_NMDA/dt = -s_NMDA/t_NMDA_decay  :1


Hello World!

Have a beautiful and productive day!