# Single Neuron Computational Modeling

## Summary and Keywords

Computational modeling is essential to understand how the complex dendritic structure and membrane properties of a neuron process input signals to generate output signals. Compartmental models describe how inputs, such as synaptic currents, affect a neuron’s membrane potential and produce outputs, such as action potentials, by converting membrane properties into the components of an electrical circuit. The simplest such model consists of a single compartment with a leakage conductance which represents a neuron having spatially uniform membrane potential and a constant conductance summarizing the combined effect of every ion flowing across the neuron’s membrane. The Hodgkin-Huxley model introduces two additional active channels; the sodium channel and the delayed rectifier potassium channel whose associated conductances change depending on the membrane potential and that are described by an additional set of three nonlinear differential equations. Since its conception in 1952, many kinds of active channels have been discovered with a variety of characteristics that can successfully be modeled within the same framework. As the membrane potential varies spatially in a neuron, the next refinement consists in describing a neuron as an electric cable to account for membrane potential attenuation and signal propagation along dendritic or axonal processes. A discrete version of the cable equation results in compartments with possibly different properties, such as different types of ion channels or spatially varying maximum conductances to model changes in channel densities. Branching neural processes such as dendrites can be modeled with the cable equation by considering the junctions of cables with different radii and electrical properties. Single neuron computational models are used to investigate a variety of topics and reveal insights that cannot be evidenced directly by experimental observation. Studies on action potential initiation and on synaptic integration provide prototypical examples illustrating why computational models are essential. Modeling action potential initiation constrains the localization and density of channels required to reproduce experimental observations, while modeling synaptic integration sheds light on the interaction between the morphological and physiological characteristics of dendrites. Finally, reduced compartmental models demonstrate how a simplified morphological structure supplemented by a small number of ion channel-related variables can provide clear explanations about complex intracellular membrane potential dynamics.

Keywords: cable equation, Hodgkin-Huxley equation, neuro-electrical circuit, compartmental modeling, conductance-based models, dendritic tree

Single Neuron Computational Models in Neuroscience

Since Cajal’s description of their detailed morphology, using a modified version of Golgi’s silver staining technique, neurons are known to communicate by receiving signals through their dendrites from the adjacent axons of other neurons. These signals change the dendritic membrane potential by means of ions flowing to and from the extracellular space. The resulting dynamical update of the membrane potential is central to understanding how neurons work and communicate with each other. In the long years that have passed since Cajal’s discovery, experimental techniques have rapidly advanced. As a result, we understand that a single neuron possesses many complex specializations involved in the integration of neural signals and leading to a rich variety of outputs. Therefore, computational modeling must supplement experiments to study how the complex structure of neurons relates to the processing of stimuli, and ultimately animal behavior.

Indeed, in the past 30 years, computational modeling has become an essential ingredient to connect the different structural components of a neuron and comprehend the mechanisms by which they allow its function. For example, neurons possess different ion channels that can be identified and characterized using experiments but the mechanisms by which these ion channels interact to generate specific membrane potential events such as action potentials is studied using modeling. Modeling also helps to investigate the functional role played by the different compartments of a neuron. Advanced imaging techniques can nowadays reveal the activation of neuronal compartments with remarkable detail, but how neuronal morphology contributes to signal propagation or attenuation is often studied using modeling.

In this article, we will start by describing how membrane properties are translated into the components of an electrical circuit and how this circuit is in turn converted into mathematical equations. By proceeding from simpler to more complex models, we will learn how different biophysical components of a neuron are integrated in a unified and comprehensive model. Finally, a selection of topical studies is discussed to demonstrate how single neuron computational modeling has been used to shed light on neuronal function.

Main Neuron Compartments and Fundamentals of Neural Signaling

As illustrated in Figure 1A, a neuron can generally be subdivided into three compartments: dendrites, axon, and soma (cell body). The incoming input is usually received by the dendrites through specialized contacts or connections called synapses. Each synaptic connection involves an axon that carries the neuronal output from other, pre-synaptic, cells to the dendrites of the neuron under consideration (post-synaptic). Within a neuron, the process of converting the synaptic input into an output signal suitable for sharing with other neurons is called synaptic integration. Figure 1B is a schematic of the lobula giant movement detector neuron (LGMD) of locusts. This neuron is located within the optic lobe, the brain structure specialized in processing visual stimuli in insects, where it receives thousands of synaptic inputs originating from the retina of the animal. The synaptic input can be subdivided into excitatory and inhibitory, which in the case of the LGMD impinge on different extended dendrites. The LGMD releases its output signals through a synapse to the descending contralateral movement detector neuron (DCMD) which in turn sends signals to motor circuits that control the legs and wings. Depending on the kind of visual stimuli it receives, the output signals of the LGMD neuron determine whether a locust will stay in place or rapidly generate an escape behavior to avoid a potential threat such as a predatory bird.

Action potentials are the major mode of communication between neurons resulting from synaptic integration. An action potential is a rapid and transient excursion of the membrane potential from its resting value to a peak, approximately 90 mV above rest, lasting about 1 ms. For an action potential to be generated, the cell’s membrane potential needs to be depolarized above a threshold level, typically 10–20 mV above rest, resulting in a positive and then regenerative feedback signal generating the action potential. As synaptic inputs can be either excitatory or inhibitory, depolarization above threshold requires the amount of excitation to exceed that of inhibition. Furthermore, the excitatory and inhibitory postsynaptic potentials (EPSPs and IPSPs) need to reach the location where the action potentials are generated, which is called the spike initiation zone (SIZ), usually located in the initial segment of the axon. The propagation of EPSPs and IPSPs is controlled by many factors including the neuronal membrane properties, dendritic morphology, and different types of ion channels or conductances. The following sections will explain how each of these factors is studied using modeling methods.

An Electrical Circuit Based on Neuronal Membrane Properties

## Basic Membrane Biophysics

