Structured equations with equation templates

[mstimberg]

Continuing the discussion from #1051 and the Discourse forum here. No need to look at the actual code yet, it is just a quick&dirty dummy implementation to get things working.

As a recap, we started out in #1051 to add a special syntax to add prefixes/suffixes to variable names in equations. The idea was that you could reuse a generic Equations object, e.g. for an ion channel or a synaptic conductance, and then generate various specializations from it. But this turned out to be a bit of a dead-end, since complex examples rather seemed to need a more flexible templating mechanism. We somewhat settled on using a template syntax using {...} which would be very similar to Python's string formatting (and not necessarily add much to it). Fast-forward to 2020 and the discussion on the forum, where @yzerlaut presented his own code to generate equations for complex equations in a structured way, based on classes and string-formatting (with generic equations using placeholders with {...} syntax as well). So I revisited this issue, and I think it could really be a nice addition to Brian. At the first sight, it does not seem to offer much over Python's builtin string formatting syntax, but I think there are a number of potential advantages. Most importantly, we can check the basic syntax already in a template which should make errors much easier to find (think of complex nested expressions with many parentheses).

So here is what my current proposal looks like. It is based on two new classes EquationTemplate and ExpressionTemplate, even though I am not 100% sure that we really should introduce these new classes. An alternative would be that an Equations /Expression object simply is a template as long as it has any {...} placeholders in it. In contrast to our discussion in #1050, I'm not using a substitute or replace method here, but make the objects callable, which leads to a bit more compact code which I think is quite nice.

So here's the example I used for testing. Note that the code is overall longer with the new method, but I think many users would appreciate the more structured representation. The original equations from the Rothman&Manis (2013) example:

# Classical Na channel
eqs_na = """
ina = gnabar*m**3*h*(ENa-v) : amp
dm/dt=q10*(minf-m)/mtau : 1
dh/dt=q10*(hinf-h)/htau : 1
minf = 1./(1+exp(-(vu + 38.) / 7.)) : 1
hinf = 1./(1+exp((vu + 65.) / 6.)) : 1
mtau =  ((10. / (5*exp((vu+60.) / 18.) + 36.*exp(-(vu+60.) / 25.))) + 0.04)*ms : second
htau =  ((100. / (7*exp((vu+60.) / 11.) + 10.*exp(-(vu+60.) / 25.))) + 0.6)*ms : second
"""

# KHT channel (delayed-rectifier K+)
eqs_kht = """
ikht = gkhtbar*(nf*n**2 + (1-nf)*p)*(EK-v) : amp
dn/dt=q10*(ninf-n)/ntau : 1
dp/dt=q10*(pinf-p)/ptau : 1
ninf =   (1 + exp(-(vu + 15) / 5.))**-0.5 : 1
pinf =  1. / (1 + exp(-(vu + 23) / 6.)) : 1
ntau =  ((100. / (11*exp((vu+60) / 24.) + 21*exp(-(vu+60) / 23.))) + 0.7)*ms : second
ptau = ((100. / (4*exp((vu+60) / 32.) + 5*exp(-(vu+60) / 22.))) + 5)*ms : second
"""

For the structured representation, I first define a few templates:

gating_var = EquationTemplate('''d{name}/dt = q10*({name}__inf - {name})/tau_{name} : 1
                                 {name}__inf = {rate_expression}                    : 1
                                 tau_{name} = {tau_scale}/({forward_rate} + 
                                                           {reverse_rate})
                                              + {tau_base}                          : second''')

pos_sigmoid = ExpressionTemplate('1./(1+exp(-(vu - {midpoint}) / {scale}))')
sqrt_sigmoid = ExpressionTemplate('1./(1+exp(-(vu - {midpoint}) / {scale}))**0.5')
neg_sigmoid = ExpressionTemplate('1./(1+exp((vu - {midpoint}) / {scale}))')
exp_voltage_dep = ExpressionTemplate('{magnitude}*exp((vu-{midpoint})/{scale})')
neg_exp_voltage_dep = ExpressionTemplate('{magnitude}*exp(-(vu-{midpoint})/{scale})')

(Potentially one could also nest the templates and put the exp((vu-{midpoint})/{scale} part into its own template that gets reused in all the ExpressionTemplates). Then we can define the currents and make use of the templates to describe the gating variables:

# Classical Na channel
m = gating_var(name='m',
               rate_expression=pos_sigmoid(midpoint=-38., scale=7.),
               forward_rate=exp_voltage_dep(magnitude=5., midpoint=-60, scale=18.),
               reverse_rate=neg_exp_voltage_dep(magnitude=36., midpoint=-60, scale=25.),
               tau_base=0.04*ms, tau_scale=10*ms)
h = gating_var(name='h',
               rate_expression=neg_sigmoid(midpoint=-65., scale=6.),
               forward_rate=exp_voltage_dep(magnitude=7., midpoint=-60., scale=11.),
               reverse_rate=neg_exp_voltage_dep(magnitude=10., midpoint=-60., scale=25.),
               tau_base=0.6*ms, tau_scale=100*ms)

ina = Equations('ina = gnabar*m**3*h*(ENa-v) : amp') + m + h

# KHT channel (delayed-rectifier K+)
n = gating_var(name='n',
               rate_expression=sqrt_sigmoid(midpoint=-15, scale=5.),
               forward_rate=exp_voltage_dep(magnitude=11., midpoint=-60, scale=24.),
               reverse_rate=neg_exp_voltage_dep(magnitude=21., midpoint=-60, scale=23.),
               tau_base=0.7*ms, tau_scale=100*ms)

p = gating_var(name='p',
               rate_expression=pos_sigmoid(midpoint=-23., scale=6.),
               forward_rate=exp_voltage_dep(magnitude=4., midpoint=-60., scale=32.),
               reverse_rate=neg_exp_voltage_dep(magnitude=5., midpoint=-60., scale=22.),
               tau_base=5*ms, tau_scale=100*ms)

ikht = Equations('ikht = gkhtbar*(nf*n**2 + (1-nf)*p)*(EK-v) : amp') + n + p

Potentially we could also have some syntax to make lines like

'''
...
dv/dt = (ileak + ina + ikht + iklt + ika + ih + ihcno + I)/C : volt
'''

independent of the channels (so that it is easy to remove/add channels), but maybe leave this discussion for another time.

Feedback welcome!

[yzerlaut]

I find it a very elegant way to deal with it: it's both simple and very effective. Readability is greatly improved (e.g. for the parameters and their meaning, compare the two versions^^). Conceptually it really makes sense to emphasize the gating_var object in the model, the EquationTemplate object allow to exactly do that. It's great !

[mstimberg]

Thanks for the feedback :) What do you think about the question of EquationTemplate vs. Equations (and same for expressions)? I.e., we could also decide against introducing new classes and simply have Equations with placeholders (which of course would raise an error if you try to use them in a NeuronGroup).

[yzerlaut]

Mmmhhh... I understand that's not a trivial design choice. I liked the Equations, EquationTemplate and EquationExpression, because that's how we intuitively build the dynamical system. Maybe this issue is more about readibility and presentation than anything else, so the above hierarchy of objects is helpful in my opinion. I would keep it like that. But that would be totally fine to just have Equations as well.

[mstimberg]

