Let’s say that I have a function which I want to use for synaptic current (I_syn) in a Hodgkin-Huxley model. Roughly, something like this:
dv/dt = 1/c * (sum(I_syn) + I_Na +...)
This current consists of two parts, one of which depends on the presynaptic spike time (t_spike)
I_syn(t) = A * ( exp(-(t-t_spike)/tau1) - exp(-(t-t_spike)/tau2) ) * (v - E_syn)
and let’s say I_syn lasts only for a certain amount of time (T).
I tried to implement this in brian2 and I think that I need the variable lastupdate for this. However, I’m not sure how to utilize it.
Also, on a different note, how am I supposed to include the currents coming from all the synapses? I mean, during the simulation I have the number of the synapses entering each neuron but I don’t know how to create a function that iterates over all the synapses and sums the currents.
I would really appreciate it if you could help me with this issue.
Hi @LBENE . Synapses are often described by equations like the ones you post, but these equations show some practical difficulties: you have to sum up the currents over all recent spikes, introduce some cutoff value when to forget them, etc. A more efficient (and finally simpler) implementation is to reformulate the currents with differential equations and instantaneous updates of variables when a spike arrives. We have a section in the documentation covering this topic: Converting from integrated form to ODEs — Brian 2 2.5.0.3 documentation It also lists the equations that you can use for your model (biexponential synapse) – you have to slightly adapt it, though, to model the current I_syn instead of the post-synaptic potential. Actually, since Brian 1 came with a small library of synapse models that included a biexponential current-based synapse, the equations that you need can be found here as well: Library models (Brian 1 –> 2 conversion) — Brian 2 2.5.0.3 documentation.
I know that this is not as discoverable as it should be, and we had some plans to make some kind of website listing all kinds of models in a nice and accessible way, but well, there are quite a number of things on our todo list…
Hi @mstimberg. Thank you very much. The way that was introduced in the documentation, if I convert the current equations to the differential equations, I’ll have:
eqs = '''
dv/dt = 1/c * (I_syn + I_Na +...) : volt
I_syn = A * g_syn * (v - E_syn) : amp
dg_syn/dt = x : siemens
dx/dt = -1/(tau1*tau2) * ((tau1+tau2)*x + g_syn) : siemens/second
...
'''
and when I define a synaptic object I write something like:
S1 = Synapses(G1, G2, on_pre = 'x += w')
S1.connect()
S1.w = C # some constant
Is this setup correct and will brian accept ‘siemens/second’ as the unit?
Hi @LBENE . Brian would accept siemens/second, but your equations do not look quite right. For a conductance-based synapse (sorry, I overlooked that part of your I_syn definition earlier), your code would follow the formulation here: Library models (Brian 1 –> 2 conversion) — Brian 2 2.5.0.3 documentation
Something like this should work, I think:
invpeak = (tau2 / tau1) ** (tau1 / (tau2 - tau1)) # more efficient to calculate it once
eqs = '''
dv/dt = 1/c * (I_syn + I_Na +...) : volt
I_syn = A * g_syn * (v - E_syn) : amp
dg_syn/dt = (invpeak*x - g_syn)/tau1 : siemens
dx/dt = -x/tau2 : siemens
'''
Thank you very much. I think they are both correct but with different approaches. In my case, I utilized the fact that the equivalent differential equation of g = exp(-t/a) - exp(-t/b)
is g + (a+b) * g' + a*b* g" = 0 (g' = dg/dt, g" = dg'/dt)
in which after defining x as x = dg/dt we get to what I’ve written above.
I checked your approach and that also holds for the differential equation. So, I suppose it’s just two different definitions for x.
Thank you very much. I would really appreciate it if you would elaborate on this. Because my reasoning behind choosing x as the derivative of g_syn is that I don’t want sudden changes like step function for a natural variable. That’s why I’m choosing to have sudden changes for the derivative ('x += w'; because this way x is not a natural variable) which I think is more realistic?
I’m not quite sure what you mean by “natural variable”, here, and why having a step function would be more realistic for a variable representing a derivative of another variable. What I meant when I said that “having a weight as a conductance” feels more natural is that in this formulation, the weight (i.e. step-change in x) has a straightforward interpretation: the maximum of g_syn for a single pulse. If for some reason you wanted to simplify your synapse to a simple exponential synapse dg_syn/dt = -g_syn/tau, then you’d make g_syn jump up by the same value to get the same maximum g_syn. Now, the maximum g_syn conductance is certainly not the best measure of the strength of a synapse under all circumstances, but I’d say it is meaningful. I’m struggling a bit with a comparable interpretation of w in your formulation. It would be the slope (i.e. first derivative) of g_syn at the beginning of the spike, not sure that this is a very intuitive value. But maybe I’m thinking about this the wrong way, happy to stand corrected.
Hi @mstimberg. Thank you very much and sorry for the delayed response. I was trying to implement both approaches and compare the results; however, I figured out that in order to compare the results, we need to calculate equivalent w for each approach.
Anyway, my approach is working alright. So, I suppose that I need to somehow check the results with a simulation of similar specs which has used your approach (which is the standard form).
Regarding what you mentioned, I don’t know the actual name of what I called “natural variable”; however, what I mean is these natural variables don’t jump to a whole different level in one time step. Actually, the main reason that I’m switching to this set-up is because on_pre = 'v_post += w' is a form which implies sudden changes in voltage in just one time step.
Now, I know that action potential has an impulse-like shape at the beginning; but even that jump in voltage is spread between some time steps. In other words, the rise time is not zero.
Lastly, I’m still having trouble figuring out how the standard approach formulation and concept goes. You mentioned the peak for g_syn; yet, I still can’t connect it to the formulation.
Well, strictly speaking, all variables jump to a different level in one time step, since we are doing a time-based simulation. But I get your dislike about the v_post += w formulation. As everything in modeling, there’s a continuum of “realism” in the different approaches. The v_post += w solution is the simplest and often called a “delta synapse”, since the effect of the synapse is modelled as an infinitesimally short event and represented in equations as a Dirac delta. More realistic models instead model a synaptic current or conductance. They in turn might consider a zero rise time, or have a more complex model (alpha function, or biexponential, typically). Whether this is a “realistic” model depends on the time scale of things you are interested in. I’d say that e.g. AMPA synapses are typically considered to have a basically instantaneous rise in conductance/current, so having them jump up immediately seems reasonable. Of course, even with such an instantaneous rise in current/conductance, the membrane potential typically take a few ms to reach its peak. In my formulation above, the synaptic conductance is not jumping up immediately but following a biexponential with a finite rise time. The x variable is jumping up instantaneously, but this is just a “helper variable” – the system of equations is mathematically equivalent to the equation you showed in your initial post.
Regarding my comment about the g_syn peak: If you solve the system of equations g_\mathrm{syn} and x for a single spike at t=0 where you set x(0) = w, you will see that g_\mathrm{syn} peaks at t_\mathrm{max}=\frac{\tau_1\tau_2}{\tau_2-\tau_1}\log\left(\frac{\tau_2}{\tau_1}\right) (ok, I copy&pasted this from Converting from integrated form to ODEs — Brian 2 2.5.0.3 documentation, but it should be correct). If you plug this into the equation for g_\mathrm{syn}(t), you should find that g_\mathrm{syn}(t_\mathrm{max}) = w, i.e. the value you set for the jump in x translates into the peak value of g_\mathrm{syn}.