The lipid bilayer composing the neuron’s membrane is a good electric insulator preventing ions from flowing in and out of the cell. This impermeability to charged molecules makes the cell membrane act as a capacitor by separating charge and allowing it to accumulate on the membrane’s interior and exterior surfaces. Ion channels and ion pumps are transmembrane proteins with selective permeability toward specific kinds of ions, for example sodium, potassium or calcium. Ion pumps, on the one hand, use energy to transport ions against their concentration gradients. Their role is to maintain the equilibrium concentration level of these ions constant within the neuron. Ion channels, on the other hand, have gates that open under specific conditions and selectively pass ions through the membrane usually following their concentration gradients (controlled by ion pumps). Thus, ion channels have a low resistance for the specific ions they let pass and the resulting conductance depends on the density and the type of ion channels inserted in the membrane. Figure 2 pictures these biophysical components embedded within the neuronal membrane.

The net effect of the imbalance in ion concentration maintained by pumps is to create, for each ion, an electrical potential difference between the inside and outside of the cell at which ions flowing in and out of the cell will cancel out. This is the membrane potential at which the electrical gradient across the membrane exactly cancels the diffusion gradient for the ion under consideration and is called the Nernst potential. It is determined from the internal and external ion concentrations by the Nernst equation. For example, in the case of sodium ions:

where $k$ is Boltzmann’s constant, $T$ the temperature, $z$ the valence of the ion under consideration (e.g., +1 for Na^{+}), e the elementary charge and ${\left[N{a}^{+}\right]}_{o/i}$ the sodium concentration outside/inside the cell, respectively. If we use concentration values typical of the giant axon neuron in squid, ${\left[N{a}^{+}\right]}_{i}=50$ mM, ${\left[N{a}^{+}\right]}_{o}=440$ mM, we obtain ${V}_{Na}=+56$ mV. In contrast, ${V}_{K}=-77$ mV. Since, typically, at steady state the electrical potential of the neuron measured inside relative to outside amounts to -70 mV, $N{a}^{+}$ ions will flow inside the cell when a $N{a}^{+}$ channel opens whereas ${K}^{+}$ ion will flow outside when a ${K}^{+}$ channel opens. In each neuron, the total membrane current is the sum of currents generated by all different types of ion channels. It is often measured per unit area to compare neurons of different sizes.

The membrane potential measured in different locations within a neuron may vary. For example, when different types of synaptic inputs are activated the resulting gradient results in ion flow along the intracellular medium to cancel these differences. While ions flow from one place to another, the resistance created by the cytoplasm against this longitudinal current, ${I}_{L}$, depends on the length of the segment between any two locations and its cross-sectional area. The associated specific resistance called the axial resistivity, ${R}_{A}$, is measured in $\text{\Omega}cm$. Using Ohm’s law, the difference of membrane potentials is the products of ${I}_{L}$ and ${R}_{A}$, scaled by segment length and cross-sectional area. In addition to the axial resistance, spatial variations originating in the morphology of dendrites or how ion channels are distributed can influence membrane potential and, consequently, how synaptic inputs are processed as they propagate to different neuronal compartments and an appropriate output is generated.

## Equivalent Electrical Circuit Components in Mathematical Terms

How much current is needed to change the membrane potential at a given rate can be calculated using the membrane capacitance, ${C}_{m}$, usually measured in $\mu F/c{m}^{2}$. If $V(t)$ represents the membrane potential (inside relative to outside, in $mV$) and $Q$ the excess charge inside, then ${C}_{m}=\frac{Q}{V}$, $I=\frac{dQ}{dt}$ and thus $I={C}_{m}\frac{dV}{dt}$ (in $\frac{\mu A}{c{m}^{2}}$), where the current $I$ is positive when flowing toward the outside.

The current generated by ion channels has the form ${g}_{ion}(V\left(t\right)-{V}_{ion})$. Here ${g}_{ion}$ is the conductance per unit area due to the channels and ${V}_{ion}$ is the Nernst or reversal potential associated with the ion flowing through the channel. The current ${g}_{ion}(V\left(t\right)-{V}_{ion})$ moves $V(t)$ toward ${V}_{ion}$ ; if $V\left(t\right)<{V}_{ion}$, then it depolarizes a neuron (negative or inward current); if $V\left(t\right)>{V}_{ion}$, then it hyperpolarizes it (positive or outward current); and if ${V}_{ion}$ is near $V(t)$ then it passes little net current.

Single Compartmental (Isopotential) Model

The model in which a whole neuron is described as having a uniform membrane potential independent of location is called the isopotential model. This model is not concerned with the longitudinal current flowing from one location to another within the neuron. When no synaptic input is received, the total current flowing across the membrane originates from ion channels; so, by using Kirchhoff’s current law:

${C}_{m}\frac{dV}{dt}$ is the capacitive current density and $\sum}_{i}{g}_{i}(V\left(t\right)-{V}_{i})$ is the sum of currents generated by identified ion channels with conductance ${g}_{i}$ and reversal potential ${V}_{i}$ for the ${i}^{th}$ channel. The term ${g}_{Leak}\left(V\left(t\right)-{V}_{Leak}\right)$ is called the leakage current. It represents the bulk of the resting conductance, *g _{Leak}*, and its associated reversal potential,

*V*, that are relatively constant and may be attributed to unspecified ion pumps and unidentified channels. The steady-state membrane potential $V$ that makes ${I}_{ss}={g}_{Leak}\left(V-{V}_{Leak}\right)+{\displaystyle \sum}_{i}{g}_{i}(V-{V}_{i})$ zero is called the resting potential and ${I}_{ss}$ is called the steady-state current. As explained above, the resting potential is approximately equal to -70 mV in most neurons.

_{Leak}To investigate the characteristics of ion channels and their associated membrane currents, experimentalists often inject into a cell through an electrode a current of fixed amplitude, ${I}_{stim}$. Since every term in Equation 1 is scaled per unit area ($\frac{\mu A}{c{m}^{2}}$), and ${I}_{stim}$ represents the current injected into the whole neuron, one needs to multiply the left-hand side of Equation 1 by the total surface area of the neuron and reapply Kirchhoff’s current law, yielding

