I am trying to simulate and study Hodgkin-Huxley dynamics for a spatial / multicompartment neuron given
variable parameter settings while emphasizing its compartments. Now there are two questions or issues which arise for me:
While the HH model seems clear for me when working with non-spatial groups of neurons, I am not quite sure how it works when using a multicompartment model: Does the whole system of equations apply to the whole neuron? Or merely certain parts? If so, which ones or how can I identify these?
How is the current being distributed among the “tree”, how the voltage? I didn’t seem to find any clues in the tutorials, but please feel free to point me in the right direction!
I am trying to see whether certain parts on the dendrites can be changed in values (e.g. the conductances and channel dynamics) and how it might impact the dynamics of the whole structure. But whether I change the parameters or not doesn’t seem to impact my results at all.
I will put a brief example of what I mean below in the “Minimal code” section.
I hope I didn’t oversee something, but if so, please show me where. Or could it simply be that this kind of
feature isn’t implemented yet? I am still new to Brian2 so I might have overlooked some simple things, so I hope you can bear with me here.
Thanks in advance!
Minimal code to reproduce problem
Crucial parts (please tell me whether I should include more parts of the code, I will do this in the next post then!):
# Morphology and specifications:
morpho = Soma(40 * um)
morpho.axon = Cylinder(diameter = 1*um, length = 300*um, n = 100)
morpho.L = Cylinder(diameter = 1*um, length = 150*um, n = 50) # Left dendrite.
morpho.R = Cylinder(diameter= 1*um, length= 150*um, n = 50) # Right dendrite
MorphNeuron = SpatialNeuron(morphology=morpho, model=Ho_Hux,
Cm=1*uF/cm**2, Ri=200*ohm*cm, threshold='v > 50*mV', refractory='v > -10*mV', reset= 'v = -60*mV', method='exponential_euler')
MorphNeuron.I[morpho.R[50*um]] = 0.2 * nA # injecting in the right dendrite
run(5 * ms)
MorphNeuron.I = 0.2 * amp
MorphNeuron.I[morpho.L[100*um]] = 5* nA
#MorphNeuron.I[0] = 0.2 * amp
# Accessing parameters on dendrites and their sections: No influence?
MorphNeuron.L[10*um:30*um].g_Na = 500 * g_L
MorphNeuron.L[30*um:150*um].g_Na = 600 * g_L
MorphNeuron.R[30*um:60*um].g_Na = 50000 * g_L
#MorphNeuron.L[60*um:230*um].g_Na = 700 * g_L # Isn't this of bounds actually?
The equations apply to the whole neuron (but since everything is expressed “per area”, the morphology of the compartments matter), but compartments can have different parameters. For example, in many models you’d only have potassium and sodium conductances in the soma (or in the initial segment of the axon) – in the dendrites, the respective values would be set to zero.
Yes, this should definitely work (see e.g. Example: spike_initiation — Brian 2 2.5.4 documentation which only has non-zero sodium/potassium conductances in the axon initial segment). Could you give a complete code that shows your problem, please? The posted code only calls run before changing the g_Na values, but I guess in your example you had another run statement after that?
So basically there is one defining equation - in my case the HH model - which “collects” the total values? To stay with your example, it calculates the total conductances of the soma plus the zero conductances of the dendritic segment, amounting to the final voltage or current value?
Sure! I hope the whole code will be visible, otherwise I’ll probably have to upload it.
So whether I am adding those values to my dendrites or not doesn’t seem to impact my overall dynamics; from what I gathered the plots remain the same. Hence why I wonder whether I made some simple mistakes somewhere.
I would also like to add synaptic weights at certain dendritic arms later too, by the way. I presume this also is possible?
I did not yet have the time to look into your example in full detail, but here are a few comments that will hopefully be helpful:
For each compartment, the equation you give defines the membrane currents, i.e. the currents flowing through the ion channels. In many models, this would only be a leak current for the passive dendritic compartments. In addition, Brian calculates the currents flowing in from the neighboring compartments. In your example, the soma would receive currents from the initial compartments of the two dendrites (if there is a difference in their respective membrane potentials), each dendritic compartment would receive currents from the compartment before and after, etc. Does that make the general approach clearer?
I think there are a number of mistakes in the equations that probably explain what you are seeing:
in the definition of Im, the currents have the wrong sign, they need negative signs or formulate them as E_Na - v instead of v - E_Na, etc.
your model equation state that each compartment has its own g_K, g_Na, and g_L values, but you never set the g_L values. You have external constants with the same names, but those are ignored as the warning you get mentions. If you want to use the same g_L values for all compartments, delete g_L : siemens * meter**-2 from the equations – Brian will then use the g_L value defined outside of the equations. Conversely, if you want to set g_Na and g_K for each compartment separately, delete the respective external constants.
some of the quantities are off by orders of magnitude, I think. E.g. g_L is defined as 0.03*msiemens*meter**-2, whereas it’d be typically ~10000 higher (e.g. 0.3*msiemens*cm**-2). Similarly, the line MorphNeuron.I = 0.2 * amp injects a very high current into all compartments – this probably breaks the numerical integration completely.
So in spatial models, the system of equations governs pretty much every compartment and also takes into account its neighbouring compartments? If that is so yes, that is clear now.
May I ask why this is reversed? If I check the literature (e.g. Dayan and Abbott 2001), the current equation shows me “V-E_Na” etc. Maybe I overlooked something?
Otherwise I took into consideration that feedback, though I still don’t see changes in the dynamics when altering values on dendritic arms…
Oh, good point, I did not realize that we don’t mention this specifically in the documentation. The convention used by Dayan & Abbott (and many others) is that positive currents depolarize the neurons, except for membrane currents where the sign convention is the other way round. I think this comes from the convention that Hodgkin&Huxley used and is probably due to the fact that you measure those currents in voltage clamp experiments, i.e. you actually measure the inverse current (the current needed to keep the voltage at the same value). In Brian, we decided to use the “more logical” convention that all currents follow the same sign convention, i.e. positive currents depolarize the cell. Some other simulators (e.g. GENESIS, if I remember correctly) do the same, but of course we should mention this clearly in the documentation!
Ah, that explains it! Good to know, and yes I agree that it might avoid confusion if this were specified in the documentation.
Sure, I will post the full current version, with plots and all (maybe I made some errors there, too?)
Also, another issue I can’t seem to avoid is that my voltage doesn’t cross -60 mV in the dynamics. I wonder if I forgot to provide some input and whether this is related to my issue of no dendritic influence on the dynamics?
Hi @DenisS I had another look at the code, and there are still (at least?) three major issues that affect the results:
the value for the leak conductance that is also used to set the Na conductance is still orders of magnitude too low (/m^2 instead of /cm^2, i.e. off by a factor of 10000)
the leak conductance of the individual compartments is not set, and therefore 0 everywhere
the current injected via MorphNeuron.I = 0.002 * amp into all compartments is much too strong, and orders of magnitudes stronger than the currents injected into parts of the dendrites.
This drowns out everything related to the active conductances you are setting. In general, I’d recommend changing one thing at a time to verify that everything is working as you expect. In your example above, after 15 ms, you add active Na channels of various densities to several places in the dendrite, inject a (way too strong) current into all compartments, and a weaker current into a specific compartment. It is unclear what to expect in this situation.
Below, a simpler code that only injects a current into parts of the right dendrite between 10-15ms, with no other current stimulation, and running the same simulation twice once without any active channels, and once with Na + Kd channels in parts of the right dendrite. As you can see, having the active channels makes a clear difference – not claiming that the specific values model anything of biological interest, of course. Please also don’t forget that the classical HH equations are for the axon of a squid, not for a neuron.
Indeed, I was planning to extend the HH model anyways and / or start to implement the Connor-Stevens model at a later time.
I think I can work with that, and will thus mark your last reply as solution. However, I still might have questions over time. Can I ask them here if they still fit under the “umbrella topic”? Or is it better to start a completely new topic?
And thank you very much for your clear answers as well as patience! Highly appreciated.
as in my code from before, does it mean that the current is injected over all the compartments on the right dendrite which has the length of 50 * um? And then at a later point, when using the monitor and trying to plot the results, I am extracting all the recorded voltages per time unit from each compartment?
Because what I actually want is to have the voltage for every compartment at every time t, not for instance an “average” voltage which arrives at the soma in the end. I would gladly explain this in greater detail if my question isn’t clear enough, though!
I think it relates to question 1, but in the definition of my morphology, I can define lengths and diameters. So for instance if I define a dendrite with say n = 50 compartments and also its length = 100 um, how does the length tie in with the compartments? Is it averaged too, as in this case one compartment having a length of 2 um?
Because what I want in the end is - as first simple example - one dendrite with fixed length (and thus fixed compartments) and for each compartment the voltage, current and channel / gate probabilities; and all of that depicted through different graphs. Currently, my interpretation is that I only see the “end result” processed at the soma rather than also seeing intermediate results which happen in every single compartment.
I also might have some follow-up questions, but maybe they will be answered afterwards already. We will see.
If you use a e.g. morpho.R[50*um], this means the compartment that includes the “50 µm from start” point in it. So if you inject a current there, it would be only injected in that specific compartment, and if you record from it, it will give you only the value for that compartment. What happens internally is that morpho.R[50*um] will return the index of the respective compartment, and something like morpho.R[50*um: 100*um] would return a list of compartment indices. If your morphology is very simple, let’s say a spherical soma and a single dendrite of 10 compartments, then these indices are straightforward (the soma would have index 0, the dendritic compartments would have indices from 1 to 10), but things get more complex when you have branching, and in either case it is nice that you can use this syntax without caring about the actual number of compartments.
You certainly have access to voltages at all compartments individually, you might be interested in the examples here: Examples — Brian 2 2.5.4 documentation (most of them actually directly specify indices, instead of indexing with the morphology and distances…).
If you use the Cylinder class, then it indeed divides up everything into compartments of the same size. If you want detailed control about each compartment, you can use the Section class instead, where you specify length&diameter for each compartment. See Multicompartment models — Brian 2 2.5.4 documentation for a visualization.
So if I had, for instance, code as “mon_dendrite[morpho.dendrite[50 * um]].m” it will exactly return the value(s) of gating probability m of the compartment(s) "from 0 up to 50um. So if I’d like to have the values of my whole dendrite (of say length = 100 um), I’d simply need to specify and put in 100 um etc. On the other hand, having only certain compartments and their respective lengths, I’d simply need to specify them via indexing (for example “50um : 100 um”).
In my specific case, however, it might be more useful to specify the specific compartments rather than the length itself.
Hope I summarized your explanation correctly.
Yes, I will play around with it more, but if the above explanation is correct it should be clear.
Perfect, yes this might be something I may use in future for more complex morphologies.
In that case, a follow up question: So far, I only created one to two dendrites which then may receive some input and deliver this to the soma (for now I don’t need any axon transporting the signal further) while studying its dynamics.
But would it be possible in Brian for two dendrites (or even more) to “branch into each other” such that, for example, their EPSPs would “meet” at an intersection, be processed, and their result propagated towards the soma? My ad hoc idea would be to model this point of intersection as another soma which then processes this and delivers it to the actual soma (via an axon), but I am not sure if the behaviour / dynamics would be the same or whether there exist more elegant solutions to this.
Not exactly, morpho.dendrite[50*um] only refers to a single compartment, the compartment that contains the 50µm point. To get the values for all compartments from 0 to 50µm, you’d need to use morpho.dendrite[:50*um].
If you recorded from the whole dendrite and you want to get all values, you’d not need any indexing at all.
Maybe things are clearer with code. In the example below, I create a soma with a single dendrite and I inject a current for 1ms at the end of the dendrite. On the left, I plot the values of all recorded compartments (no indexing needed) at different points in time. On the right, I plot the membrane potential over time for three compartments (each indexed with their distance from the start). If you change the numbers of compartments on top, nothing should change in the plot, since all indices are based on distances and therefore independent of the number of compartments.
from brian2 import *
n_compartments = 11
morpho = Soma(30*um)
morpho.dendrite = Cylinder(diameter=5*um, length=100*um, n=n_compartments)
gL = 1e-4*siemens/cm**2
EL = -70*mV
eqs = '''Im = gL * (EL - v) : amp/meter**2 # simple, passive model
I : amp (point current)'''
neuron = SpatialNeuron(morphology=morpho, model=eqs, method="exponential_euler")
neuron.v = EL
# Record from dendrite
mon = StateMonitor(neuron, 'v', record=morpho.dendrite)
neuron.dendrite.I[-1] = 1*nA # Inject current at the end of the dendrite for 1 ms
run(1*ms)
neuron.I = 0*nA # Switch off stimulation
run(4*ms)
fig, axs = plt.subplots(1, 2, figsize=(10, 4), sharey=True, layout='constrained')
# plot the membrane potential across the dendrite for different times
cmap = matplotlib.cm.get_cmap('winter')
for t in np.arange(1, 5)*0.5*ms:
timestep = int(round(t / defaultclock.dt))
color = cmap(t/(2*ms))
axs[0].plot(neuron.dendrite.distance / um, mon.v[:, timestep] / mV, label=f't = {t}', color=color)
axs[0].set_xlabel('Distance (µm)')
axs[0].set_ylabel('v (mV)')
axs[0].legend()
# plot the membrane potential at specific positions over time
positions = [30, 60, 90]*um
for position in positions:
axs[1].plot(mon.t/ms, mon[morpho.dendrite[position]].v / mV, label=f'{position/um} µm')
axs[1].set_xlabel('Time (ms)')
axs[1].legend()
plt.show()
So I also assume that if I wouldn’t use lengths, but simple indexing of compartments, it would work analogously: morpho.dendrite[ x ] would get all the value(s) from compartment no. x, whereas [y:x] for instance would contain all the values from y to x. And if I’d be interested in some changes or differences between two compartments, I could simply extract values from these compartments and calculate the differences. Would this be a correct explanation / interpretation?
I see, yes, this is clear!
And one last question for this thread, if you allow: In my HH model above - or rather the Spatial neuron where I define it - (see Simulating (semi-)morphological dynamics in the Hodgkin-Huxley model - #7 by DenisS) the refractory, reset and threshold didn’t seem to have much of an impact when I removed their conditions from the spatial neuron. I believe they must be adjusted to a proper voltage given my HH model which I Initially achieved by increasing R_i, but I am still not sure about this way to “solve” it. Any ideas on what I might have overlooked?
Thank you once again for the explanations and examples.
A neuron model like the HH model that models the membrane potential during an action potential shouldn’t have a reset statement. This is meant for models like the integrate-and-fire model, where a spike artificially resets the membrane potential to model the spike. In a HH model (if you remove the reset, as you should), the threshold and refractory argument can be used, but they do not affect the simulation directly. There only use case would be to register the spiking activity with a SpikeMonitor, or have other things triggered by discrete spiking events (e.g. synapses).
Yes, I have taken another deeper look at the model itself and understand - to put it in simple terms - how there is no “real” threshold but rather an emulation of one via the equations themselves (given certain input current profiles) which can probably be defined more precisely via those two arguments.
In any case, I have worked up my code and also fitted the parameters to move away from the Giant squid axon. Everything works - technically!
I still have one decisive issue regarding the interpretation of my results (I will post my code below): The seemingly instantaneous nature of the dynamics. For example, I would expect the “trigger” of an action potential at the soma at a later point and not relatively close around the start where the current is injected in the first compartment of my dendrite - which I interpreted to be at the top, with the signal traveling through the other compartments and finally arriving at the soma.
I would also expect some sort of delays at the single compartments themselves rather than seeing a spike at each compartment I have tested. So for instance, I would expect to see some higher activity at compartment x with different conductances and then another one possibly barely active dynamics at compartment y. So the results don’t tell me much, really and in this particular case there is simply too much alignment with “homogeneous” channels.
I will also continue to search for the errors myself and would report them, as I suspect it might be due to some simple plotting errors, erroneous parameter settings or something simple I just happened to overlook.
Hi @DenisS . There seems to be an issue with your equations: if you look at the dynamics without any input, it will fire an action potential and then stay at a high membrane potential. This even happens in the absence of input, so it is unrelated to any propagation (and therefore “instantaneous” everywhere). Maybe check the rate functions with the technique for plotting functions described in the documentation . There’s at least one mistake (but fixing it doesn’t seem to be enough):
The initial action potential unrelated to input masks everything else a bit, you do actually see a small delay in the propagation of the current that you inject at 10ms.
The first compartment is the one closest to the soma, you would want to inject into the last compartment instead.
As I mentioned above, the action potential you see is actually unrelated to the current injection and happens everywhere at the same time.
Note that in the inhomogneous case, you only change the sodium conductance and therefore affect the dynamics mostly when the membrane potential is above the threshold of action potential generation. Close to the resting potential, the Na channels are closed, so increasing their maximal conductance won’t have much effect.
Also note that the legend of your plots contains a mistake: the m and h gating variables are for the Na channels, and the n is for K.
Hey Marcel, thanks! I have fixed the smaller issues and the only culprit for everything working properly still is the occurrence of the AP without any sort of input current / pulse paradigm.
After fixing some of the smaller issues, I tried to see whether having a starting voltage of -65mV would have some meaningful impact, though I only find a similar picture: In absence of any pulse, there is still some sort of spike occuring (albeit not one with a high amplitude). Checking the Sodium dynamics also revealed that the m-Channels have an unusually high probability directly at the beginning, while the others are close to 0. This can’t be the case in any HH-oriented scenario and I believe there must lie the issue; so there may still be an error within the equations, but I have yet find where …
EDIT: After testing / using the rate functions and parameters from here, the dynamics is definitely different and aligns more with what I would expect (without any input current), corroborating the suspicion that my current rating functions have some issues.