I worked a bit more on this, and I think I am now quite happy with the general approach (remember that there's still a lot of work to do to make the code robust, have proper error checking, etc.). I had another though about the issue of "scoping", and I think for now we don't have to do much about it. What I think is important, is that it is possible to change the name of e.g. a channel in a single place. This is now possible with a formulation like this:

# Classical Na channel
ina = Equations('ina = gnabar*{m}**3*{h}*(ENa-v) : amp',
                m=gating_var(name='m',
                             rate_expression=pos_sigmoid(midpoint=-38., scale=7.),
                             forward_rate=exp_voltage_dep(magnitude=5., midpoint=-60, scale=18.),
                             reverse_rate=neg_exp_voltage_dep(magnitude=36., midpoint=-60, scale=25.),
                             tau_base=0.04*ms, tau_scale=10*ms),
                h=gating_var(name='h',
                             rate_expression=neg_sigmoid(midpoint=-65., scale=6.),
                             forward_rate=exp_voltage_dep(magnitude=7., midpoint=-60., scale=11.),
                             reverse_rate=neg_exp_voltage_dep(magnitude=10., midpoint=-60., scale=25.),
                             tau_base=0.6*ms, tau_scale=100*ms))

Here, if you are having a clash with the names m and h, you only have to change the name='...' in a single place, e.g. to m_na. I think this is already a big step, we can think about a more elaborate approach at a later point (the nice thing about the current approach is that it is completely on the level of Equations, i.e. nothing about the variable access in NeuronGroup has to change, for example). Note that the gating_var calls here are complete Equations. What does it mean to insert them into another equation? It will replace the placeholder by the name of the first variable defined in the equations, and will add the equations to the original equations. I find this is quite readable/understandable, and avoids that I have to take care to be consistent in the name that I give in the equations and the name that I use in another equation referring to that equation.

The final addition to the syntax is this:

eqs =Equations("""dv/dt = ({currents} + I)/C : volt
                  vu = v/mV : 1  # unitless v
                  I : amp""",
               currents=[ileak, ina, ikht, ih, iklt, ika, ihcno], voltage='vu')

Here, currents is a list of Equations, which is interpreted in the same way as described above, i.e. using the first variable from the equations as its name to insert into the equations, and adding all the equations. The names are combined using the + operator, since I think currents are the main use case for this technique. So this now gives us finally an easy and elegant way to remove/add a channel to a model, no need to manually change the part of the equations where the currents are added up!

As you might have noticed, I have merged EquationTemplate and Equations (same for ExpressionTemplate and Expression). The main reason is: if we only allow {...} placeholders in EquationTemplate, this would mean that all the above examples would have to be templates, but then in the call of the initializer we are already filling in all the placeholders. So we'd either use an ugly

eqs = EquationTemplate("""dv/dt = ({currents} + I)/C : volt
                  vu = v/mV : 1  # unitless v
                  I : amp""")(currents=[ileak, ina, ikht, ih, iklt, ika, ihcno], voltage='vu')

call, or do some dirty Python magic so that calling

eqs = EquationTemplate("""dv/dt = ({currents} + I)/C : volt
                  vu = v/mV : 1  # unitless v
                  I : amp""", currents=[ileak, ina, ikht, ih, iklt, ika, ihcno], voltage='vu')

would actually create an Equations object...

One problem of my current approach is that we now have two mechanisms to replace variables in equations. Before the changes in this branch, you could already do:

eqs = Equations('dx/dt = -x/tau : 1', x='v')

to replace the x by v. The current branch adds the equivalent syntax

eqs = Equations('d{x}/dt = -{x}/tau : 1', x='v')

Not sure whether we should deprecate the old syntax, it might still be useful.

As always: comments/feedback welcome! I converted the full Rothman&Manis example so that it now uses templates for all currents: https://github.com/brian-team/brian2/blob/equation_template/examples/frompapers/Rothman_Manis_2003_with_templates.py

[rat-h]

@mstimberg Could you please give us an example of how to resolve a name clash? For example, if iklt in Rothman&Manis model has standard m and h instead of w and z?

[mstimberg]

@mstimberg Could you please give us an example of how to resolve a name clash? For example, if iklt in Rothman&Manis model has standard m and h instead of w and z?

There is no mechanism to deal with this automatically at the moment, so you'd have to use a non-ambiguous name yourself:

iklt = Equations('ih = gkltbar*{w}**4*{z}*(EK-v) : amp',
                 w=gating_var(name='m_klt',
                              rate_expression=power_sigmoid(midpoint=-48., scale=6., power=0.25),
                              forward_rate=exp_voltage_dep(magnitude=6., midpoint=-60., scale=6.),
                              reverse_rate=neg_exp_voltage_dep(magnitude=16, midpoint=-60, scale=45.),
                              tau_scale=100*ms, tau_base=1.5*ms),
                 z=gating_var(name='h_klt',
                              rate_expression=shifted_neg_sigmoid(shift=zss, midpoint=-71., scale=10.),
                              forward_rate=exp_voltage_dep(magnitude=1., midpoint=-60., scale=20.),
                              reverse_rate=neg_exp_voltage_dep(magnitude=1., midpoint=-60., scale=8.),
                              tau_scale=1*second, tau_base=50*ms))

What I find important is that in the above, you only had to only change the name in a single place, there is no room for an inconsistency between the ih definition and the definition of the gating variables.

Note that if you want to use the suffix approach more generally, you could also define things slightly different in the first place:

gating_var = Equations('''d{name}_{suffix}/dt  = q10*({name}_{suffix}__inf - {name}_{suffix})/tau_{name}_{suffix} : 1
                          {name}_{suffix}__inf = {rate_expression}                                               : 1
                          tau_{name}_{suffix}  = {tau_scale}/({forward_rate} +
                                                    {reverse_rate})
                                                 + {tau_base}                                                    : second''')

# KLT channel (low threshold K+)
iklt = Equations('i_{suffix} = gkltbar*{m}**4*{h}*(EK-v) : amp',
                 m=gating_var(name='m',
                              rate_expression=power_sigmoid(midpoint=-48., scale=6., power=0.25),
                              forward_rate=exp_voltage_dep(magnitude=6., midpoint=-60., scale=6.),
                              reverse_rate=neg_exp_voltage_dep(magnitude=16, midpoint=-60, scale=45.),
                              tau_scale=100*ms, tau_base=1.5*ms),
                 h=gating_var(name='h',
                              rate_expression=shifted_neg_sigmoid(shift=zss, midpoint=-71., scale=10.),
                              forward_rate=exp_voltage_dep(magnitude=1., midpoint=-60., scale=20.),
                              reverse_rate=neg_exp_voltage_dep(magnitude=1., midpoint=-60., scale=8.),
                              tau_scale=1*second, tau_base=50*ms),
                 suffix='klt')

Here, all variables are suffixed and there is only a single place where you have to define the suffix.

As I mentioned earlier, I do not want to make this suffix (or prefix) approach mandatory for now. Otherwise, we'd have to introduce a new syntax (e.g. how do you initialize m of the iklt current specifically) and there are a papers (like the Rothman&Manis that we are implementing here!) that do not use suffixes but unique variable names like w and z.

I get the whole idea of scoping/modularity and there is certainly potential to do something here. Since the other equations never have to refer to the gating variables etc. of a current, they could be treated differently from other variables. But I'd prefer to leave this for another time, since we'd have to come up with new syntax for recording/intializing/accessing such variables. I guess a dot-syntax could work (e.g. iklt.m), but this would be a major change and I think the proposed template syntax here would already be a big step.

[mstimberg]

I've worked on this a bit more and the basic implementation is now mostly finished. It will now look at the dependencies between substitutions and replace things in the correct order. For example, you could have

eqs = Equations('dv/dt = ((E_{leak_suffix} - v) + {currents})/C : 1', currents=['I_Na', 'I_{leak_suffix}'], leak_suffix='L')

and it would use L both in E_L and I_L.

What's missing is mostly error checking, testing, and documentation.

Assuming everyone is ok with the decision to use Equations instead of EquationTemplate (see reasoning in #1244 (comment)) ), there are still a few open questions.

1. What do we do if we can replace both normal names and placeholders? E.g. currently you could have the following:

>>> eqs = Equations('dx/dt = -x/tau_{x} : volt', x='v')
>>> print(eqs)
dv/dt = -v/tau_v : V

i.e. x would be both replaced directly and as a placeholder {x}. I'm not sure whether this is the right approach or whether we should rather raise an error here (i.e. force the user to write d{x}/dt = -{x}/tau_{x} : volt).

2. How general do we want the substitution mechanism to be? One of the advantages of this proposed approach over simple string manipulation is that we can check errors early, and do not have to wait for the equations to be complete. But if we make the system completely general, then we cannot check any syntax. E.g. imagine we allow equations with placeholders for an operator: x {operator} y, where we then replace operator by + or -. I think I would restrict our replacements to names / parts of names, so that we can make basic syntax checks even when there are still placeholders.

As always, comments welcome!

[rat-h]

Thinking about two questions above, it seems that eqs = Equations('dx/dt = ({x}/mV-x)/tau : 1', x='v') should be OK, just extremely confusing. Brian can(must) flag it with warning, but should allow passing through.
The second question is a good one. I agree that placeholders are the substitution for variables/parameters (names), not for operators or functions. We can leave a more general inline equation formation to Python formatting. I'm concerned that if this mechanism will be too general it can create more mess in code than the current version.

[mstimberg]

Thinking about two questions above, it seems that eqs = Equations('dx/dt = ({x}/mV-x)/tau : 1', x='v') should be OK, just extremely confusing. Brian can(must) flag it with warning, but should allow passing through.

Thanks for your opinion on that question, I agree that a warning is probably the best approach (I implemented this in the current PR).

The second question is a good one. I agree that placeholders are the substitution for variables/parameters (names), not for operators or functions. We can leave a more general inline equation formation to Python formatting. I'm concerned that if this mechanism will be too general it can create more mess in code than the current version.

I have the same concern, and I agree that something like replacing operators would really go too far into simple string replacement territory. I have now implemented the following approach: a string has to be either a name, or a full expression. The latter could be done by wrapping it in Expression(...) as well, but I think not making this mandatory makes some use cases a bit simpler. What I like about it is that it makes it possible to use basic expressions in a flexible way. E.g., with the definition

pos_sigmoid = Expression('1./(1+exp(-({voltage} - {midpoint}) / {scale}))')

you can either use v directly and define everything else with units as well:

pos_sigmoid(midpoint=-60*mV, scale=7*mV, voltage='v')

Or you can use a unitless v by directly giving it as an expression (which in this formalism might be clearer than referring to a unitless vu variable):

pos_sigmoid(midpoint=60., scale=7., voltage='v/mV')

The requirement that strings need to be valid expressions means that we can still check for basic syntax correctness (e.g. matching parentheses), without worrying about weird replacements.

I'm quite happy with the way this is turning out. Would be grateful if some of you could try out the current code a bit and see whether there are still things that are cumbersome to do or confusing. If not, I'd go ahead with more code polishing, testing, and documentation.

[rat-h]

Aha, That is a great point @mstimberg! I use unitless or a mixture of "unitfull" and unitless models quite often. This ability to hide units in sub-expressions may be really handy. (^-^)

I cannot commit to doing beta testing for 3-4 weeks, because my shoulder is broken and I'm quite behind the schedule. However, I have a model of a network with multicompartment neurons. They have gradients of channel densities along the main track of apical dendrite, collaterals, and tufts. I can try to move this model from NEURON to Brian.

Honestly, up to now, my hands instinctively picked NEURON as soon as more than a few channels are needed for the model or as soon as complex morphology with channel distribution should is involved. Hopefully, this is going to change and Brian won't be a niche simulator with clean and simple code only for a subset of neural models.

[thesamovar]

Just trying to get back into this after having not thought about it for a while. There's a lot of things pulling in different directions here, and I find it particularly tricky to work out what's the right thing to do.

I think it's important that we have some feature like this, but I'm a bit concerned about coming up with an implementation that we'll have to support going forward but that isn't sufficient for all the things we'd like it to do and so adds technical debt, so I think it's worth taking a little bit more time to discuss. My other concern is that we avoid adding a syntax that has non-obvious or hidden behaviour.

With those in mind, I'm immediately anxious about two parts of the proposed syntax:

• {current} and currents=[i1,i2,i3] to replace the first by i1+i2+i3 because the + is implicit/hidden
• Equations(..., m=gating_var(...), ...) to replace m with the first defined variable in gating_var because the "first defined" rule is implicit/hidden

The second one could potentially be fixed by having a "primary variable" or some such, where in the gating_var=Equations(...) you somehow declare one of the variables to be the primary variable that will be used when you subsequently write something=gating_var(...). We could declare this either with a new flag or with a new keyword to Equations (flag might be better). Or we could do something like m=gating_var(...).name but I'm not a huge fan of that.

For the first one, I'm not sure what would be a good solution. Possibly a specific function for this, e.g. something like equation_sum('I_tot', i1, i2, i3) which returns s+i1+i2+i3, where the +i1+i2+i3 is just the concatenated equations of the currents, and s is an Equations('I_tot = primary var from i1 + primary var from i2 + primary var from i3 (primary)') or something like that. In other words, it creates a new variable with a given name, sets it equal to the sum of the primary variables passed to it, and sets that as the primary variable for the created Equations. The original syntax {currents} would then be replaced by I_tot (or more generally when you do {x} and x is an Equations, then you replace the {x} with the primary variable of those Equations, and concatenate.

OK I didn't exactly mean to do it but I seem to have come up with a concrete proposal. Not sure what I think of it myself yet, so just throwing that out there for discussion for the moment.

I have another point to discuss, but I think I'll put it in a second post to keep things simple.

[thesamovar]

The other point I wanted to discuss was that it would be nice to keep track of the structured representation in the final equations. We could use this for (a) debugging/error messages, (b) for generating better representations of the equations for the automatic model descriptions package. I wonder if we could do this by somehow storing metadata about where substitutions came from when using the substitutions mechanisms.

[thesamovar]

I'm also wondering if we should just bite the bullet and implement something like the dotted syntax. If we end up wanting to do it anyway, maybe it's better to do it now? I'm not sure I agree with this, but if we end up wanting to do it, and we don't do it now, we'll have to support two systems that might interfere with each other, which could turn into a nightmare. This possibility would be reduced if we don't allow all this as part of Equations, but instead add new classes which have to produce an Equations object in the end (but on the other hand, that's much less elegant and composable).

[thesamovar]

On the specific questions:

I think dx/dt=-x/tau_{x} should probably lead to an error. I don't see any use case for it and it could certainly introduce confusion.

I would agree with only looking at name substitutions, not operator etc.

[mstimberg]

Thanks for the feedback @thesamovar , all very reasonable suggestions. Let me start with the easiest point

I think dx/dt=-x/tau_{x} should probably lead to an error.

What I had in mind were some slightly more complex situations where you mix up things from equations and templates. E.g. you'd have one equation that only uses x, and then another equation that only uses {x}, and these get added together before you substitute it. But this is probably overthinking it, I'd be happy to raise an error. I'd probably already raise it for Equations('dx/dt = -x/tau_{x} : 1'), i.e. before even attempting to substitute x.

it would be nice to keep track of the structured representation in the final equations.

Yes, I agree! I'll have to think a bit more about how to do this in practice, but keeping this information around might certainly come handy.

I'm also wondering if we should just bite the bullet and implement something like the dotted syntax.

IMHO, we shouldn't do this now. The nice thing about the proposed approach is that it is completely on the level of Equations, nothing else changes. I don't necessarily see that the two systems would be incompatible. Even with the "dot syntax", I'd still want equation templates to avoid repetition. What I had in mind was something like the following (the (local: ...) flag is the new addition, but never mind the details of the syntax for now):

gating_var = EquationTemplate('''d{name}/dt = q10*({name}__inf - {name})/tau_{name} : 1 (local: {current})
                                 {name}__inf = {rate_expression}                    : 1(local: {current})
                                 tau_{name} = {tau_scale}/({forward_rate} + 
                                                           {reverse_rate})
                                              + {tau_base}                          : second (local: {current})''')

ina = Equations('{current} = gnabar*{m}**3*{h}*(ENa-v) : amp', m=gating_var(name='m',
                  rate_expression=pos_sigmoid(midpoint=-38., scale=7.),
                  forward_rate=exp_voltage_dep(magnitude=5., midpoint=-60, scale=18.),
                  reverse_rate=neg_exp_voltage_dep(magnitude=36., midpoint=-60, scale=25.),
                  tau_base=0.04*ms, tau_scale=10*ms),
                h=gating_var(name='h',
                  rate_expression=neg_sigmoid(midpoint=-65., scale=6.),
                  forward_rate=exp_voltage_dep(magnitude=7., midpoint=-60., scale=11.),
                  reverse_rate=neg_exp_voltage_dep(magnitude=10., midpoint=-60., scale=25.),
                  tau_base=0.6*ms, tau_scale=100*ms))

The idea would be that the "local" variables would not be accessible from elsewhere in the equations (and they could therefore have the same name), internally they would be renamed to something like _local_ina_m, and from outside you could access them via neurons.ina.m, or record them as ina.m, etc. I feel that would make quite use of the dot syntax without making things too complicated internally (but we'd have take care to make things not too confusing in error messages, etc.).

[mstimberg]

About the questions regarding the summing and the "primary" variable. I also see the danger of hiding away details here, but at the same time we want to be concise and not introduce too much new syntax. I feel there is a danger of ending up with equations that have a lot of idiosyncratic stuff like (primary) in it, and we can no longer claim that our model descriptions are "human-readable mathematical descriptions"...

For the "primary equation" question, I tend towards saying that using the first equation would be ok and consistent with the way you'd write these equations naturally. Adding new syntax in the equations would also come with new problems, e.g. what if you add two equations together that both have a (primary) flag?

For the summing, I agree that it is maybe a bit too magic. I find the equation_sum syntax a bit too wordy though and feel that it would take away from the simplicity of the template syntax. How about embedding this information in the template placeholder itself? E.g. something like {+currents} or {sum(currents)}. Both would also allow for a straightforward extension to a product which I could think as potentially being useful (something like Equations(... {product(gatings)}*g*(E - v)...', gatings=[m, m, m, h])`).

[rat-h]

neurons.dendrite.ina.m is what I think can make notation less messy. However, it isn't clear to me why you want to introduce local: flag? The structure of definitions should be enough to figure it out: all {placeholders} in a template, which are defined in an equation are local to the equation/template (i.e there names _equation_name). All {placeholders} in an equation which is included into another equation and defined there are local to this set of equations. Let's try:

gating_var = EquationTemplate('''d{name}/dt = q10*({name}__inf - {name})/tau_{name} : 1
                                 {name}__inf = {rate_expression}                    : 1
                                 tau_{name} = {tau_scale}/({forward_rate} + {reverse_rate})
                                              + {tau_base}                          : second ''')

ina = Equations('{name} = gnabar*{m}**3*{h}*(ENa-v) : amp/meter**2', 
               name = 'ina',
               include =[
                 gating_var(name='m',
                   rate_expression=pos_sigmoid(midpoint=-38., scale=7.),
                   forward_rate=exp_voltage_dep(magnitude=5., midpoint=-60, scale=18.),
                   reverse_rate=neg_exp_voltage_dep(magnitude=36., midpoint=-60, scale=25.),
                   tau_base=0.04*ms),
                gating_var(name='h',
                  rate_expression=neg_sigmoid(midpoint=-65., scale=6.),
                  forward_rate=exp_voltage_dep(magnitude=7., midpoint=-60., scale=11.),
                  reverse_rate=neg_exp_voltage_dep(magnitude=10., midpoint=-60., scale=25.),
                  tau_base=0.6*ms, tau_scale=100*ms)
              ]
      )

i =  Equations('{name} = ina+ik+ileak: amp/meter**2', 
               name = 'totcur',
               include=[ina,ik,ileak]
)

neuron = SpatialNeuron(morphology=morpho,
         model="""
                  {totalcurrent}
                  totcur.ina.m.tau_scale : second 
                 """, 
         Cm=1*uF/cm**2, Ri=100*ohm*cm,
          totalcurrent = i
)
neuron.totcur.ina.m.tau_scale = 1*ms
neuron.basal.totcur.ina.m.tau_scale = 10*ms
neuron.apical[10:].totcur.ina,m.tau_scale = 3*ms

It may look a bit too much, but as soon as we depart from single compartment model, like Rothman&Manis, structural approach is really needed.

[mstimberg]

I see your point, it would indeed be a possibility to make local variables the default when inserting an equation into another equation. I feel that this would limit things quite a bit to this specific way of constructing equations, though. For example, you could not concatenate together two equations that both are supposed to have local variables. Maybe this would be ok, not sure. But as I said in my comment to @thesamovar, I feel that we should improve things more incrementally. For me this change feels too big. Even if you have tens of different channels, giving them each a suffix is quite a straightforward solution and also what you'd do to describe the model in a paper, for example.

I'm not too happy about your tocur.ina.m.tau_scale proposal in the final equations. I guess you want to use this so this becomes a per-compartment parameter instead of a constant, but this treats a placeholder {tau_scale} as a variable name which I feel makes things messy.

It may look a bit too much, but as soon as we depart from single compartment model, like Rothman&Manis, structural approach is really needed.

I don't quite see why multi-compartment vs. single-compartment makes a difference here, to be honest. The equations will not look any different.

[rat-h]

I don't quite see why multi-compartment vs. single-compartment makes a difference here, to be honest. The equations will not look any different.

Because we may have different subunits of ion channel in basal and apical dendrites (for example) and it is necessary to assess parameters of equations and set different taus or steady-state voltages /calcium consecrations in different parts of a neuron. Classical example is K_A in proximal and apical dendrites of CA1 pyramidal neuron (see Jarsky 2005 for example).

I'm not too happy about your tocur.ina.m.tau_scale proposal in the final equations.

I tried to make notation homogeneous. If placeholder isn't found on local level it is recursively sent up with automatically generates dot prefix until it will be define. But this may not be a good way to do so. Another version is

ina = Equations(
             '''{name} = gnabar*{m}**3*{h}*(ENa-v) : amp/meter**2
               m.tau_scale : seconds
              ''', 
               name = 'ina',
               include =[
                 gating_var(name='m',
                   rate_expression=pos_sigmoid(midpoint=-38., scale=7.),
                   forward_rate=exp_voltage_dep(magnitude=5., midpoint=-60, scale=18.),
                   reverse_rate=neg_exp_voltage_dep(magnitude=36., midpoint=-60, scale=25.),
                   tau_base=0.04*ms),
                gating_var(name='h',
                  rate_expression=neg_sigmoid(midpoint=-65., scale=6.),
                  forward_rate=exp_voltage_dep(magnitude=7., midpoint=-60., scale=11.),
                  reverse_rate=neg_exp_voltage_dep(magnitude=10., midpoint=-60., scale=25.),
                  tau_base=0.6*ms, tau_scale=100*ms)
              ]
      )

[rat-h]

I feel that this would limit things quite a bit to this specific way of constructing equations, though.

Maybe in this case you can introduce a flag global if something shouldn't be in hierarchical tree.

For example, you could not concatenate together two equations that both are supposed to have local variables.

Yes should be an error only if both equations have the same local variable. However, one can concatenate two equations if variables are different, right? It will be just a new equation with all local variables from both:

ihh = ina+ik
ihh.m = ...
ihh.h = ...
ihh.n = ...

[mstimberg]

Because we may have different subunits of ion channel in basal and apical dendrites (for example) and it is necessary to assess parameters of equations and set different taus or steady-state voltages /calcium consecrations in different parts of a neuron.

Right, I get that. But this simply means that you have things that are parameters in the equations, not constants, and again you can handle this with suffixes, say calling it tau_scale_m_Na. This is the same for single-compartment models if things vary across neurons.

I still find that your proposed notation mixes up the placeholder {tau_scale} with the name tau_scale. My preference (at least for now), would still be to write it as something like this, which is already supported by the current PR:

ina = Equations('i_{suffix} = gnabar*{m}**3*{h}*(E_{suffix}-v) : amp/meter**2',
                suffix='na',
                m=gating_var(name='m_{suffix}',
                             rate_expression=pos_sigmoid(midpoint=-38., scale=7.),
                             forward_rate=exp_voltage_dep(magnitude=5., midpoint=-60, scale=18.),
                             reverse_rate=neg_exp_voltage_dep(magnitude=36., midpoint=-60, scale=25.),
                             tau_base=0.04*ms, tau_scale=Equations('tau_scale_{name} : second (constant)')),
                h=gating_var(name='h_{suffix}',
                             rate_expression=neg_sigmoid(midpoint=-65., scale=6.),
                             forward_rate=exp_voltage_dep(magnitude=7., midpoint=-60., scale=11.),
                             reverse_rate=neg_exp_voltage_dep(magnitude=10., midpoint=-60., scale=25.),
                             tau_base=0.6*ms, tau_scale=100*ms)
                )

Apologies for being impressed by my own code, but look at the beauty of tau_scale_{name} → tau_scale_m_{suffix} → tau_scale_m_na and how not a single variable name or suffix has to be specified more than once and how natural having tau_scale as a parameter in m_na but a constant in h_na is 😄 !

[mstimberg]

(Sorry, I wrote my reply in parallel to your second post. Will sleep over things and get back to you tomorrow!)

[rat-h]

I like your tau_scale=Equations('tau_scale_{name} : second (constant)') - brilliant!

IMHO, prefixes are better than suffixes. As in any unix system, /var/log/system.log looks better than system.log/log/var/ . Also in dot notation A.B.C it is more intuitively that B is a part of A and C is a part of B. It is clear that placeholder can do both name="{prefix}_m".

However, here why I'm a bit concern of just free-handed placeholders: look at your perfect code, where is the suffix defined and where it is used? The same for the {name} placeholder. It seems legit to do this definition:

ina = Equations('i_{suffix} = gnabar*{m}**3*{h}*(E_{suffix}-v) : amp/meter**2',
                suffix='na',
                name='ina',
                m=gating_var(name='m_{suffix}', suffix='xxx',
                             rate_expression=pos_sigmoid(midpoint=-38., scale=7.),
                             forward_rate=exp_voltage_dep(magnitude=5., midpoint=-60, scale=18.),
                             reverse_rate=neg_exp_voltage_dep(magnitude=36., midpoint=-60, scale=25.),
                             tau_base=0.04*ms, tau_scale=Equations('tau_scale_{name} : second (constant)')),
                h=gating_var(name='h_{suffix}',
                             rate_expression=neg_sigmoid(midpoint=-65., scale=6.),
                             forward_rate=exp_voltage_dep(magnitude=7., midpoint=-60., scale=11.),
                             reverse_rate=neg_exp_voltage_dep(magnitude=10., midpoint=-60., scale=25.),
                             tau_base=0.6*ms, tau_scale=100*ms)
                )

or even worse name="{name}_m". I agree, that dot notation may be quite restrictive, but I think what you are proposing is too flexible. There should be an optimal looseness 🤪

[rat-h]

What isn't grate about your perfect code is that there is m and there is m_{suffix} here

m=gating_var(name='m_{suffix}'

and it is easy to mess what is there.

In contrast, in my naive suggestion, everything is straightforward: each Equations() has name and, if needed, include and all of this creates hierarchical tree. It doesn't need to hide names, because their locality is defined by how you define them. User needs to put an additional flag gloabal if a placeholder needs to be out of tree. With your brilliant idea, all placeholders must be substituted at the moment when the object is created creation ... so my proposal should look like this:

gating_var = EquationTemplate('''d{name}/dt = q10*({name}__inf - {name})/tau_{name} : 1
                                 {name}__inf = {rate_expression}                    : 1
                                 tau_{name} = {tau_scale}/({forward_rate} + {reverse_rate})
                                              + {tau_base}                          : second ''')

xxx = Equations('{name} = {gbar}*{m}**3*{h}*(ENa-v) : amp/meter**2', 
               name = 'ina',
               include =[
                 gating_var(name='m',
                   rate_expression=pos_sigmoid(midpoint=-38., scale=7.),
                   forward_rate=exp_voltage_dep(magnitude=5., midpoint=-60, scale=18.),
                   reverse_rate=neg_exp_voltage_dep(magnitude=36., midpoint=-60, scale=25.),
                   tau_base=0.04*ms, tau_scale=Equations('{name} : second ', name="tau_scale")),
                 gating_var(name='h',
                   rate_expression=neg_sigmoid(midpoint=-65., scale=6.),
                   forward_rate=exp_voltage_dep(magnitude=7., midpoint=-60., scale=11.),
                   reverse_rate=neg_exp_voltage_dep(magnitude=10., midpoint=-60., scale=25.),
                   tau_base=0.6*ms, tau_scale=100*ms),
                Equations('{name} : semens/meter**2', name="gbar"))
              ]
      )

i =  Equations('{name} = ina+ik+ileak: amp/meter**2', 
               name = 'totcur',
               include=[xxx,ik,ileak]
)

neuron = SpatialNeuron(morphology=morpho,
         model="{totalcurrent}", 
         Cm=1*uF/cm**2, Ri=100*ohm*cm,
         totalcurrent = i
)
neuron.totcur.ina.m.tau_scale = 1*ms
neuron.basal.totcur.ina.m.tau_scale = 10*ms
neuron.apical[10:].totcur.ina.m.tau_scale = 3*ms
neuron.axon.totcur.ina.gbar = 5000*nS/cm**2

(Sorry for putting comments, while you are asleep, now is my turn ;)

[mstimberg]

There is really a fundamental question here how radical we want the current system to change: at the moment, equations are strings, and all my proposed system does is help you create this string. There is no change needed to anything else in Brian, and you can easily inspect the final equations.
In your (@rat-h) proposal, equations mor closely reflect the hierarchical approach of NEURON, etc. and are represented by a more complex object. There is no straightforward string representation of the final equations any more. Or rather, we could still have a fairly straightforward string representation like dina.m/dt = (ina.m_inf - ina.m)/ina.tau_m : 1, but this would not correspond to the actual code we run (which would rather use something like _ina_m) and is slightly confusing as mathematical notation since some people use a simple dot to mean multiplication. Just for displaying it we could of course have work arounds (showing some block structure with indenting, etc.), but e.g. dumping the full equations to disk in a way that can be reused later would need a more complex serialization than a flat string (and therefore additional work).

What isn't grate about your perfect code is that there is m and there is m_{suffix} here
m=gating_var(name='m_{suffix}'

I don't think I agree with this as a problem: the first m is the name {m} from the template, it is a placeholder and not linked to the actual name used in equations. If you were concerned about this, you could use something like '{name} = {gbar}*{gating1}**3*{gating2}*(ENa-v) and then gating1=gating_var(name='m_{suffix}'). What I find important that if you wanted to change the name of the gating variable from m_{suffix} to something else like simply m or p or whatever, there is only a single place where you have to do it (m=gating_var(name='p', ...) might be confusing at first sight but there is nothing ambiguous about it).

Regarding your example with multiple names and suffixes all over the place (which currently would not work, but it might actually make sense to remove the limitation that makes it fail), I agree that this would be confusing. But any reasonable complex system will make it possible for the user to mess things up. I hope you don't take this as defensive/returning the blame (I am sincerely grateful for this discussion!), but there are also things that I find confusing about your proposed system: there are two different ways to know what fills in a placeholder: a placeholder can either be filled directly from an argument as in Equations('{name} = ...', name = 'ina'), or indirectly like {m} which gets filled because somewhere in includes a template has the attribute name='m'. So the link goes the other way round, instead of specifying what m is, we search for things named m. The dot notation also potentially gives rise to some redundancy e.g. in names like ina.m.m_inf.

Sorry for jumping around between topics, but one more comment:

IMHO, prefixes are better than suffixes. As in any unix system, /var/log/system.log looks better than system.log/log/var/ . Also in dot notation A.B.C it is more intuitively that B is a part of A and C is a part of B.

I agree with this in general and from the computer science/programming language perspective. But one of the guiding principles of Brian's equation mechanism has always been that we strive to make the link between equations in a paper and the equations you specify in Brian as clear as possible. Equations in a paper traditionally use suffixes in variable names like tau_m, alpha_m, etc. instead of, say, m.tau, m.alpha, even though the latter could be considered "better" according to your reasoning.

Taking a step back and looking back at your proposal, I wonder about a general change that gives us more flexibility in the long-term: currently, Equations can have replacements as keyword arguments to the initializer directly, and my proposal extends this usage to replacements for placeholders. Maybe it would indeed be cleaner to discontinue this (at the risk of some backwards-incompatibility) and only have standard arguments at this place. So things like Equations(..., voltage='...', suffix='...') would no longer be allowed, and instead you'd have to use something like

eqs = Equations('d{x}/dt = -{x}/tau_{x}', replacements={'x': 'v'})

or

eqs = Equations('d{x}/dt = -{x}/tau_{x}')
v_eqs = eqs(x='v')

This would free up the Equations initializer for standardized keywords like your name and include, but also for things like primary or description/comment/metadata/reference or whatnot. That would mean breaking backwards-compatibility for the current Equations substitution syntax, but we could deal with this with warnings mostly.

My current proposal would look less clean with this unfortunately, so you could no longer do:

ina = Equations('ina = gnabar*{m}**3*{h}*(ENa-v) : amp',
                m=gating_var(name='m',
                             rate_expression=pos_sigmoid(midpoint=-38., scale=7.),
                             forward_rate=exp_voltage_dep(magnitude=5., midpoint=-60, scale=18.),
                             reverse_rate=neg_exp_voltage_dep(magnitude=36., midpoint=-60, scale=25.),
                             tau_base=0.04*ms, tau_scale=10*ms),
                h=...)

but would have to do either something like

 ina = Equations('ina = gnabar*{m}**3*{h}*(ENa-v) : amp',
                replacements={'m': gating_var(name='m',
                             rate_expression=pos_sigmoid(midpoint=-38., scale=7.),
                             forward_rate=exp_voltage_dep(magnitude=5., midpoint=-60, scale=18.),
                             reverse_rate=neg_exp_voltage_dep(magnitude=36., midpoint=-60, scale=25.),
                             tau_base=0.04*ms, tau_scale=10*ms),
                                        'h':...)

or

```Python
 ina = Equations('ina = gnabar*{m}**3*{h}*(ENa-v) : amp')
 ina = ina(m=gating_var(name='m',
                             rate_expression=pos_sigmoid(midpoint=-38., scale=7.),
                             forward_rate=exp_voltage_dep(magnitude=5., midpoint=-60, scale=18.),
                             reverse_rate=neg_exp_voltage_dep(magnitude=36., midpoint=-60, scale=25.),
                             tau_base=0.04*ms, tau_scale=10*ms),
                                       h=...)

which then renders my previous argument over Equations vs. EquationsTemplate somewhat invalid 😬

Ok, I need to think about this a bit more and put some more examples for the various syntax proposals together. Let's ignore the currents/equation_sum question for now and focus on the insertion of equations into equations and the related topics.

[thesamovar]

I'd probably already raise it for Equations('dx/dt = -x/tau_{x} : 1'), i.e. before even attempting to substitute x.

Agreed!

For the "primary equation" question, I tend towards saying that using the first equation would be ok and consistent with the way you'd write these equations naturally. Adding new syntax in the equations would also come with new problems, e.g. what if you add two equations together that both have a (primary) flag?

This precisely highlights the problem with just using the first defined variable to me. If there is an implicit rule (first variable) then if you write eq1+eq2 it will mean something different than eq2+eq1 which in my view is very unexpected behaviour! But I agree that we need to ensure we come up with something human readable. I think I'm happy to delay further discussion of this until you do your concrete syntax proposals post.

{sum(currents)}

I like this!

Overall, I think I feel more convinced by Marcel's argument for making a new system that fits with the current system rather than introducing dotted notation. We do want to make the correspondence between how you'd write it in a paper and the code as close as possible. It seems to me that there are two ways you might write it in a paper, either writing out all the equations in one huge set, with unique names for all variables, or as a seqeuence of separate sets of equations with explicit points of coupling between them. For large sets of equations, this latter approach is probably more readable and more or less corresponds to Marcel's construction. If we have metadata that shows how the equations were originally structured, then this correspondence could be maintained and pretty explicit I feel. Introducing dotted names breaks this connection somewhat.

[rat-h]

I hope you don't take this as defensive/returning the blame

@mstimberg, you don't need to say this! It is nice and hopefully useful discussion between colleagues, good Brianstorm 😄

In your (@rat-h) proposal, equations mor closely reflect the hierarchical approach of NEURON, etc.

I bet you are incorrect! In the original hoc language, all gating variables and parameters in equations had module suffix dend(0.5) m_hh = 0.5 (actually SUFFIX field in NEURON structure is still in mod format). It was quite a good way for making bugs in the code, because names were all over the place. I was in the room of NEURON course at SfN 2014 (it was in DC) when Michael Hines admitted that suffix notation wasn't a good idea, allowing user confusion. I was very happy when they switch to dot notation in python. It helps to avoid mistakes and bugs in complex models. I just afraid that Brian is going to run into the same wall!
Now you know my motivation 😆

Returning to your proposal...

1. please show how to change the tau_scale of an m-gating variable in the sodium channel along the dendrite. It isn't clear to me how placeholder {prefix}_tau_scale which should be substituted at the object creation (as far as I understand, please correct me if I'm wrong), will be converted into a parameter, available for modifications in different compartments.
2. it feels that placeholders and their replacements should have some kind of locality. Like in C or python, always use local replacements (defined in current object) and if there is not a local replacement for a placeholder - request a parent class.

[mstimberg]

Good Brianstorm indeed 😁

I bet you are incorrect! In the original hoc language, all gating variables and parameters in equations had module suffix

You are of course right – it has been a while since I wrote anything in hoc... But I wasn't so much focussed on the name syntax used in NEURON, but rather the fact that you typically define each channel/mechanism in a separate mod file. So while you can set all the (exposed) parameters from the main hoc script, all the mechanisms stay separated whereas Brian works with a single big equation in the end, i.e. something that would be more comparable to a merged mod file.

1. please show how to change the tau_scale of an m-gating variable in the sodium channel along the dendrite. It isn't clear to me how placeholder `{prefix}_tau_scale` which should be substituted at the object creation (as far as I understand, please correct me if I'm wrong), will be converted into a parameter, available for modifications in different compartments.

This was what my earlier example did (using a suffix instead of a prefix):

ina = Equations('i_{suffix} = gnabar*{m}**3*{h}*(E_{suffix}-v) : amp/meter**2',
                suffix='na',
                m=gating_var(name='m_{suffix}',
                             rate_expression=pos_sigmoid(midpoint=-38., scale=7.),
                             forward_rate=exp_voltage_dep(magnitude=5., midpoint=-60, scale=18.),
                             reverse_rate=neg_exp_voltage_dep(magnitude=36., midpoint=-60, scale=25.),
                             tau_base=0.04*ms, tau_scale=Equations('tau_scale_{name} : second (constant)')),
                h=gating_var(name='h_{suffix}',
                             rate_expression=neg_sigmoid(midpoint=-65., scale=6.),
                             forward_rate=exp_voltage_dep(magnitude=7., midpoint=-60., scale=11.),
                             reverse_rate=neg_exp_voltage_dep(magnitude=10., midpoint=-60., scale=25.),
                             tau_base=0.6*ms, tau_scale=100*ms)
                )

With the additional voltage='v/mV' replacement, this leads to these equations:

h_na__inf = (1./(1+exp(((v/mV) - (-65.0)) / 6.0))) : 1
m_na__inf = (1./(1+exp(-((v/mV) - (-38.0)) / 7.0))) : 1
tau_m_na = tau_scale_m_na/((5.0*exp(((v/mV)-(-60))/18.0)) + (36.0*exp(-((v/mV)-(-60))/25.0))) + (40. * usecond) : s
i_na = gnabar*m_na**3*h_na*(E_na-v) : A/m^2
tau_h_na = (100. * msecond)/((7.0*exp(((v/mV)-(-60.0))/11.0)) + (10.0*exp(-((v/mV)-(-60.0))/25.0))) + (0.6 * msecond) : s
dh_na/dt = q10*(h_na__inf - h_na)/tau_h_na : 1
dm_na/dt = q10*(m_na__inf - m_na)/tau_m_na : 1
tau_scale_m_na : s (constant)

Here, tau_scale_m_na is now a parameter that can vary across neurons or compartments.

I did not yet get around preparing an overview over the different syntax proposals but hopefully I'll finish this soon.

[rat-h]

It seems you are proposing "flat" placeholder replacement. There is something confusing in your example (but it may be just me 0_0): in the line

 m=gating_var(name='m_{suffix}',

you are creating an object of gating_var, which has an argument name, which requires a replacement for {suffix}, which in turn will be defined much later as an argument for Equations. So this replacement should wait until Equations will be created and grab from where {suffix}, right? Now consider two channels, with the same m and h but different suffixes:

ina = Equations('i_{suffix} = gnabar*{m}**3*{h}*(E_{suffix}-v) : amp/meter**2',
                suffix='na',
                m=gating_var(name='m_{suffix}', .... )
                h=gating_var(name='h_{suffix}', ... )
                )
ical = Equations('i_{suffix} = gcal*{m}**3*{h}*(E_{suffix}-v) : amp/meter**2',
                suffix='cal',
                m=gating_var(name='m_{suffix}', .... )
                h=gating_var(name='h_{suffix}', ... )
                )

I'm sure, you have a solution :), but if you don't, I would rather require that all placeholders must be substituted by replacements in the arguments unless the flag (parent) is set before the placeholder, like this:

ina = Equations('i_{suffix} = gnabar*{m}**3*{h}*(E_{suffix}-v) : amp/meter**2',
                suffix='na',
                m=gating_var(name='m_{(parent):suffix}', .... )
                h=gating_var(name='h_{(parent):suffix}', ... )
                )
ical = Equations('i_{suffix} = gcal*{m}**3*{h}*(E_{suffix}-v) : amp/meter**2',
                suffix='cal',
                m=gating_var(name='m_{(parent):suffix}', .... )
                h=gating_var(name='h_{(parent):suffix}', ... )
                )

Well, don't kill me please, now we can set a placeholder to the replacement which will be defined in a specific object, for example name="{ical:suffix}"! You can see a use of this notation for example in summing all calcium currents to compute inner calcium concentration, or to access calcium concentrations inside/outside to compute GHK current.

[mstimberg]

I think you are overcomplicating things a bit here, there is nothing problematic in your first example:

ina = Equations('i_{suffix} = gnabar*{m}**3*{h}*(E_{suffix}-v) : amp/meter**2',
                suffix='na',
                m=gating_var(name='m_{suffix}', .... )
                h=gating_var(name='h_{suffix}', ... )
                )

Python will first evaluate gating_var(name='m_{suffix}', ...) which will return a new Equations object like this:

m_{suffix}__inf = (1./(1+exp(-({voltage} - (-38.0)) / 7.0))) : 1
tau_m_{suffix} = (100. * msecond)/((5.0*exp(({voltage}-(-60))/18.0)) + (36.0*exp(-({voltage}-(-60))/25.0))) + (40. * usecond) : s
dm_{suffix}/dt = q10*(m_{suffix}__inf - m_{suffix})/tau_m_{suffix} : 1

So m_{suffix} has been used as the name everywhere. The same happens for the h gating variable. Now, these two Equations objects are provided as arguments to the Equations object defining ina. It sees that the gating variables refer to {suffix}, so it will first insert the equations to lead to something like this:

i_{suffix} = gnabar*m_{suffix}**3*h_{suffix}*(E_{suffix}-v) : A/m^2
dm_{suffix}/dt = q10*(m_{suffix}__inf - m_{suffix})/tau_m_{suffix} : 1
tau_m_{suffix} = (100. * msecond)/((5.0*exp(({voltage}-(-60))/18.0)) + (36.0*exp(-({voltage}-(-60))/25.0))) + (40. * usecond) : s
m_{suffix}__inf = (1./(1+exp(-({voltage} - (-38.0)) / 7.0))) : 1
... # h

And then it will replace {suffix} by na:

i_na = gnabar*m_na**3*h_na*(E_na-v) : A/m^2
dm_na/dt = q10*(m_na__inf - m_na)/tau_m_na : 1
tau_m_na = (100. * msecond)/((5.0*exp(({voltage}-(-60))/18.0)) + (36.0*exp(-({voltage}-(-60))/25.0))) + (40. * usecond) : s
m_na__inf = (1./(1+exp(-({voltage} - (-38.0)) / 7.0))) : 1
# ...

So ina does not contain any reference to {suffix} any more. That ical then also uses {suffix} isn't a problem.

If you don't believe me you can try it out in this PR branch :)

[rat-h]

I got it.

Now can we figure out how to write in new notation for (say) 3 Ca2+ ion channels CaL, CaH, and CaT, The equations for these channels are GHK and they use inner calcium consecration, which is computed in another equation. The later one has to sum all calcium currents + subportion of NMDA current to compute inner concentration? I wonder if placeholder/replacements will make this system more elegant and simpler.

[mstimberg]

It would be definitely useful to have a look at something like this. Do you have a reference that uses the kind of equations you have in mind?

[rat-h]

Here an example Multitarget Multiscale Simulation for Pharmacological Treatment of Dystonia in Motor Cortex with ModelDB record. However, it is a very good but heavy model. I'll try to extract equations from there (will take a bit of time) and put them here, so you can work with equations, not a heavy NEURON model.

[mstimberg]

That would be great, thanks! No need to hurry.

[rat-h]

n-calcium channel

gcan = gcanbar*m*m*h* ki/(ki+cai)
KTF = ((25./293.15)*(celsius + 273.15))
f = KTF/2
nu = v/f
efun = lambda z: 1 - z/2 if fabs(z) < 1e-4 else z/(exp(z) - 1)
ghk=-f*(1. - (ci/co)*exp(nu))*efun(nu)
ica = gcan*ghk

alph = 1.6e-4*exp(-v/48.4)
beth = 1/(exp((-v+39.0)/10.)+1.)
alpm = 0.1967*(-1.0*v+19.88)/(exp((-1.0*v+19.88)/10.0)-1.0)
betm = 0.046*exp(-v/20.73)
alpmt = exp(0.0378*zetam*(v-vhalfm)) 
betmt = exp(0.0378*zetam*gmm*(v-vhalfm)) 
qt=q10^((celsius-25)/10)
a = alpm(v)
b = 1/(a + betm(v))
minf = a*b
taum = betmt(v)/(qt*a0m*(1+alpmt(v)))
if (taum<mmin/qt) {taum=mmin/qt}
a = alph(v)
b = 1/(a + beth(v))
hinf = a*b
tauh= 80

dm/dt = (minf - m)/taum
dh/dt = (hinf - h)/tauh

l-calcium channel

gcal = gcalbar*m*m*ki/(ki+cai)

KTF = ((25./293.15)*(celsius + 273.15))
f = KTF/2
nu = v/f
ghk=-f*(1. - (ci/co)*exp(nu))*efun(nu)

ica = gcal*ghk

alp = 15.69*(-1.0*v+81.5)/(exp((-1.0*v+81.5)/10.0)-1.0)
bet = 0.29*exp(-v/10.86)
alpmt = exp(0.0378*zetam*(v-vhalfm)) 
betmt = exp(0.0378*zetam*gmm*(v-vhalfm)) 
dm/dt = (minf - m)/tau
qt=q10^((celsius-25)/10)
a = alp(v)
b = 1/((a + bet(v)))
minf = a*b
tau = betmt(v)/(qt*a0m*(1+alpmt(v)))
if (tau<mmin/qt) {tau=mmin/qt}

similar equations for t-calcium channel

Calcium influx through NMDA and AMPA receptors

mgblock = 1.0 / (1.0 + 0.28 * exp(-0.062(/mV) * v) )
ampai =  ampag * (v - e) *           (1-AMPAfracca)
nmdai =  nmdag * (v - e) * mgblock * (1-NMDAfracca)
i_ca2p = 0
if(AMPAfracca>0.0){ : calcium current through AMPA
    i_ca2p = i_ca2p +  ampag * ghkg(v,cai,cao,2) *           AMPAfracca
}
if(NMDAfracca>0.0){ : calcium current through NMDA
    i_ca2p = i_ca2p +  nmdag * ghkg(v,cai,cao,2) * mgblock * NMDAfracca
}
ica = i_ca2p
i = ampai + nmdai

Calcium dynamics in simplest way

ica = ica_n+ica_t+ica_l+ica_syn
cai' = -(1000) * alpha * ica - cai/tau 

but in a more complex model calcium have to be able to diffuse into neighboring compartments and deep into dendrite. There is a problem here. When we accurately model cable properties of a dendrite, we usually end up with pretty small segments and therefore, if a segment has a synapse, it will be overflown by calcium. To avoid this, longitudinal diffusion is needed.

LONGITUDINAL_DIFFUSION DCa*PI*(diam*depth-depth*depth) {cai}
~ cai << (-ica*PI*diam/(2*FARADAY)- (cai - minCai)/decay) : ica is Ca infflux	, and simple first order pumping
~ cai + Buffer <-> CaBuffer (k1buf*PI*(diam*depth-depth*depth), k2buf*PI*(diam*depth-depth*depth)) :Calmodulin bufferening RxD

[mstimberg]

Thanks for the models, I will have a closer look soon.

but in a more complex model calcium have to be able to diffuse into neighboring compartments and deep into dendrite.

This is a general limitation in Brian at the moment, the only thing that can "diffuse" is the membrane potential. There are ways to work around the limitations but not with the same numerical accuracy as for the longitudinal currents. In principle, it should be possible to reuse the same mechanism, but it's some non-trivial work and our human resources to work on this are quite limited unfortunately...

[rat-h]

This is a general limitation in Brian at the moment, the only thing that can "diffuse" is the membrane potential.

Wait, I thought we have access to the dynamical variables in the nearby segments! Am I missing something?

If we don't have access to dynamical variables, i.e. one cannot write a diffusion implementation, I would suggest to return to @thesamovar post and "bite the bullet". It seems you, guys can apply for a grant to make Brian universal tools, not a niche product. Brian is very popular software and good tool, so make it better is an excellent point for a proposal.

[mstimberg]

Wait, I thought we have access to the dynamical variables in the nearby segments! Am I missing something?

A segment is like a neuron, it does not have access to variables outside itself. The only way to get access is to connect it via Synapses (which would be the ugly workaround to implement diffusion). Or did I misunderstand what you meant?

If we don't have access to dynamical variables, i.e. one cannot write a diffusion implementation, I would suggest to return to @thesamovar post and "bite the bullet". It seems you, guys can apply for a grant to make Brian universal tools, not a niche product. Brian is very popular software and good tool, so make it better is an excellent point for a proposal.

I don't quite get the reference to Dan's post, but we'd certainly be happy to get money and people to work on various aspects of Brian. We know from experience that getting money for this kind of work is incredibly hard, almost all related grants require you to propose new software that does something that wasn't possible before... Usability improvements, maintenance, etc. is almost impossible to fund via grants. I've seen some NSF grants in that direction, but these are of course US-only.