Alternatively, if synaptic input is considered, its current, ${I}_{syn}$, can be implemented with the conductance, ${g}_{syn}$, and the reversal potential, ${V}_{syn}$, yielding

Figure 3 depicts the electrical circuit associated with a one compartment model comprising a leakage conductance and one ion conductance. Depending on the source of input, current injection or synaptic input, the circuit needs to be modified to reflect either Equation 2 or Equation 3.

## Different Types of Ion Channels

Three broad types of ion channels can be distinguished based on the factors that affect the channel conductance, ${g}_{i}:$ passive, active, and ligand-gated. Passive channels/conductances are independent of the membrane potential and thus follow Ohm’s law. One example already mentioned is the leak conductance. Active channels are also called voltage-gated channels because their conductance is a function of the membrane potential. Active conductances are modeled by a combination of functions that describe the probability of activation and/or inactivation of the gates controlling the opening of the channels. Furthermore, the open probability and time course of opening depend on the gates’ opening and closing rates that are a function of $V$. In the next sections, explicit descriptions of active conductances associated with different kinds of channels are provided. As their name implies, ligand-gated channels require the binding of a ligand for opening. Neurotransmitter released following the action potential of a presynaptic neuron is one example of a ligand.

Active Isopotential Model: Hodgkin-Huxley (HH) Equations

A single channel fluctuates between the open and closed states in a probabilistic manner and the probability of opening depends on $V$ at a given time. The probability that one channel is in the open state, $P(V)$, approximates the fraction of the population of the channel type that is open by the law of large numbers, if we assume that their random fluctuations are independent of each other. As mentioned earlier, the maximum conductance per unit area, $\overline{g}$, is the conductance when every channel of the given type is open, $P\left(V\right)=1$, so the bulk conductance is given by $g\left(V\right)=\overline{g}P(V)$.

## How Active Channels Work

Hodgkin and Huxley described how the permeability of sodium and potassium channels changes as the membrane potential changes over time using the voltage clamp technique on the giant axon of the squid and expressed it as a mathematical model called the Hodgkin-Huxley (HH) model (Hodgkin & Huxley, 1952),

where

In this equation, ${g}_{K}$ and ${g}_{Na}$ are two different types of conductances, a persistent potassium conductance and a transient sodium conductance. Conceptually speaking, persistent conductances have one gate that opens and closes according to the value of $V(t)$. The opening of the gate is called activation and its closing is called de-activation. The open probability of ${g}_{K}$, ${\left[n\left(V\right)\right]}^{p}$, implies that to open the gate, $p$ independent events with probability, $n(V)$, are required. In the Hodgkin-Huxley model $p$ was chosen as 4 to best fit experimental data and should not be interpreted as reflecting the real substructure of the gate.

In the model for the sodium channel, the term ${\left[m\left(V\right)\right]}^{q}$ plays a similar role representing activation. However, transient conductances have an extra gate that can be closed during activation. This gate is represented by the term ${\left[h(V)\right]}^{r}$. The closed state of this extra gate is called inactivation and its open state is called de-inactivation. Thus, ${\left[m\left(V\right)\right]}^{q}$ is the probability of activation and ${\left[h(V)\right]}^{r}$ that of de-inactivation and the product of the two, the open probability, is the fraction of the channel population which passes $N{a}^{+}$. The exponents $q=3$ and $r=1$ are also subject to data fitting. Figure 4A and 4B illustrates in a cartoon example how such gates work.

Each activation and inactivation variable can be constructed from voltage-dependent rates of transitioning between open and closed states. For example, let us consider the activation variable $n(V)$ since $m\left(V\right)$ and $h(V)$ are constructed using the same conceptual reasoning. Let ${\alpha}_{n}(V)$ be the rate of gate opening from its closed state and ${\beta}_{n}(V)$ be the rate of gate closing from its open state. The gate’s rate of change over time, $\frac{dn}{dt}$, is the difference of the rate of opening and closing so that:

By multiplying Equation 5 by $\frac{1}{{\alpha}_{n}\left(V\right)+{\beta}_{n}(V)}$, we get another useful form:

Equation 6 implies that for fixed $V$, $n$ approaches its steady state value ${n}_{\infty}$ exponentially with the time constant ${\tau}_{n}\left(V\right)$. Figure 5A illustrates these functions ${n}_{\infty}$ and ${\tau}_{n}\left(V\right)$ as a function of membrane potential. The steady-state activation and inactivation of the sodium channel are illustrated in Figure 5B and 5C, respectively, together with their time constants (see Hodgkin & Huxley, 1952, for the corresponding formulas).

For transient conductances such as ${g}_{Na}$, the relative magnitude of the time constants of the activation variable, $m(V)$, and the inactivation variable, $h(V)$, play an important role especially for fast changing $V\left(t\right)$. The scale of ${\tau}_{m}$ is about 1/20 of the scale of ${\tau}_{h}$ which means that the transition of activation from closed to open is much faster than the transition of inactivation from open to closed for a sudden change of $V(t)$. For example, if the membrane potential suddenly jumps from -75 mV to -25 mV, this will result in a large and transient increase in ${g}_{Na}$ since activation proceeds much faster than inactivation. As a result, sodium ions will flow into the cell, causing a transient depolarization that activates ${g}_{K}$, which in turn results in ${K}^{+}$ ions flowing out of the cell, thus returning the membrane potential toward its resting value. Such a rapid sequence of events leading to transient opening of $N{a}^{+}$ channels followed by sustained opening of ${K}^{+}$ channels is precisely the cause of the action potential.

## Numerical Solution of the HH Model

Equation 4 together with the equations governing the gates $n,m$ and $h$ (Equations 5 or 6) are collectively named the HH equations. As the HH equations are nonlinear differential equations (at least the first one of them), they need to be solved numerically. Hines (1984) modified the backward Euler method to achieve second-order accuracy and yet remove the great cost of solving a nonlinear system for every iteration. His method is called a staggered Euler scheme because the discretization and update of the variables $n,m$ and $h$ is staggered by half a time step with respect to the discretization and update of $V$. For $j=1,\dots N$, this amounts to the following scheme, where for conciseness we omit $m$ and $h$ which are treated similarly to $n$:

