So far I've presented three different models for simulating neural spiking dynamics. A key assumption that was made but never stated in each of these examples was that the neuron we were modeling had no defined morphology. In other words, we were looking at models that assumed the neuron was a dimensionless sphere or point. These point neurons can be very effective for studying the behavior of large-scale spiking neural networks (e.g. Izhikevich), but are impractical if you want to investigate how anatomical features of a neuron contribute to signal propagation. For this we return to the Hodgkin-Huxley model and cable theory.
Cable theory and compartmental modeling
Multi-compartmental neuron models attempt to add a level of anatomical basis to the observed physiological dynamics. By approximating dendrites and axons as a series of discrete cylinders, it is possible to build up fairly complex and morphologically accurate neural models. This can then allow for sophisticated simulation experiments that evaluate how signals propagate throughout the neuron based on where a virtual stimulating electrode is placed. One way to evaluate this computationally is to use cable theory, which specifies how electrical currents propagate down a cylindrical section of a neuron.
There are two kinds of signal propagation methods that can be considered -- passive and active. Passive propagation is typically used to model signals arriving in the dendrites and traveling to the soma. As there is nothing reinforcing the signals as they travel, they will slowly degrade depending on the various passive properties of the neural membrane. Active propagation, on the other hand, is typically used for modeling action potentials as they travel down an axon. This can be accomplished by having each axon segment evaluate the Hodgkin-Huxley equations, using the currents received from either side as inputs.
In this example, we're going to demonstrate how action potentials travel down a seven-compartment length of axon using active propagation. Each segment is identical in terms of length, radius, and Hodgkin-Huxley parameters, and a single external stimulus pulse is applied to the center segment.
The action potential traces appear almost like shark fins due to the scaling of the plots. The spike is first elicited in the central compartment, then travels in both directions with consistent velocity due to the uniformity of the compartments. Moving the stimulating electrode will produce a different radiating pattern, while altering the the properties of the segments can affect trasmission speed and efficacy.
The code for this simulation is actually quite similar to the single Hodgkin-Huxley model code. The main trick employed in this example is the usage of linear algebra to greatly simplify the code needed to evaluate the cable equation, so that will be the primary focus.
The connection matrix represents the contribution factor of each segment's current to every other segment. For example, each end segment (represented by the first and last rows of the matrix) adds itself (current going out) and subtracts its only neighbor (current coming in), represented by a 1 and -1, respectively. Each middle segment has two neighbors, so it counts itself twice (current travelling out in both directions) and subtracts each neighbor, represented by -1, 2, -1. A 0 in the matrix indicates that no current flows directly between those two segments. We will see shortly how the addition of this simple representation can have a significant impact.
## connection matrix
Sc = zeros([S,S])
for i in range(S):
if i == 0:
Sc[i,0:2] = [1,-1]
elif i == S-1:
Sc[i,i-1:S] = [-1,1]
Sc[i,i-1:i+2] = [-1,2,-1]
Simulating the cable
The simulation loop in full:
## simulate model
for i in range(1,len(time)):
g_Na = gbar_Na*(m**3)*h
g_K = gbar_K*(n**4)
g_l = gbar_l
m += dt*(alpha_m(Vm[:,i-1])*(1 - m) - beta_m(Vm[:,i-1])*m)
h += dt*(alpha_h(Vm[:,i-1])*(1 - h) - beta_h(Vm[:,i-1])*h)
n += dt*(alpha_n(Vm[:,i-1])*(1 - n) - beta_n(Vm[:,i-1])*n)
dV = -hh(Vm[:,i-1], g_Na, g_K, g_l) - Sc.dot(Vm[:,i-1]) / Ra
dV[elec] += I[i-1]
Vm[:,i] = Vm[:,i-1] + dt * dV / Cm
The first several lines we've seen before in the single compartment case, only this time extended to include multiple compartments. Line 11 here is where we calculate the contributions of all internal currents to all the segments. A single line! There's no mention of the number of compartments either as this formalism can handle an arbitrary number as long as that connection matrix is properly specified. Let's take a closer look at how this works.
First, the ionic channel-based currents are determined by using the hh helper function, which simply evaluates the Hodgkin-Huxley equation for all segments. While conveniently compact, this isn't as interesting as the second half of the expression. In order to factor in the contribution of currents traveling in and out of the compartment, we take the dot product of the connection matrix with the vector of all the current membrane potentials and divide by the intracellular resistance. This is the equivalent of evaluating (Vm[i] - Vm[j])/Ra for every neighboring combination of i and j without the need for a loop or special boundary cases!
Moving on to line 12, we can easily apply the external stimulus to any segment (or multiple, varying segments!) by simply specifying the appropriate index. The last line of code here completes the numerical integration.
Compartmental modeling is an effective tool for exploring morphologically detailed neuron models in simulation. While the example shown here was a simple, uniform, unbranched length of axon, detailed models can have hundreds or thousands of varying compartments representing highly complex dendritic arbors and long branching axons. Attempting to simulate such a model from scratch in Python would be daunting at best. Fortunately, there exist several quality simulator tools, such as NEURON and GENESIS, that take care of the heavy lifting and allow researchers a robust, feature-rich environment for investigating these kinds of models. While these tools are not written in Python, they now have Python-based interfaces that streamline the modeling process even further. I will explore some of these in future posts.