Theta model

It is possible to describe a multitude of parabolic bursting cells by deriving a simple mathematical model, called a canonical model. Derivation of the Ermentrout and Kopell canonical model begins with the general form for parabolic bursting, and notation will be fixed to clarify the discussion. The letters ${\displaystyle f}$, ${\displaystyle g}$, ${\displaystyle h}$, ${\displaystyle I}$ are reserved for functions; ${\displaystyle x}$, ${\displaystyle y}$, ${\displaystyle \theta }$ for state variables; ${\displaystyle \varepsilon }$, ${\displaystyle p}$, and ${\displaystyle q}$ for scalars.

In the following generalized system of equations for parabolic bursting, the values of ${\displaystyle f}$ describe the membrane potential and ion channels, typical of many conductance-based biological neuron models. Slow oscillations are controlled by ${\displaystyle h}$, and ultimately described by ${\displaystyle y}$. These slow oscillations can be, for example, slow fluctuations in calcium concentration inside a cell. The function ${\displaystyle g}$ couples ${\displaystyle {\dot {y}}}$ to ${\displaystyle {\dot {x}}}$, thereby allowing the second system, ${\displaystyle {\dot {y}}}$, to influence the behavior of the first system, ${\displaystyle {\dot {x}}}$. In more succinct terms, "${\displaystyle x}$ generates the spikes and ${\displaystyle y}$ generates the slow waves".[7] The equations are:

${\displaystyle {\dot {x}}=f(x)+\varepsilon ^{2}g(x,y,\varepsilon ),}$

${\displaystyle {\dot {y}}=\varepsilon h(x,y,\varepsilon ),}$

where ${\displaystyle x}$ is a vector with ${\displaystyle p}$ entries (i.e. ${\displaystyle x\in \mathbb {R} ^{p}}$), ${\displaystyle y}$ is a vector with ${\displaystyle q}$ entries (i.e. ${\displaystyle y\in \mathbb {R} ^{q}}$), ${\displaystyle \varepsilon }$ is small and positive, and ${\displaystyle f}$, ${\displaystyle g}$, ${\displaystyle h}$ are smooth (i.e. infinitely differentiable).[7] Additional constraints are required to guarantee parabolic bursting. First, ${\displaystyle {\dot {x}}=f(x)}$ must produce a circle in phase space that is invariant, meaning it does not change under certain transformations. This circle must also be attracting in ${\displaystyle \mathbb {R} ^{2}}$ with a critical point located at ${\displaystyle x=0}$. The second criterion requires that when ${\displaystyle {\dot {y}}=h(0,y,0)}$, there exists a stable limit cycle solution. These criteria can be summarized by the following points:

1. When ${\displaystyle \varepsilon =0}$, ${\displaystyle {\dot {x}}=f(x)}$ "has an attracting invariant circle with a single critical point", with the critical point located at ${\displaystyle x=0}$, and
2. When ${\displaystyle x=0}$, ${\displaystyle {\dot {y}}=h(0,y,0)}$ has a stable limit cycle solution.[7]

The theta model can be used in place of any parabolic bursting model that satisfies the assumptions above.

The theta model is a reduction of the generalized system from the previous section and takes the form,

${\displaystyle {\frac {d\theta }{dt}}=1-\cos \theta +(1+\cos \theta )I(t),\qquad \theta \in S^{1}.}$

This model is one of the simplest excitable neuron models.[18] The state variable ${\displaystyle \theta }$ represents the angle in radians, and the input function, ${\displaystyle I(t)}$, is typically chosen to be periodic. Whenever ${\displaystyle \theta }$ reaches the value ${\displaystyle \theta =\pi }$, the model is said to produce a spike.[8][18]