with initial conditions:${V}^{1}={V}_{r}$ and ${n}^{1}={n}_{\infty}({V}_{r})$, ${m}^{1}={m}_{\infty}({V}_{r})$, ${h}^{1}={h}_{\infty}({V}_{r})$.

Every time-dependent variable except ${V}^{j}$ is updated following the trapezoid rule,

and solved for at the ${j}^{th}$ iteration:

Then ${V}^{j-1/2}$ is calculated using backward Euler’s method from ${V}^{j-1}$ using updated ${n}^{j},{m}^{j},{h}^{j}$, and ${I}^{j}$ values with a step size $\frac{dt}{2}$. Finally, the next half step is obtained using forward Euler, ${V}^{j}=2{V}^{j-1/2}-{V}^{j}$ (see Gabbiani & Cox, 2017 for details).

## The HH Model Generalizes to a Variety of Active Channels

Following the successful prediction of the action potential by Hodgkin and Huxley, similar conductance-based models were shown to capture the characteristics of subsequently discovered channels. Indeed, a combination of genetics and electrophysiological methods has demonstrated that neurons can potentially express a vast repertoire of ion channels. Even for one type of ion, many channels with different properties are known. For instance, the Connor-Stevens model (Connor & Stevens, 1971) added a transient ${K}^{+}$ current called the A-type current (${I}_{A}$) to the HH model. The main characteristic of this current is that it activates very fast at subthreshold potentials and can thus delay $N{a}^{+}$ channel opening and the action potential, provided ${I}_{A}$ is not inactivated. As a result, the relationship between the amplitudes of current injection steps and the neuron’s firing rate transitions from a discontinuous jump (e.g., 0 to 15 spk/s) at a given threshold current (Type II response) in the HH model to a linear increase (Type I response) in the Connor-Stevens model. In another example, a transient $C{a}^{2+}$ conductance paired with a fast $N{a}^{+}$ conductance could reproduce the characteristic clustered spike sequence, also called a burst of spikes, generated by thalamic relay cells following a hyperpolarizing current step (Wang, 1994). There is also a variety of channels that are voltage-gated but depend as well as on other factors such as the concentration of $C{a}^{2+}$. For example, the voltage and $C{a}^{2+}$ activated ${K}^{+}$ channels (BK channels) are sensitive to low $C{a}^{2+}$ concentration as well as the membrane potential (Stuart, Spruston, & Häusser, 2016).

Synaptic Conductances

When an action potential reaches the presynaptic terminal, neurotransmitter is released into the synaptic cleft between the pre- and post-synaptic sides of the synapse. With neurotransmitter binding to the postsynaptic neuron, the end-point is reached for the process of passing a signal from one neuron to the next and signal propagation in the postsynaptic neuron starts. The ligand-gated receptors in the postsynaptic neuron are activated by the neurotransmitter either through the transient opening of a synaptic conductance (ionotropic receptor) or indirectly by triggering a second messenger cascade (metabotropic receptor) before the ligand unbinds and is cleared from the synaptic cleft. Synaptic input can either be excitatory, depolarizing, or inhibitory, hyperpolarizing. Glutamate and GABA ($\gamma $-aminobutyric acid) are the principal excitatory and inhibitory neurotransmitters in the mammalian central nervous system, respectively. The main ionotropic receptors for glutamate are AMPA (α-amino-3-hydroxyl-5-methyl-4-isoxazole-propionate) and NMDA (N-methyl D-aspartate) whose conductance is voltage-dependent. Both receptors are $N{a}^{+}$ and ${K}^{+}$ selective with a reversal potential near 0 mV, but NMDA receptors also selectively pass $C{a}^{2+}$. Inhibitory receptors activated by GABA include the A-type ionotropic receptor (GABA_{A}) that opens a conductance selective to $C{l}^{-}$ and has a reversal potential around -68 mV.

Synaptic conductances are generally modeled in two ways. First, with a simple time dependence: ${g}_{syn}={\overline{g}}_{syn}f(t)$ where $f\left(t\right)$ is a function of time implementing a transient rise and decay of the conductance over a short time interval (5-20 ms) such as a step function or an $\alpha $-function. Second, synaptic conductances can be modeled using the law of mass action by tracking the concentrations of open and closed receptors, as well as the neurotransmitter binding, unbinding and clearance from the synaptic cleft.

Cable Equation and Passive Compartmental Model (Dendrite)

As experimental techniques using microelectrodes made it possible to obtain intracellular recordings from the soma of neurons in the 1950s, Rall pointed out that an isopotential neuronal model could not explain the experimental data and that a neuron needs to be treated as a cylindrically-shaped electrical cable. The idea is that electrical stimulation spreading in a neuron behaves like the current passing through a cable that may be broken down into nearly isopotential compartments made up of their own electrical circuits (Rall, 1959). Therefore, even in the simplest case where a neuron has only uniform passive properties in every compartment (i.e., a capacitance and leakage conductance), including the soma and the axon, the axial resistance due to the resistivity of the cytoplasm will attenuate current passing down from its injection site. The same remark holds for synaptic input currents. Figure 6 illustrates the simplest discrete form of a model that has a uniform capacitance and leakage conductance. Since we now consider the location on the cable, we use a spatial variable $x$ and partial derivatives with respect to $t$. There are two main differences with respect to Figure 3:

1. In its discrete form, the cable equation assumes each compartment has length $dx$. Hence, we need to scale physical quantities expressed per unit area in the isopotential case by the area of each compartment. For a cable of radius $r$, for example, $C={C}_{m}\times \left(2\pi r\right)\times dx$ and ${G}_{Leak}={g}_{Leak}\times \left(2\pi r\right)\times dx$ (see Figure 6).

