It is possible to create neuron models with a spatially extended morphology, using the SpatialNeuron class. A SpatialNeuron is a single neuron with many compartments. Essentially, it works as a NeuronGroup where elements are compartments instead of neurons. A SpatialNeuron is specified mainly by a set of equations for transmembrane currents (ionic channels) and a morphology.
Morphologies are stored and manipulated with Morphology objects. A new morphology can be loaded from a .swc file (a standard format for neuronal morphologies):
morpho = Morphology('corticalcell.swc')There is a large database of morphologies in the swc format at http://neuromorpho.org/neuroMorpho.
Morphologies can also be created manually by combining different standard geometrical objects:
morpho = Soma(diameter = 30*um)
morpho = Cylinder(length = 100*um, diameter = 1*um, n = 10)The first statement creates a sphere (as a single compartment), the second one a cylinder with 10 compartments (default is 1 compartment). It is also possible to specify the type of process, which is soma, axon or dendrite, and the relative coordinates of the end point. For example:
morpho = Cylinder(diameter = 1*um, n = 10, type = 'axon', x = 50*um, y = 100*um, z = 0*um)In this case, length must not be specified, as it is calculated from the coordinates. These coordinates are mostly helpful for visualization. If they are not specified, 3D direction is chosen at random.
A tree is created by attaching Morphology objects together:
morpho = Soma(diameter = 30*um)
morpho.axon = Cylinder(length = 100*um, diameter = 1*um, n = 10)
morpho.dendrite = Cylinder(length = 50*um, diameter = 2*um, n = 5)These statements create a morphology consisting of a cylindrical axon and a dendrite attached to a spherical soma. Note that the names axon and dendrite are arbitrary and chosen by the user. For example, the same morphology can be created as follows:
morpho = Soma(diameter = 30*um)
morpho.output_process = Cylinder(length = 100*um, diameter = 1*um, n = 10)
morpho.input_process = Cylinder(length = 50*um, diameter = 2*um, n = 5)The syntax is recursive, for example two branches can be added at the end of the dendrite as follows:
morpho.dendrite.branch1 = Cylinder(length = 50*um, diameter = 1*um, n = 3)
morpho.dendrite.branch2 = Cylinder(length = 50*um, diameter = 1*um, n = 3)Equivalently, one can use an indexing syntax:
morpho['dendrite']['branch1'] = Cylinder(length = 50*um, diameter = 1*um, n = 3)
morpho['dendrite']['branch2'] = Cylinder(length = 50*um, diameter = 1*um, n = 3)Finally there is a special shorter syntax for quickly creating trees, using L (for left), R (for right), and digits from 1 to 9. These can be simply concatenated (without using the dot):
morpho.L=Cylinder(length = 10*um, diameter = 1 *um, n = 3)
morpho.L1=Cylinder(length = 5*um, diameter = 1 *um, n = 3)
morpho.L2=Cylinder(length = 5*um, diameter = 1 *um, n = 3)
morpho.L3=Cylinder(length = 5*um, diameter = 1 *um, n = 3)
morpho.R=Cylinder(length = 10*um, diameter = 1 *um, n = 3)
morpho.RL=Cylinder(length = 5*um, diameter = 1 *um, n = 3)
morpho.RR=Cylinder(length = 5*um, diameter = 1 *um, n = 3)These instructions create a dendritic tree with two main branches, 3 subbranches on the first branch and 2 on the second branch. After these instructions are executed, morpho.L contains the entire subtree. To retrieve only the primary branch of this subtree, use the main attribute:
mainbranch = morpho.L.mainThe number of compartments in the entire tree can be obtained with len(morpho). Finally, the morphology can be displayed as a 3D plot:
morpho.plot()Neuronal morphologies can be created by assembling cylinders and spheres, but also more complex processes with variable diameter. This can be done by directly setting the attributes of a Morphology object. A Morphology object stores the diameter, x, y, z, length and area of all the compartments of the main branch (i.e., not children) as arrays. For example, morpho.length is an array containing the length of each of its compartments. When creating a cylinder, the length of each compartment is obtained by dividing the total length provided at creation by the number of compartments. The area is calculated automatically. Complex processes can be created manually by directly specifying the diameter and length of each compartment:
morpho.axon = Morphology(n = 5)
morpho.axon.diameter = ones(5) * 1 * um
morpho.axon.length = [1, 2, 1, 3, 1] * um
morpho.axon.set_coordinates()
morpho.axon.set_area()Note the last two statements: set_coordinates() creates x-y-z coordinates and is required for plotting; set_area() calculates the area of each compartment (considered as a cylinder) and is required for using the morphology in simulations. Alternatively the coordinates can be specified, instead of the lengths of compartments, and then set_length() must be called. Note that these methods only apply to the main branch of the morphology, not the children (subtrees).
A SpatialNeuron is a spatially extended neuron. It is created by specifying the morphology as a Morphology object, the equations for transmembrane currents, and optionally the specific membrane capacitance Cm and intracellular resistivity Ri:
gL=1e-4*siemens/cm**2
EL=-70*mV
eqs='''
Im=gL*(EL-v) : amp/meter**2
I : amp (point current)
'''
neuron = SpatialNeuron(morphology=morpho, model=eqs, Cm=1 * uF / cm ** 2, Ri=100 * ohm * cm)
neuron.v = EL+10*mVSeveral state variables are created automatically: all the variables of the morphology object are linked to state variables of the neuron (diameter, x, y, z, length and area). Additionally, a state variable Cm is created. It is initialized with the value given at construction, but it can be modified on a compartment per compartment basis (which is useful to model myelinated axons). Finally the membrane potential is stored in state variable v. The integration method can be specified as for a NeuronGroup with the method keyword. In general, for models with nonlinear conductances, the exponential Euler method should be used: method = "exponential_euler".
The key state variable, which must be specified at construction, is Im. It is the total transmembrane current, expressed in units of current per area. This is a mandatory line in the definition of the model. The rest of the string description may include other state variables (differential equations or subexpressions) or parameters, exactly as in NeuronGroup. At every timestep, Brian integrates the state variables, calculates the transmembrane current at every point on the neuronal morphology, and updates v using the transmembrane current and the diffusion current, which is calculated based on the morphology and the intracellular resistivity. Note that the transmembrane current is a surfacic current, not the total current in the compartement. This choice means that the model equations are independent of the number of compartments chosen for the simulation. The space constant can obtained for any point of the neuron with the space_constant attribute:
l = neuron.space_constant[0]The calculation is based on the local total conductance (not just the leak conductance). Therefore, it can potentially vary during a simulation (e.g. decrease during an action potential).
To inject a current I at a particular point (e.g. through an electrode or a synapse), this current must be divided by the area of the compartment when inserted in the transmembrane current equation. This is done automatically when the flag point current is specified, as in the example above. This flag can apply only to subexpressions or parameters with amp units. Internally, the expression of the transmembrane current Im is simply augmented with +I/area. A current can then be injected in the first compartment of the neuron (generally the soma) as follows:
neuron.I[0]=1*nAState variables of the SpatialNeuron include all the compartments of that neuron (including subtrees). Therefore, the statement neuron.v=EL+10*mV sets the membrane potential of the entire neuron at -60 mV.
Subtrees can be accessed by attribute (in the same way as in Morphology objects):
neuron.axon.gNa = 10*gLNote that the state variables correspond to the entire subtree, not just the main branch. That is, if the axon had branches, then the above statement would change gNa on the main branch and all the subbranches. To access the main branch only, use the attribute main:
neuron.axon.main.gNa = 10*gLA typical use case is when one wants to change parameter values at the soma only. For example, inserting an electrode current at the soma is done as follows:
neuron.main.I = 1*nAA part of a branch can be accessed as follows:
initial_segment = neuron.axon[10*um:50*um]There are two methods to have synapses on SpatialNeuron. The first one to insert synaptic equations directly in the neuron equations:
eqs='''
Im = gL*(EL-v) : amp/meter**2
Is = gs*(Es-v) : amp (point current)
dgs/dt = -gs/taus : siemens
'''
neuron = SpatialNeuron(morphology=morpho, model=eqs, Cm=1 * uF / cm ** 2, Ri=100 * ohm * cm)Note that, as for electrode stimulation, the synaptic current must be defined as a point current. Then we use a Synapses object to connect a spike source to the neuron:
S = Synapses(stimulation,neuron,pre = 'gs += w')
S.connect(0,50)
S.connect(1,100)This creates two synapses, on compartments 50 and 100. One can specify the compartment number with its spatial position by indexing the morphology:
S.connect(0,morpho[25*um])
S.connect(1,morpho.axon[30*um])In this method for creating synapses, there is a single value for the synaptic conductance in any compartment. This means that it will fail if there are several synapses onto the same compartment and synaptic equations are nonlinear. The second method, which works in such cases, is to have synaptic equations in the Synapses object:
eqs='''
Im = gL*(EL-v) : amp/meter**2
Is = gs*(Es-v) : amp (point current)
gs : siemens
'''
neuron = SpatialNeuron(morphology=morpho, model=eqs, Cm=1 * uF / cm ** 2, Ri=100 * ohm * cm)
S = Synapses(stimulation,neuron,model='''dg/dt = -g/taus : siemens
gs_post = g : siemens (summed)''',pre = 'g += w')Here each synapse (instead of each compartment) has an associated value g, and all values of g for each compartment (i.e., all synapses targeting that compartment) are collected into the compartmental variable gs.
To detect and record spikes, we must specify a threshold condition, essentially in the same way as for a `NeuronGroup`:
neuron = SpatialNeuron(morphology=morpho, model=eqs, threshold = "v > 0*mV", refractory = "v > -10*mV")Here spikes are detected when the membrane potential v reaches 0 mV. Because there is generally no explicit reset in this type of model (although it is possible to specify one), v remains above 0 mV for some time. To avoid detecting spikes during this entire time, we specify a refractory period. In this case no spike is detected as long as v is greater than -10 mV. Another possibility could be:
neuron = SpatialNeuron(morphology=morpho, model=eqs, threshold = "m > 0.5", refractory = "m > 0.4")where m is the state variable for sodium channel activation (assuming this has been defined in the model). Here a spike is detected when half of the sodium channels are open.
With the syntax above, spikes are detected in all compartments of the neuron. To detect them in a single compartment, use the threshold_location keyword:
neuron = SpatialNeuron(morphology=morpho, model=eqs, threshold = "m > 0.5", threshold_location = 30,
refractory = "m > 0.4")In this case, spikes are only detecting in compartment number 30. Reset then applies locally to that compartment (if a reset statement is defined). Again the location of the threshold can be specified with spatial position:
neuron = SpatialNeuron(morphology=morpho, model=eqs, threshold = "m > 0.5",
threshold_location = morpho.axon[30*um],
refractory = "m > 0.4")