The theta model is capable of a single saddle-node bifurcation and can be shown to be the "normal form for the saddle-node on a limit cycle bifurcation."[8] When ${\displaystyle (I<0)}$, the system is excitable, i.e., given an appriate perturbation the system will produce a spike. Incidentally, when viewed in the plane (${\displaystyle \mathbb {R} ^{2}}$), the unstable critical point is actually a saddle point because ${\displaystyle S^{1}}$ is attracting in ${\displaystyle \mathbb {R} ^{2}}$. When ${\displaystyle (I>0)}$, ${\displaystyle {\dot {\theta }}}$ is also positive, and the system will give rise to a limit cycle. Therefore, the bifurcation point is located at ${\displaystyle I(t)=0}$.

Near the bifurcation point, the theta model resembles the quadratic integrate and fire model:

${\displaystyle {\frac {dx}{dt}}=x^{2}+I.}$

For I > 0, the solutions of this equation blow up in finite time. By resetting the trajectory ${\displaystyle x(t)}$ to ${\displaystyle -\infty }$ when it reaches ${\displaystyle +\infty }$, the total period is then

${\displaystyle T={\frac {\pi }{\sqrt {I}}}.}$

Therefore, the period diverges as ${\displaystyle I\rightarrow 0^{+}}$ and the frequency converges to zero.[8]

When ${\displaystyle I(t)}$ is some slow wave which can be both negative and positive, the system is capable of producing parabolic bursts. Consider the simple example ${\displaystyle I(t):=\sin(\alpha t)}$, where ${\displaystyle \alpha }$ is relatively small. Then for ${\displaystyle \alpha t\in (0,\pi )}$, ${\displaystyle I(t)}$ is strictly positive and ${\displaystyle \theta }$ makes multiple passes through the angle ${\displaystyle \pi }$, resulting in multiple bursts. Note that whenever ${\displaystyle \alpha t}$ is near zero or ${\displaystyle \pi }$, the theta neuron will spike at relatively a low frequency, and whenever ${\displaystyle \alpha t}$ is near ${\displaystyle \alpha t=\pi /2}$ the neuron will spike with very high frequency. When ${\displaystyle \alpha t=\pi }$, the frequency of spikes is zero since the period is infinite since ${\displaystyle \theta }$ can no longer pass through ${\displaystyle \theta =\pi }$. Finally, for ${\displaystyle \alpha t\in (\pi ,2\pi )}$, the neuron is excitable and will no longer burst. This qualitative description highlights the characteristics that make the theta model a parabolic bursting model. Not only does the model have periods of quiescence between bursts which are modulated by a slow wave, but the frequency of spikes at the beginning and end of each burst is high relative to the frequency at the middle of the burst.

The derivation comes in the form of two lemmas in Ermentrout and Kopell (1986). Lemma 1, in summary, states that when viewing the general equations above in a subset ${\displaystyle S^{1}\times \mathbb {R} ^{2}}$, the equations take the form:

${\displaystyle {\dot {x_{1}}}={\overline {f}}(x_{1})+\varepsilon ^{2}{\overline {g}}(x_{1},y,\varepsilon )\qquad x_{1}\in S^{1},}$

${\displaystyle {\dot {y}}=\varepsilon {\overline {h}}(x_{1},y,\varepsilon )\;\;\;\;\;\;y\in \mathbb {R} ^{q}.}$

By lemma 2 in Ermentrout and Kopell 1986, "There exists a change of coordinates... and a constant, c, such that in new coordinates, the two equations above converge pointwise as ${\displaystyle \varepsilon \rightarrow 0}$ to the equations

${\displaystyle {\dot {\theta }}=(1-\cos \theta )+(1+\cos \theta ){\overline {g}}(0,y,0),}$

${\displaystyle {\dot {y}}={\frac {1}{c}}{\overline {h}}(0,y,0),}$

for all ${\displaystyle \theta \neq \pi }$. Convergence is uniform except near ${\displaystyle \theta =\pi }$." (Ermentrout and Kopell, 1986). By letting ${\displaystyle I(t):={\overline {g}}(0,y,0)}$, resemblance to the theta model is obvious.