2. A current stimulating the cell, such as ${I}_{stim}$ (or ${I}_{syn}$) will not only flow through the membrane, but also to downstream compartments, which is represented by the longitudinal current, ${I}_{L}$, computed from Ohm’s law and the axial resistivity of the compartments (${R}_{A}$ in $\text{\Omega}cm$). As the distance from the stimulation site increases, the amount of longitudinal current flow decreases due to capacitive loading and leakage across the membrane. The inter-compartment resistance is obtained from ${R}_{A}$ through the formula $R=\frac{{R}_{A}}{\pi {r}^{2}}dx$. In other words, resistance is proportional to the distance between the compartments ($dx$) and inversely proportional to the cross-sectional area of the cable ($\pi {r}^{2}$).

From Figure 6 we therefore derive for the first compartment

If we let ${v}_{n}={q}_{n}-{q}_{0}-{V}_{Leak}$, we obtain:

Defining $\tau =C/{G}_{Leak}$, which is the membrane time constant associated with the leakage conductance, and $\lambda =\sqrt{\frac{r}{2{R}_{A}{g}_{Leak}}}$, which is called the space constant, we arrive at:

Similarly, the equations for the intermediate compartments and the last one are written as:

This spatially discretized system yields the membrane potentials for the $N$ compartments by solving the system of linear differential equations:

where

In this equation, ${e}_{1}$ is a column vector that is 1 for the first element and 0 for the rest of them. We have assumed that ${I}_{stim}$ is injected in the first compartment but if it is injected in the ${k}^{th}$ compartment, ${e}_{1}$ would simply have to be replaced by ${e}_{k}$.

We now consider the continuous cable equation by substituting ${v}_{n}=v\left(\left(n-\frac{1}{2}\right)dx,t\right)$, in a cable of length $l=Ndx$, and taking the limit $dx\to 0$. Equations 12–14 then become

If the entire cable is initially at rest, then the initial condition is $v\left(x,0\right)=0$.

Active Uniform Cable

When the active HH channels described in Equations 4–6 are added to Equation 17, an active cable system with uniform channel density is formulated:

This time we consider the general situation where ${I}_{stim}$ is injected at an arbitrary position ${x}_{s}$ along the cable with an amplitude, ${I}_{0}$, during the time interval ${t}_{1}\le t\le {t}_{2}$ so that ${I}_{stim}\left(x,t\right)={I}_{0}{\chi}_{\left[{t}_{1},{t}_{2}\right]}\delta (x-{x}_{s})$, where $\chi $ is the indicator function and $\delta $ the Dirac delta function. In this case the boundary conditions are $\frac{\partial V}{\partial x}\left(0,t\right)=\frac{\partial V}{\partial x}\left(l,t\right)=0$ and

$V\left(x,0\right)={V}_{r},m\left(x,0\right)={m}_{\infty}\left({V}_{r}\right),h\left(x,0\right)={h}_{\infty}\left({V}_{r}\right),n\left(x,0\right)={n}_{\infty}\left({V}_{r}\right)$ when $V$ is initially at the resting state. The staggered Euler scheme can be used to get a numerical solution for this active uniform cable equation.

In place of the time discretization used in the isopotential active model, we use a time-space grid with step sizes $dt$ and $dx$, respectively. Let $M$ be the total number of iterations over time, i.e., $M=\frac{T}{dt},0\le t\le T$. Then,

The same procedure leading from Equation 8 to Equation 9 is applied to the gating variables, yielding

To implement the half-step backward Euler rule, we collect in a single vector all spatial compartment terms:

resulting in a system of differential equations which is the discrete version of Equation 19:

where $S$ is the matrix defined in Equation 15. Finally, the half step update, ${V}^{j}=2{V}^{j-1/2}-{V}^{j-1}$, is needed to advance the membrane potential to the next half step and obtain the numerical solution of the uniform active cable.

The Active Non-Uniform Cable

The density of channels inserted in the membrane of neurons is generally non-uniform. Changes in channel densities are modeled by varying their associated maximum conductances at different locations within the neuron. One obvious example is the spike initiation zone which typically possesses much higher densities of sodium and potassium channels and therefore has the lowest threshold for spike initiation. Another consequence of spatially varying channel densities is that the resting potential may vary spatially. This in turn implies different steady states for the gating variables of channels as a function of location and thus different dynamics in response to external stimuli. In an active cable with spatially varying densities of sodium and potassium channels, the resting membrane potential is solution of the vector equation:

where ${g}_{Leak},{\overline{g}}_{Na}$, and ${\overline{g}}_{K}$ are vectors that contain different maximum conductance values for each compartment. The integration method described in the previous paragraph (Equation 22) carries over without changes in this case.

Passive Dendritic Trees

A uniform cable is not a realistic abstraction of most neurons as their shape is usually considerably more complex than a cylinder. Hines in 1984 extended the numerical methods used to efficiently solve the cable equation to passive and active dendritic trees. His numerical method is based on the generalization of the cable equation from Equation 11 to branched dendritic trees. This generalization is straightforward, mainly requiring to take into account that this leads to multiple space constants, since each cable has potentially a different radius. Consider the simple case of a soma connected to a mother branch and two daughters illustrated in Figure 7. The longitudinal currents flowing out at the junction of the mother branch with its two daughter branches illustrates such a case, as the three cables with lengths ${l}_{1},{l}_{2}$, and ${l}_{3}$ have different radiuses, ${r}_{1},{r}_{2}$, and ${r}_{3}$. We select a fixed discretization step $dx$ so that the number of compartments for each branch is calculated as ${N}_{i}={l}_{i}/dx$ for $i$ = 1, 2, and 3.

Let ${v}_{i,j}={V}_{i,j}-{V}_{leak}$ and ${v}_{i,j}^{\prime}=\frac{\partial {v}_{i,j}}{\partial t}$ then current conservation leads to

By defining the spatial constants ${\lambda}_{i}^{2}={r}_{i}/2{R}_{A}{g}_{Leak}$, Equation 15 becomes

Similar to Equation 13,

Moreover, the soma of the neuron is modeled as a single compartment with membrane potential ${V}_{4,1}$ and surface area ${A}_{4}$. If we assume that current is injected through an electrode located at the soma,

or

where ${A}_{3}=2\pi {r}_{3}dx$.

The continuous passive dendrite equation can now be derived as for the passive cable, by taking the limit, $dx\to 0$, and substituting ${v}_{i,n}={v}_{i}((n-1/2)dx,t)$. The requirement of continuity of the potential at the junction node translates into ${v}_{1}\left({l}_{1},t\right)={v}_{2}\left({l}_{2},t\right)={v}_{3}(0,t)$ and ${v}_{3}\left({l}_{3},t\right)={v}_{4}(0,t)$. Hence Equations 25 and 27 become

Equation 29 is equivalent to the following equation, using the continuity condition for the potential,

The differential equations for the remaining compartments of the linear cables, including Equation 26, are formulated as

Under the assumption that the daughter branches have sealed ends, the boundary condition is

Active Dendritic Trees

The extension of the passive dendritic tree model to the active one proceeds in a similar manner as the extension of the passive cable to the active cable. In principle, active dendritic tree models can reproduce the detailed morphologies, voltage recordings and active channel localizations evidenced by experiments. However, the most detailed models are not always the most useful ones. Moreover, results that were found using simpler model like passive ones are not necessarily useless or wrong. For example, the more active channels we incorporate into a model, the higher the number of parameters that need to be selected, potentially leading to more than one set of parameters that can reproduce the experimental data.

In 1962, Rall used a passive dendritic model of a single motoneuron to explain the two phases of extracellular voltage transients observed during antidromic activation of the soma. Although active propagation of action potentials in motoneuron dendrites were not implemented in the model, the passive dendrites were sufficient to reproduce experimental observations, as Rall recalled in a more recent publication (Stuart et al., 2016)

Until now, we reviewed how membrane properties and dendritic structure are used to construct computational models of single neurons. In the following sections, we review models that have been used to reveal how different properties of a neuron impact its function, such as action potential initiation and synaptic integration.

Action Potential Initiation

The mechanisms of action potential initiation have implications both for the downstream propagation of spikes and the upstream integration of synaptic inputs. The axon initial segment (AIS) has been known to have the lowest threshold for action potential initiation due to its high density of voltage-gated $N{a}^{+}\phantom{\rule{0.2em}{0ex}}\left(N{a}_{V}\right)$ channels. Kole et al. (2008) showed, using a morphologically realistic model of cortical layer 5 pyramidal neurons, that high $N{a}_{V}$ channel density is required to reproduce experimental data on the rising rate of the axonal action potential, its initiation, and backpropagation. Cortical layer 5 pyramidal neurons are most likely to initiate action potentials from the distal end of the AIS rather than the proximal AIS, although both exhibit similarly high levels of $N{a}_{V}$ channels (Moore, Stockbridge, & Westerfield, 1983; Palmer & Stuart, 2006). Hu et al. (2009) found experimentally that there are two kinds of $N{a}_{V}$ channels in the AIS, low threshold $N{a}_{V}1.6$ channels and high threshold $N{a}_{V}1.2$ channels, which are partially concentrated on the distal side and the proximal side of the AIS, respectively. By compartmental modeling, they found that distal $N{a}_{V}1.6$ channels help action potential initiation while proximal $N{a}_{V}1.2$ channels participate to their backpropagation toward the soma.

Synaptic Integration

Multiple synaptic inputs need to be integrated and propagated to the spike initiation zone (SIZ) for an action potential to be fired. In this context, Katz et al. (2009) investigated the functional role of different spine density distributions of apical oblique dendrites in CA1 pyramidal neurons that can generate dendritic spikes. Their experimental observations revealed that the spine density in the branches near the origin of apical obliques was about twice as high as in terminal branches and that larger post-synaptic densities were more prevalent there as well. They constructed a model based on these experimental observations where the synaptic density and synaptic conductance both decrease from the origin of a branch toward its end. This model was compared with one having uniform synaptic density and normalized synaptic conductance per branch so that the total synaptic conductance was kept the same as in the experimentally derived model. The experimentally derived model required more inputs near oblique branch ends to trigger dendritic spikes, implying a spatially dependent normalization of synaptic efficacy in real neurons. Moreover, the observed synapse distribution increased the ability of dendritic spikes to fire an action potential. For a fixed number of active synapses on a branch, the highest depolarization occurred at the soma in the model where synaptic density decreased from proximal to distal.

The role played by the large-scale spatial distribution of inputs during synaptic integration was also investigated in layer 5 cortical pyramidal neurons by Jadi, Polsky, Schiller, and Mel (2012) and Behabadi, Polsky, Jadi, Schiller, and Mel (2012). Both studies first used active multicompartmental models that described layer 5 pyramidal neurons able to generate dendritic spikes, as observed experimentally. They then tested if the results were due to interactions of synaptic inputs and active membrane properties by comparison with a two-compartment passive model representing soma and dendrites with only a leakage conductance in each compartment and synaptic currents localized in the dendritic compartment.

Jadi et al. (2012) found that, on the one hand, dendritic inhibition at the same distal location as excitation resulted in a total leak and inhibitory conductance that increased the dendritic spike threshold whereas somatic inhibition did so only slightly. On the other hand, somatic inhibition caused higher voltage attenuation (${V}_{dend}/{V}_{soma}$) whereas dendritic inhibition did not, because of the relationship ${V}_{soma}={V}_{dend}-{I}_{dend}/{g}_{A}$, where ${g}_{A}$ is the axial conductance.

Behabadi et al. (2012) studied how neocortical pyramidal neurons integrate two different classes of excitatory synaptic inputs: “driver inputs” originating from other cortical layers and associated with the basic visual receptive field properties of the neurons, and “contextual inputs” modulating the neural responses, as may be expected, for example, from changes in attention. They tested how the resulting excitatory synaptic interactions affected somatic depolarization depending on the distance between the two inputs, and their distance from the soma. When two inputs were co-localized, the relationship between the number of synapses and the peak somatic response had a consistent sigmoidal shape and as the inputs grew closer to the soma, the threshold and the amplitude of dendritic spikes increased. When two inputs were farther away from each other, the pattern of responses was more complex and asymmetric. For both subthreshold and suprathreshold responses, a proximal modulatory input paired with a distal driver lowered the threshold for spikes and increased the input and output response for the driver input. Distal modulation paired with proximal drive, however, only lowered the threshold. In both cases, similar results were obtained in the detailed and the two-compartment model.