The phase response curve of the theta model with K = 1. Since perturbations always result in a phase advance, this is a type 1 PRC.

In general, given a scalar phase model of the form

${\displaystyle {\dot {\theta }}=f(\theta )+g(\theta )S(t),}$

where ${\displaystyle S(t)}$ represents the perturbation current, a closed form solution of the phase response curve (PRC) does not exist.

However, the theta model is a special case of such an oscillator and happens to have a closed-form solution for the PRC. The theta model is recovered by defining ${\displaystyle f}$ and ${\displaystyle g}$ as

${\displaystyle f(\theta )=(1-\cos \theta )+I(1+\cos \theta ),}$

${\displaystyle g(\theta )=(1+\cos \theta ).}$

In the appendix of Ermentrout 1996, the PRC is shown to be ${\displaystyle Z(\theta )=K(1+\cos \theta )}$.[40]

The authors of Soto-Treviño et al. (1996) discuss in great detail the similarities between Plant's (1976) model and the theta model. At first glance, the mechanisms of bursting in both systems are very different: In Plant's model, there are two slow oscillations – one for conductance of a specific current and one for the concentration of calcium. The calcium oscillations are active only when the membrane potential is capable of oscillating. This contrasts heavily against the theta model in which one slow wave modulates the burst of the neuron and the slow wave has no dependence upon the bursts. Despite these differences, the theta model is shown to be similar to Plant's (1976) model by a series of coordinate transformations. In the process, Soto-Trevino, et al. discovered that the theta model was more general than originally believed.

The quadratic integrate-and-fire (QIF) model was created by Latham et al. in 2000 to explore the many questions related to networks of neurons with low firing rates.[12] It was unclear to Latham et al. why networks of neurons with "standard" parameters were unable to generate sustained low frequency firing rates, while networks with low firing rates were often seen in biological systems.

According to Gerstner and Kistler (2002), the quadratic integrate-and-fire (QIF) model is given by the following differential equation:

${\displaystyle \tau {\dot {u}}=a_{0}(u-u_{\text{rest}})(u-u_{c})+R_{m}I,}$

where ${\displaystyle a_{0}}$ is a strictly positive scalar, ${\displaystyle u}$ is the membrane potential, ${\displaystyle u_{\text{rest}}}$ is the resting potential ${\displaystyle u_{c}}$ is the minimum potential necessary for the membrane to produce an action potential, ${\displaystyle R_{m}}$ is the membrane resistance, ${\displaystyle \tau }$ the membrane time constant and ${\displaystyle u_{c}>u_{\text{rest}}}$.[41] When there is no input current (i.e. ${\displaystyle I=0}$), the membrane potential quickly returns to rest following a perturbation. When the input current, ${\displaystyle I}$, is large enough, the membrane potential (${\displaystyle u}$) surpasses its firing threshold and rises rapidly (indeed, it reaches arbitrarily large values in finite time); this represents the peak of the action potential. To simulate the recovery after the action potential, the membrane voltage is then reset to a lower value ${\displaystyle u_{r}}$. To avoid dealing with arbitrarily large values in simulation, researchers will often set an upper limit on the membrane potential, above which the membrane potential will be reset; for example Latham et al. (2000) reset the voltage from +20 mV to −80 mV.[12] This voltage reset constitutes an action potential.

The theta model is very similar to the QIF model since the theta model differs from the QIF model by means of a simple coordinate transform.[10][12] By scaling the voltage appropriately and letting ${\displaystyle \Delta I}$ be the change in current from the minimum current required to elicit a spike, the QIF model can be rewritten in the form

${\displaystyle {\dot {u}}=u^{2}+\Delta I.}$

Similarly, the theta model can be rewritten as

${\displaystyle {\dot {\theta }}=1-\cos \theta +(1+\cos \theta )\Delta I.}$

The following proof will show that the QIF model becomes the theta model given an appropriate choice for the coordinate transform.