Model Reduction and Alternative Approaches

The more features are added to a model, the more experimental observations the model can capture and explain. However, a more detailed model is also more expensive to simulate. Additionally, the problem of finding adequate parameters may become ill-posed as many combinations of parameter values could produce outcomes resembling those observed experimentally. In such cases, it may be hard to identify the real qualitative mechanism behind an experimental observation. Reduced models have thus been introduced and studied in this context. A good reduced model has a simpler structure than a realistic one, allowing one to understand easily its features, yet it should still capture sufficiently well the intrinsic properties of the subject neuron to have explanatory power. At present, there is no general method to reduce model complexity, due to the highly non-linear structure of active compartmental models. One successful example of a reduced model is that of CA3 hippocampal pyramidal neurons by Pinsky and Rinzel (1994), based on the detailed model originally introduced by Traub, Wong, Miles, and Michelson (1991).

Traub et al. developed a model of CA3 hippocampal pyramidal neurons of guinea pig based on an active cable model with nineteen compartments representing soma, basal dendrites, and apical dendrites (one, eight, and ten compartments, respectively). The model has on the order of 100 variables. Pinsky and Rinzel adopted the active currents and their gating kinetics from Traub’s model but the number of compartments was reduced to two: a “soma-like compartment” (SOMA) and a “dendrite-like compartment” (DENDRITE). Moreover, instead of having most channels in every compartment as in the Traub model shown in Figure 8A, they segregated the channels between SOMA and DENDRITE as displayed in Figure 8B. As a result, the number of variables was reduced to eight. In addition to the leakage current in each compartment, the SOMA has a sodium current, ${I}_{Na}$, and a delayed-rectifier potassium current, ${I}_{K-DR}$; the DENDRITE has a calcium current, ${I}_{Ca}$, a calcium-activated potassium current (fast activating), ${I}_{K-C}$, and a potassium afterhyperpolarization current (slow activating), ${I}_{K-AHP}$, which is not voltage-dependent but $C{a}^{+}$-dependent. This setup allows the SOMA to fire action potentials and the DENDRITE to generate calcium spikes on a much slower time scale. Finally, instead of the axial resistivity to current flow between compartments used in regular cable equations, “the coupling conductance” value was designed to reduce the attenuation of somatic spikes allowing them to initiate dendritic voltage spikes. This was achieved by adjusting the compartment surface area so that each compartment represented half of the total area, which is not based on the geometry of CA3 pyramidal neurons but was determined empirically to result in an adequate total influence of currents from each compartment for achieving the model’s objectives. Thus, the reduced model could produce simulation results qualitatively equivalent to those of Traub’s model in terms of modes of firing observed in the model and the frequency of firing as a function of injected current. Additionally, the model provided a biophysical explanation of how complex interactions between active conductances in the two compartments generate firing in CA3 pyramidal neurons. Specifically, for low somatic current injections the model generates low frequency bursts whereas for high somatic injections the model exhibits regular spiking after an initial transient period, exactly as in the Traub model.

Future Direction of Single Neuron Computation Modeling

Single neuron modeling has been developed based on anatomical and physiological observations. For example, advanced techniques including intracellular recordings, anatomical stainings, and high-resolution microscopy allow to measure the impact of channel density changes over the spatial structure of a neuron and thus help construct more realistic non-uniform cable models (Migliore & Shepherd, 2002). Moreover, imaging techniques such as two-photon microscopy, have helped to constrain in models the role played by morphological structures such as dendrites in the processing of sensory synaptic inputs (Zhu & Gabbiani, 2016).

As refined observations are made based on improving experimental techniques, simplifying assumptions often have to be reconsidered to help design more realistic models. For example, traditional single neuron compartmental models suppose that the extracellular media is isopotential, whereas in most cases it has been shown to acts as an Ohmic resistor (Miceli, Ness, Einevoll, & Schubert, 2017). As noticed early on, this assumption has little effect on single neuron modeling (Rall, 1959, 1962). However, when considering tightly packed neuronal processes in vivo, extracellular electric fields may induce coupling and synchronization of neuronal populations and thus affect the intracellular membrane potential (Anastassiou, Perin, Markram, & Koch, 2011; Fröhlich & McCormick, 2010; Jefferys, 1995). In other cases, particularly when tight, the extracellular space might not even act as an Ohmic resistor (Weckström & Laughlin, 2010), which would require generalizing cable theory to include its effects (Bédard & Destexhe, 2013). Finally, much data used to constrain single neuron models comes from experiments carried out in vitro. The extent to which these neuronal properties might differ from those observed in vivo remains unknown in many cases (Sarid, Randy, Sakmann, Segev, & Fedmeyer, 2007).

The availability of increasingly powerful computational tools pushes back the limits of compartmental modeling, as ever more detailed descriptions of neurons require increased computational power. For example, the Blue Brain Project aims to construct a biologically detailed model of an entire rodent’s brain, and eventually a model of the human brain. The model is built by incorporating large amounts of experimental data simulated using supercomputers. The first main report from this project focused on “the microcircuitry of the somatosensory cortex of juvenile rat” simulating 31,000 neurons of 207 subtypes based on a detailed morphological classification and patch clamp recordings (Markram et al., 2015). This is an example where exhaustive single neuron models are used as building blocks toward a comprehensive representation of the whole brain.

Although advanced tools can be used to build ever more realistic and complex models that better predict experimental data, they may not be sufficient to provide deeper insights about how neurons and neural networks function. Work in the opposite direction, to construct simpler models embodies the process of picking out key factors that are essential to explain a complex system. The Traub model and its reduced model by Pinsky and Rinzel illustrate this point: the full model more closely resembled the experimental observations and was useful for predicting complex network dynamics, whereas the reduced model provided a mechanistic explanation of the complex interaction between different channels in soma and dendrites underlying bursting. In the future, single neuron modeling will need to be developed in both directions to help explain how neurons and neural networks process information.