Define ${\displaystyle u(t)=\tan(\theta /2)}$. Recall that ${\displaystyle d\tan(x)/dx=1/\cos ^{2}(x)}$, so taking the derivative yields

${\displaystyle {\dot {u}}={\frac {1}{\cos ^{2}\left({\frac {\theta }{2}}\right)}}{\frac {1}{2}}{\dot {\theta }}=u^{2}+\Delta I.}$

An additional substitution and rearranging in terms of ${\displaystyle \theta }$ yields

${\displaystyle {\dot {\theta }}=2\left[\cos ^{2}\left({\frac {\theta }{2}}\right)\tan ^{2}\left({\frac {\theta }{2}}\right)+\cos ^{2}\left({\frac {\theta }{2}}\right)\Delta I\right]=2\left[\sin ^{2}\left({\frac {\theta }{2}}\right)+\cos ^{2}\left({\frac {\theta }{2}}\right)\Delta I\right].}$

Using the trigonometric identities ${\displaystyle \cos ^{2}(x/2)={\frac {1+\cos(x)}{2}}}$, ${\displaystyle \sin ^{2}(x/2)={\frac {1-\cos(x)}{2}}}$ and ${\displaystyle {\dot {\theta }}}$ as defined above, we have that

${\displaystyle {\dot {\theta }}=2\left[{\frac {1-\cos \theta }{2}}+\left({\frac {1+\cos \theta }{2}}\right)\Delta I\right]=1-\cos \theta +(1+\cos \theta )\Delta I.}$

Therefore, there exists a change of coordinates, namely ${\displaystyle u(t)=\tan(\theta /2)}$, which transforms the QIF model into the theta model. The reverse transformation also exists, and is attained by taking the inverse of the first transformation.

Though the theta model was originally used to model slow cytoplasmic oscillations that modulate fast membrane oscillations in a single cell, Ermentrout and Kopell found that the theta model could be applied just as easily to systems of two electrically coupled cells such that the slow oscillations of one cell modulates the bursts of the other.[7] Such cells serve as the central pattern generator (CPG) of the pyloric system in the lobster stomatograstic ganglion.[42] In such a system, a slow oscillator, called the anterior burster (AB) cell, modulates the bursting cell called the pyloric dilator (PD), resulting in parabolic bursts.[7]

A group led by Boergers,[16] used the theta model to explain why exposure to multiple simultaneous stimuli can reduce the response of the visual cortex below the normal response from a single (preferred) stimulus. Their computational results showed that this may happen due to strong stimulation of a large group of inhibitory neurons. This effect not only inhibits neighboring populations, but has the extra consequence of leaving the inhibitory neurons in disarray, thus increasing the effectiveness of inhibition.

Osan et al. (2002) found that in a network of theta neurons, there exist two different types of waves that propagate smoothly over the network, given a sufficiently large coupling strength.[43] Such traveling waves are of interest because they are frequently observed in pharmacologically treated brain slices, but are hard to measure in intact animals brains.[44] The authors used a network of theta models in favor of a network of leaky integrate-and-fire (LIF) models due to two primary advantages: first, the theta model is continuous, and second, the theta model retains information about "the delay between the crossing of the spiking threshold and the actual firing of an action potential". The LIF fails to satisfy both conditions.

The theta model can also be applied to research beyond the realm of biology. McKennoch et al. (2008) derived a steepest gradient descent learning rule based on theta neuron dynamics.[18] Their model is based on the assumption that "intrinsic neuron dynamics are sufficient to achieve consistent time coding, with no need to involve the precise shape of postsynaptic currents..." contrary to similar models like SpikeProp and Tempotron, which depend heavily on the shape of the postsynaptic potential (PSP). Not only could the multilayer theta network perform just about as well as Tempotron learning, but the rule trained the multilayer theta network to perform certain tasks neither SpikeProp nor Tempotron were capable of.