## References

Anastassiou, C. A., Perin, R., Markram, H., & Koch, C. (2011). Ephaptic coupling of cortical neurons. *Nature Neuroscience*, *14*, 217–223.Find this resource:

Bédard, C., & Destexhe, A. (2013). Generalized cable theory for neurons in complex and heterogeneous media. *Physical Review E*, *88*, 022709.Find this resource:

Behabadi, B. F., Polsky, A., Jadi, M., Schiller, J., & Mel, B. W. (2012). Location-dependent excitatory synaptic interactions in pyramidal neuron dendrites. *PLoS Computational Biology*, *8*(7), e1002599.Find this resource:

Connor, J. A., & Stevens, C. F. (1971). Prediction of repetitive firing behaviour from voltage clamp data on an isolated neurone soma. *The Journal of Physiology*, *213*(1), 31.Find this resource:

Fröhlich, F., & McCormick, D. A. (2010). Endogenous electric fields may guide neocortical network activity. *Neuron*, *67*(1), 129–143.Find this resource:

Gabbiani, F., & Cox, S. (2017). *Mathematics for neuroscientists* (2nd ed.). San Diego, CA: Academic Press.Find this resource:

Hines, M. (1984). Efficient computation of branched nerve equations. *International Journal of Bio-Medical Computing*, *15*(1), 69–76.Find this resource:

Hodgkin, A. L., & Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. *The Journal of Physiology*, *117*(4), 500.Find this resource:

Hu, W., Tian, C., Li, T., Yang, M., Hou, H., & Shu, Y. (2009). Distinct contributions of Nav1.6 and Nav1.2 in action potential initiation and backpropagation. *Nature Neuroscience*, *12*(8), 996–1002.Find this resource:

Jadi, M., Polsky, A., Schiller, J., & Mel, B. W. (2012). Location-dependent effects of inhibition on local spiking in pyramidal neuron dendrites. PLoS Computational Biology, 8(6), e 1002550.Find this resource:

Jefferys, J. G. (1995). Nonsynaptic modulation of neuronal activity in the brain: Electric currents and extracellular ions. *Physiological Reviews*, *75*(4), 689–723.Find this resource:

Katz, Y., Menon, V., Nicholson, D. A., Geinisman, Y., Kath, W. L., & Spruston, N. (2009). Synapse distribution suggests a two-stage model of dendritic integration in CA1 pyramidal neurons. *Neuron*, *63*(2), 171–177.Find this resource:

Kole, M. H., Ilschner, S. U., Kampa, B. M., Williams, S. R., Ruben, P. C., & Stuart, G. J. (2008). Action potential generation requires a high sodium channel density in the axon initial segment. *Nature Neuroscience*, *11*, 178–186.Find this resource:

Markram, H., Muller, E., Ramaswamy, S., Reimann, M., Abdellah, M., Sanchez, C., . . . Schürmann, F. (2015). Reconstruction and simulation of neocortical microcircuitry. *Cell*, *163*(2), 456–492.Find this resource:

Miceli, S., Ness, T. V., Einevoll, G. T., & Schubert, D. (2017). Impedance Spectrum in cortical tissue: Implications for propagation of LFP signals on the microscopic level. *eNeuro*, *4*(1), 0291.Find this resource:

Migliore, M., & Shepherd, G. M. (2002). Emerging rules for the distributions of active dendritic conductances. *Nature Reviews Neuroscience*, *3*, 362–370.Find this resource:

Moore, J. W., Stockbridge, N., & Westerfield, M. (1983). On the site of impulse initiation in a neurone. *The Journal of Physiology*, *336*, 301–311.Find this resource:

Palmer, L. M., & Stuart, G. J. (2006). Site of action potential initiation in layer 5 pyramidal neurons. *The Journal of Neuroscience*, *26*(6), 1854–1863.Find this resource:

Pinsky, P. F., & Rinzel, J. (1994). Intrinsic and network rhythmogenesis in a reduced Traub model for CA3 neurons. *Journal of Computational Neuroscience*, *1*(1), 39–60.Find this resource:

Rall, W. (1959). Branching dendritic trees and motoneuron membrane resistivity. *Experimental Neurology*, *1*(5), 491–527.Find this resource:

Rall, W. (1962). Electrophysiology of a dendritic neuron model. *Biophysical Journal*, *2*, 145–167.Find this resource:

Sarid, L., Randy, B., Sakmann, B., Segev, I., & Fedmeyer, D. (2007). Modeling a layer 4-to-layer 2/3 module of a single column in rat neocortex: Interweaving in vitro and in vivo experimental observations. *PNAS*, *104*(41), 16353–16358.Find this resource:

Stuart, G., Spruston, N., & Häusser, M. (2016). *Dendrites* (3d ed.). Oxford, U.K.: Oxford University Press.Find this resource:

Traub, R. D., Wong, R. K., Miles, R., & Michelson, H. (1991). A model of a CA3 hippocampal pyramidal neuron incorporating voltage-clamp data on intrinsic conductances. *Journal of Neurophysiology*, *66*, 635–650.Find this resource:

Wang, X. J. (1994). Multiple dynamical modes of thalamic relay neurons: Rhythmic bursting and intermittent phase-locking. *Neuroscience*, *59*(1), 21–31.Find this resource:

Weckström, M., & Laughlin, S. (2010). Extracellular potentials modify the transfer of information at photoreceptor output synapses in the blowfly compound eye. *The Journal of Neuroscience*, *30*(28), 9557–9566.Find this resource:

Zhu, Y., & Gabbiani, F. (2016). Fine and distributed subcellular retinotopy of excitatory inputs to the dendritic tree of a collision-detecting neuron. *Journal of Neurophysiology*, *115*(6), 3101–3112.Find this resource: