Modeling Interference-Free Neuron Spikes with Optogenetic Stimulation

TL;DR

This paper demonstrates how to predict and control the firing times of neurons using optogenetics and the Izhikevich model, enabling interference-free stimulation with predictable firing frequencies.

Contribution

It introduces a method to predict neuron charging and recovery times based on model parameters, improving control accuracy in optogenetic stimulation.

Findings

Charging and recovery times are predictable from model parameters.

Predictions provide lower bounds on achievable firing frequencies.

Simple functions can estimate neuron response times despite parameter deviations.

Abstract

This paper predicts the ability to externally control the firing times of a cortical neuron whose behavior follows the Izhikevich neuron model. The Izhikevich neuron model provides an efficient and biologically plausible method to track a cortical neuron's membrane potential and its firing times. The external control is a simple optogenetic model represented by an illumination source that stimulates a saturating or decaying membrane current. This paper considers firing frequencies that are sufficiently low for the membrane potential to return to its resting potential after it fires. The time required for the neuron to charge and for the neuron to recover to the resting potential are numerically fitted to functions of the Izhikevich neuron model parameters and the peak input current. Results show that simple functions of the model parameters and maximum input current can be used to predict the charging and recovery times, even when there are deviations in the actual parameter values. Furthermore, the predictions lead to lower bounds on the firing frequency that can be achieved without significant distortion.

Tables

Table 1. TABLE I: Selection of Nominal Parameter Values for the Izhikevich Neuron Model (from [ 36 ] )

Neuron Type (Acronym)	$a$	$b$	$c$	$d$
Regular Spiking (RS)	0.02	0.2	-65	8
Fast Spiking (FS)	0.1	0.2	-65	2
Low-Threshold Spiking (LTS)	0.02	0.25	-65	2
Chattering (CH)	0.02	0.2	-50	2
Intrinsically Bursting (IB)	0.02	0.2	-55	4

Table 2. TABLE II: Fitting charging behavior to a single parameter. Default parameter values are from Table  I and I max = 6 subscript 𝐼 max 6 I_{\mathrm{max}}=6 . Increment over the parameter range is 0.005 for a 𝑎 a , 0.005 for b 𝑏 b , and 0.5 for I max subscript 𝐼 max I_{\mathrm{max}} .

Neuron Type	Parameter Range	Fit	Fitting Function	$R^2$	RMSE  [ms]	Max Error [ms]
RS	$a\in [0.02,0.1]$	poly1	$4.003a+7.834$	$0.9999$	$3.72\times {10}^{-4}$	$9.61\times {10}^{-4}$
	$b\in [0.2,0.25]$	power1	$0.2922b^{-2.050}$	$0.9997$	$1.545\times {10}^{-2}$	$3.769\times {10}^{-2}$
	$I_{\max }\in [4,12]$	power2	$69.28I_{\max }^{-1.512}+3.317$	$0.9995$	$4.584\times {10}^{-2}$	$8.790\times {10}^{-2}$
FS	$a\in [0.02,0.1]$	poly1	$4.003a+7.834$	$0.9999$	$3.72\times {10}^{-4}$	$9.61\times {10}^{-4}$
	$b\in [0.2,0.25]$	power1	$0.2505b^{-2.169}$	$0.9998$	$1.246\times {10}^{-2}$	$2.632\times {10}^{-2}$
	$I_{\max }\in [4,12]$	power2	$150.6I_{\max }^{-1.993}+4.002$	$0.9980$	$10.67\times {10}^{-2}$	$0.2001$
LTS	$a\in [0.02,0.1]$	poly1	$0.8328a+4.958$	$0.9998$	$2.94\times {10}^{-4}$	$5.1\times {10}^{-4}$
	$b\in [0.2,0.25]$	power1	${0.2922}^{-2.050}$	$0.9997$	$1.549\times {10}^{-2}$	$3.769\times {10}^{-2}$
	$I_{\max }\in [2,12]$	power1	$11.37I_{\max }^{-0.4602}$	$0.9999$	$7.259\times {10}^{-3}$	$1.507\times {10}^{-2}$
IB	$a\in [0.02,0.1]$	poly1	$4.003a+7.834$	$0.9999$	$3.72\times {10}^{-4}$	$9.61\times {10}^{-4}$
	$b\in [0.18,0.21]$	power1	$0.2256I^{-2.212}$	$0.9989$	$3.175\times {10}^{-2}$	$4.816\times {10}^{-2}$
	$I_{\max }\in [4,12]$	power1	$46.60I_{\max }^{-0.9916}$	$0.9951$	$1.068\times {10}^{-1}$	$0.1376$

Table 3. TABLE III: Fitting recovery behavior to a single parameter. Default parameter values are from Table  I and I max = 6 subscript 𝐼 max 6 I_{\mathrm{max}}=6 . Increment over the parameter range is 0.005 for a 𝑎 a , 0.005 for b 𝑏 b , 1 for c 𝑐 c , 0.5 for d 𝑑 d , and 0.5 for I max subscript 𝐼 max I_{\mathrm{max}} .

Neuron Type	Parameter Range	Fit	Fitting Function	$R^2$	RMSE  [ms]	Max Error [ms]
RS	$a\in [0.02,0.1]$	power1	$3.371a^{-0.9587}$	$1.000$	$0.2696$	$0.5037$
	$b\in [0.2,0.25]$	poly3	$-1.831\times {10}^5{b}^3+1.161\times {10}^5{b}^2-2.451\times {10}^{4}b+1869$	$0.9976$	$0.1646$	$0.2718$
	$c\in [-65,-50]$	exp2	$147.0e^{3.341\times {10}^{-4}c}+9.6847\times {10}^{13}{e}^{0.6390c}$	$0.9971$	$2.799\times {10}^{-2}$	$8.074\times {10}^{-2}$
	$d\in [2,10]$	power1	$78.30b^{0.2945}$	$0.9994$	$0.3756$	$0.7580$
	$I_{\max }\in [4,12]$	power2	$23.14I_{\max }^{-2.4628}+143.6$	$0.9992$	$5.703\times {10}^{-3}$	$1.352\times {10}^{-2}$
FS	$a\in [0.02,0.1]$	power1	$3.367a^{-0.8522}$	$0.9996$	$0.4221$	$0.8628$
	$b\in [0.2,0.25]$	poly3	$-1.233\times {10}^4{b}^3+7614b^2-1548b+128.2$	$0.9998$	$2.544\times {10}^{-3}$	$3.650\times {10}^{-3}$
	$c\in [-65,-55]$	poly4	$1.559\times {10}^{-3}{c}^4+0.3819c^3$ $+35.07c^2+1431c+2.193\times {10}^4$	$0.9973$	$5.336\times {10}^{-2}$	$0.1212$
	$d\in [2,10]$	power1	$21.54d^{0.1736}$	$0.9966$	$0.1170$	$0.2609$
	$I_{\max }\in [4,12]$	exp2	$8.572e^{-0.4449I_{\max }}+22.75e^{8.625\times {10}^{-3}{I}_{\max }}$	$0.9961$	$1.418\times {10}^{-2}$	$2.463\times {10}^{-2}$
LTS	$a\in [0.02,0.1]$	power1	$3.298a^{-0.8500}$	$0.9986$	$0.7094$	$1.338$
	$b\in [0.2,0.25]$	exp2	$-8.864\times {10}^{-10}{e}^{92.44b}+71.04e^{1.1474b}$	$0.9999$	$1.613\times {10}^{-2}$	$3.052\times {10}^{-2}$
	$c\in [-70,-61]$	poly3	$1.030\times {10}^{-2}{c}^3+2.079c^2+140.3c+3252$	$0.9958$	$8.735\times {10}^{-2}$	$0.1786$
	$d\in [2,10]$	power1	$77.51d^{0.2603}$	$1.000$	$7.412\times {10}^{-2}$	$0.1997$
	$I_{\max }\in [4,10]$	poly3	$4.151\times {10}^{-2}{I}_{\max }^3-0.7104I_{\max }^2+4.274I_{\max }+84.02$	$0.9928$	$0.1189$	$0.2278$
IB	$a\in [0.02,0.1]$	power1	$3.478a^{-0.9040}$	$0.9997$	$0.4039$	$0.7936$
	$b\in [0.18,0.21]$	exp1	$97.54e^{1.048b}$	$0.9989$	$4.104\times {10}^{-2}$	$7.065\times {10}^{-2}$
	$c\in [-65,-55]$	poly3	$2.667\times {10}^{-3}{c}^3+0.4995c^2+31.27c+773.0$	$0.9981$	$2.675\times {10}^{-2}$	$5.120\times {10}^{-2}$
	$d\in [3.5,8]$	power1	$83.66d^{0.2623}$	$0.9998$	$0.1395$	$0.3388$
	$I_{\max }\in [4,6.5]$	poly2	$0.2143I_{\max }^2-2.190I_{\max }+125.7$	$0.9986$	$5.295\times {10}^{-3}$	$7.629\times {10}^{-3}$

Table 4. TABLE IV: Fitting charging behavior to b 𝑏 b and I max subscript 𝐼 max I_{\mathrm{max}} . Default parameter values are from Table  I . Increment over the parameter range is 0.005 for b 𝑏 b and 0.5 for I max subscript 𝐼 max I_{\mathrm{max}} .

Neuron Type	Parameter Range	Fitting Function	$R^2$	RMSE  [ms]	Max Error [ms]
RS	$b\in [0.2,0.25]$ $I_{\max }\in [4,12]$	$187.3-17.58I-1553b+0.6887I^2+90.70I_{\max }b+4686b^2$ $-9.523\times {10}^{-3}{I}_{\max }^3-1.765I_{\max }^{2}b-118.5I_{\max }{b}^2-5115b^3$	$0.9962$	$9.117\times {10}^{-2}$	$0.6975$
FS	$b\in [0.2,0.25]$ $I_{\max }\in [4,12]$	$269.0-24.42I_{\max }-2353b+0.9272I_{\max }^2+131.7I_{\max }b+7360b^2$ $-1.265\times {10}^{-2}{I}_{\max }^3-2.437I_{\max }^{2}b-181.2I_{\max }{b}^2-8159b^3$	$0.9913$	$0.1546$	$1.329$
LTS	$b\in [0.2,0.25]$ $I_{\max }\in [4,10]$	$228.1-23.99I_{\max }-1922b+1.033I_{\max }^2+127.3I_{\max }b+5820b^2$ $-1.577\times {10}^{-2}{I}_{\max }^3-2.719I_{\max }^{2}b-171.5I_{\max }{b}^2-6306b^3$	$0.9975$	$7.345\times {10}^{-2}$	$0.5075$
IB	$b\in [0.18,0.21]$ $I_{\max }\in [4,6.5]$	$3.900\times {10}^4-6257I_{\max }-6.137\times {10}^{5}b+477.6I_{\max }^2+6.808\times {10}^4{I}_{\max }b$ $+3.702\times {10}^6{b}^2-18.93I^3-3233I_{\max }^{2}b-2.533\times {10}^5{I}_{\max }{b}^2-1.011\times {10}^7{b}^3$ $+0.3230I_{\max }^4+59.29I^{3}b+5665I_{\max }^2{b}^2+3.205\times {10}^5{I}_{\max }{b}^3+1.052\times {10}^7{b}^4$	$0.9933$	$0.2793$	$0.8903$

Table 5. TABLE V: Fitting recovery behavior to a ∈ [ 0.02 , 0.1 ] 𝑎 0.02 0.1 a\in[0.02,0.1] and d ∈ [ 2 , 10 ] 𝑑 2 10 d\in[2,10] (except for IB neurons, where d ∈ [ 3.5 , 8 ] 𝑑 3.5 8 d\in[3.5,8] ). Default parameter values are from Table  I and I max = 6 subscript 𝐼 max 6 I_{\mathrm{max}}=6 . Increment over the parameter range is 0.005 for a 𝑎 a and 0.5 for d 𝑑 d .

Neuron Type	Fitting Function	$R^2$	RMSE  [ms]	Max Error [ms]
RS, FS	$153.5+19.16d-5491a-1.011d^2-301.7ad+7.784\times {10}^4{a}^2$ $+2.552\times {10}^{-2}{d}^3+6.646a{d}^2+1389a^{2}d-3.658\times {10}^5{a}^3$	$0.9950$	$1.924$	$5.673$
LTS	$154.8+14.51d-5358a-0.6033d^2-248.3ad+7.490\times {10}^4{a}^2$ $+9.958\times {10}^{-3}{d}^3+5.148a{d}^2+1154a^{2}d-3.484\times {10}^5{a}^3$	$0.9949$	$1.801$	$5.260$
IB	$180.5+11.60d-6048a-0.1827d^2-229.3ad+8.375\times {10}^4{a}^2$ $-6.411\times {10}^{-3}{d}^3+3.422a{d}^2+1126a^{2}d-3.892\times {10}^5{a}^3$	$0.9953$	$1.952$	$5.416$

Table 6. TABLE VI: Fitting behavior of RS neuron to maximum current I max ∈ [ 4 , 12 ] subscript 𝐼 max 4 12 I_{\mathrm{max}}\in[4,12] . Model parameter values are { a , b , c , d } = { 0.02 , 0.2 , − 65 , 8 } 𝑎 𝑏 𝑐 𝑑 0.02 0.2 65 8 \{a,b,c,d\}=\{0.02,0.2,-65,8\} . Increment of I max subscript 𝐼 max I_{\mathrm{max}} over the parameter range is 0.5.

Behavior	Fit	Fitting Function	$R^2$	RMSE [ms]	Max Error [ms]
Charging	poly1	$-0.7515I_{\max }+13.00$	$0.8704$	$0.7102$	$1.9256$
	poly2	$0.1229I_{\max }^2-2.718I_{\max }+20.13$	$0.9810$	$0.2717$	$0.6966$
	poly3	$-2.209\times {10}^{-2}{I}_{\max }^3+0.6532I_{\max }^2-6.723I_{\max }+29.54$	$0.9971$	$0.1058$	$0.2326$
	poly4	$4.133\times {10}^{-3}{I}_{\max }^4-0.1544I_{\max }^3+2.177I_{\max }^2-14.18I_{\max }+42.55$	$0.9996$	$4.123\times {10}^{-2}$	$8.338\times {10}^{-2}$
	power1	$35.70I_{\max }^{-0.8232}$	$0.9891$	$0.2064$	$0.5178$
	power2	$69.28I_{\max }^{-1.512}+3.317$	$0.9995$	$4.584\times {10}^{-2}$	$8.790\times {10}^{-2}$
	exp1	$16.97e^{-0.1161I_{\max }}$	$0.9370$	$0.4954$	$1.261$
	exp2	$69.34e^{-0.7032I_{\max }}+9.784e^{-5.901\times {10}^{-2}{I}_{\max }}$	$0.9998$	$2.455\times {10}^{-2}$	$5.011\times {10}^{-2}$
Recovery	poly1	$-7.007\times {10}^{-2}{I}_{\max }+144.4$	$0.7614$	$9.607\times {10}^{-2}$	$0.2563$
	poly2	$1.679\times {10}^{-2}{I}_{\max }^2-0.3387I_{\max }+145.4$	$0.9690$	$3.461\times {10}^{-2}$	$8.838\times {10}^{-2}$
	poly3	$-2.825\times {10}^{-3}{I}_{\max }^3+8.458\times {10}^{-2}{I}_{\max }^2-0.8506I+146.6$	$0.9955$	$1.324\times {10}^{-2}$	$2.907\times {10}^{-2}$
	poly4	$5.180\times {10}^{-4}{I}_{\max }^4-1.940\times {10}^{-2}{I}_{\max }^3+0.2756I^2-1.785I+148.2$	$0.9993$	$5.111\times {10}^{-3}$	$1.041\times {10}^{-2}$
	power1	$145.0I^{-3.878\times {10}^{-3}}$	$0.8724$	$7.027\times {10}^{-2}$	$0.1783$
	power2	$23.14I^{-2.463}+143.6$	$0.9992$	$5.703\times {10}^{-3}$	$1.352\times {10}^{-2}$
	exp1	$144.4e^{-4.875\times {10}^{-4}{I}_{\max }}$	$0.7618$	$9.599\times {10}^{-2}$	$0.2559$
	exp2	$7.081e^{-0.5910I}+143.7e^{-2.868\times {10}^{-5}{I}_{\max }}$	$0.9995$	$4.280\times {10}^{-3}$	$8.278\times {10}^{-3}$

Table 7. TABLE VII: Fitting charging behavior of RS neurons to b ∈ [ 0.2 , 0.25 ] 𝑏 0.2 0.25 b\in[0.2,0.25] and I max ∈ [ 4 , 12 ] subscript 𝐼 max 4 12 I_{\mathrm{max}}\in[4,12] . Default parameter values are from Table  I . Increment over the parameter range is 0.005 for b 𝑏 b and 0.5 for I max subscript 𝐼 max I_{\mathrm{max}} .

Fit	Fitting Function	$R^2$	RMSE  [ms]	Max Error [ms]
poly11	$20.23-0.4692I-48.34b$	$0.8687$	$0.5366$	$3.187$
poly22	$54.68-3.526I-249.5b+6.301\times {10}^{-2}{I}_{\max }^2+9.105Ib+285.1b^2$	$0.9788$	$0.2156$	$1.540$
poly33	$187.3-17.58I-1553b+0.6887I^2+90.70Ib+4686b^2$ $-9.523\times {10}^{-3}{I}_{\max }^3-1.765I^{2}b-118.5I{b}^2-{5115}^3$	$0.9962$	$9.117\times {10}^{-2}$	$0.6975$
poly44	$591.3-82.86I-6603b+4.757I^2+685.3Ib+2.838\times {10}^4{b}^2-0.1351I^3-24.82I^{2}b-1962I{b}^2$ $-5.414\times {10}^4{b}^3+1.525\times {10}^{-3}{I}_{\max }^4+0.3412I^{3}b+33.03I^2{b}^2+1948I{b}^3+3.716\times {10}^4{b}^4$	$0.9993$	$3.962\times {10}^{-2}$	$0.2983$

Equations

$$I(t)=I(t_{\mathrm{on}})+(I_{\max }-I(t_{\mathrm{on}}))(1-e^{-(t-t_{\mathrm{on}})/{\tau }_{\mathrm{on}}}),$$

$$I(t)=I(t_{\mathrm{off}})e^{-(t-t_{\mathrm{off}})/{\tau }_{\mathrm{off}}},$$

$$\frac{\mathrm{d}v(t)}{\mathrm{d}t}=$$

$$\frac{\mathrm{d}u(t)}{\mathrm{d}t}=$$

$$\text{if }v(t)\ge$$

$$v_{\mathrm{new}}=v_{\mathrm{old}}+\Delta t\frac{\mathrm{d}v(t)}{\mathrm{d}t}$$

$$v_{\mathrm{rest}}=12.5b-62.5\pm 12.5\sqrt{b^2-10b+2.6}.$$

$$R^2=1-\frac{{\sum }_{n=1}^N(y_n-{\hat{y}}_n{)}^2}{{\sum }_{n=1}^N(y_n-\overline{y}{)}^2},$$

$$\mathrm{RMSE}=\sqrt{\frac{{\sum }_{n=1}^N(y_n-{\hat{y}}_n{)}^2}{N}},$$

$$\mathrm{Max}\mathrm{Error}=\underset{n}{\max }|y_n-{\hat{y}}_n|.$$

Full text

I Introduction

Over the past decade, developments in optogenetics have given researchers the ability to directly stimulate neurons [2, 3]. Using this technique, neurons are modified with a gene that encodes a light-sensitive protein (i.e., an opsin), causing the neurons to express opsins on their surface. Certain opsins, such as channelrhodopsin [4], open an ion channel in response to light. When the channels are open, an ion current flows through the neuron’s membrane, changing its electric potential and causing it to fire. Thus, if an optogenetically-modified neuron is stimulated with a strong light source, such as a laser, then the neuron will eventually fire in response.

Dramatic advances in the study of the brain, as well as revolutionary new therapies for neurological disorders, are expected to follow from precise optogenetic control over neural circuits [5]. So far, research has often focused on the control of large groups of neurons in experimental settings [6]; e.g., studies of seizures in the mouse brain [7] or of spinal cord injury in rats [8]. However, targeted control of individual neural circuits are of considerable interest, and recent experimental results have demonstrated the feasibility of this approach [9, 10, 11]. It is widely expected that this control will one day lead to optogenetics-based therapies for neurological problems [12], such as epilepsy [13] or recovery from neural injury [14]. Within the domain of communication and networking, controlled synaptic communication has been proposed as a potential communication technique for nanonetworks [15, 16], including potential interfaces between neurons and nanomachines [17]. Communication models of synaptic systems have also been produced [18, 19, 20]. Furthermore, optogenetic techniques have been discussed in the context of therapeutic nanonetworks [21, 22] and for brain-machine interfaces [23].

For these applications, an interesting problem is to precisely control the firing time of an individual neuron, as shown conceptually in Fig. 1. Consider a neuron illuminated by a light source, where $\ell(t)$ is the time-varying light intensity. Let $\mathbf{t}=[t_{1},t_{2},\ldots,t_{n}]$ represent a vector of times at which the neuron fires. Then the neuron may be viewed as a functional $\mathsf{n}(\cdot)$, taking $\ell(t)$ as input and returning $\mathbf{t}$. The control problem is to invert $\mathsf{n}(\cdot)$: that is, given a desired vector $\mathbf{t}$, find $\ell(t)$ as a solution for $\mathbf{t}=\mathsf{n}(\ell(t))$.

The solution to this problem strongly depends on the neuron model $\mathsf{n}(\cdot)$, for which different models exist.

There have been various approaches to this problem in the recent literature. Some approaches treat the optogenetic stimulus and response as a control system [24, 25], an approach that has led to designs for therapeutic medical devices [26]. Other approaches have focused on detailed neuron models with optogenetic ion channels (particularly channelrhodopsin) [27], or models based on photoconversion [28].

The simple, yet tractable, integrate-and-fire (IF) model is an important model for neurons. It has been considered for optogenetic systems in populations of coupled neurons [29], and (in our own previous work) for individual neurons subject to a distortion criterion on the output [30, 31]. In addition, optogenetic control strategies for ensembles of neurons using this model were articulated in [32]. The IF model considers neurons as capacitors, where the current is integrated over time to find the neuron’s potential; once the potential exceeds a threshold, the neuron fires. IF is a first-order linear differential equation model, but its simplicity hides much of the complexity of real neurons. In particular, there are practical neuron behaviors that cannot be readily observed using the IF model; see [33]. Various other neuron models include linear models that address issues with IF, such as the leaky IF model, and nonlinear models, of which the Hodgkin-Huxley model [34] is likely the best known. Control strategies for neurons represented by a class of one-dimensional phase models were presented in [35], including a simplified version of the Hodgkin-Huxley model. In this paper, we use a simplified (but realistic) nonlinear model known as the Izhikevich neuron model [36] (which we hereafter simply refer to as the Izhikevich model).

This paper considers how to control the optogenetic stimulation of neurons that follow the Izhikevich model, as summarized in Fig. 2. The Izhikevich model is relatively simple to describe and simulate, but is biologically plausible because the range of neuron firing patterns that can be observed is consistent with all known types of cortical neurons, as demonstrated in [33] by tuning the model parameters. This is unlike other simple models, such as the IF model and its variants. The spiking patterns that can be generated using the Izhikevich model include the following: regular spiking (RS) neurons, in which spikes occur less frequently as stimulation is maintained; fast spiking (FS) neurons, where spiking at a high frequency can be maintained; low-threshold spiking (LTS) neurons, which are an intermediate between RS and FS; chattering (CH) neurons, in which spikes can occur in multiple bursts; and intrinsically bursting (IB), which can produce both regular spikes as well as irregular bursts. While the details of the Izhikevich model do not directly align with the biophysical mechanisms that underlie membrane potential dynamics (e.g., refractory periods are not clearly observed), its simplicity makes it amenable to the analysis that we undertake in this work. Models that include biophysical parameterization, e.g., the Hodgkin-Huxley family of models [34], can be considered in future work.

The specific contributions of this work are as follows:

We use curve fitting to estimate the illumination period required for an optogenetically-modified neuron to fire and recover to its resting potential, as a function of the Izhikevich model parameters. As in [1], our examples focus on RS neurons, but in this work we also give corresponding results for FS, LTS, and IB neurons111We note that “chattering” neurons do not align well with the methodology in this paper because, by design, they are prone to spiking multiple times even after the stimulating current is removed.. Furthermore, we adopt a more realistic optogenetic current model that is an approximation of results presented in [37], and we also fit the neuron behaviour to the peak optogenetically-induced current. Our results show that our approach leads to accurate estimation of both the charging and recovery time, as measured by metrics including the mean squared error. This enables the generation of arbitrary spike sequences when there is sufficient time between consecutive spikes. 2. 2.

We illustrate control of spike sequence generation by observing the distortion as a function of spike frequency. This expands the brief investigation of generating different spike frequencies in [1]. We show how our numerical fits enable us to predict a lower bound on the achievable frequency without significant distortion. If additional distortion can be tolerated, then our results demonstrate that we can generate spikes at a target frequency that is up to about $50\,\%$ greater than that predicted by our numerical method.

The rest of this paper is organized as follows. Section II describes the optogenetic and membrane potential models. We couple the two models in Section III. We numerically fit the times for both charging and recovery, and observe the distortion as a function of a target firing frequency, in Section IV. We conclude in Section V.

II Physical Models

In this section, we briefly describe the two physical models that we integrate to describe the neuron stimulation and membrane potential. These are the optogenetic model for the external stimulation and the Izhikevich model for the membrane potential dynamics.

II-A Optogenetic System Model

Neurons, like all animal cells, maintain an electric potential difference across their membranes. This membrane potential can be varied through the selective opening and closing of ion channels on the cell surface, allowing ions such as Na+, Ca2+, K+, and Cl2- to flow across the membrane. Neurons have voltage-gated ion channels, which open in response to changes in the membrane potential. This sets up a positive feedback loop. For example, in depolarization, a stimulus causes Na+ channels to open, thus raising the membrane potential, which causes more Na+ channels to open, further raising the membrane potential, and so on. The resulting rapid change in membrane potential causes the neuron to “fire”; see [38].

Ion channels can also be light-gated, such that they open in response to light. A well-studied example of this is channelrhodopsin (ChR); see [4, 39, 40]. An optogenetically-modified neuron expresses light-gated channels in addition to voltage-gated channels. Thus, illuminating the neuron (for example with a laser) can initiate the firing of the neuron by triggering the initial flow of ions.

While the ion channel is open, the ion current passing through the channel is dependent on a number of environmental factors, including pH and ion concentration [39]. It can also depend on the precise number and location of receptors on the surface of the neuron, which is usually unknown. Moreover, the dwell time in each channel state is a random variable. Works that model the states in detail include [37], and experimental results [4, 39] suggest that a neuron will experience a stable steady-state current in response to a constant illumination intensity $\ell(t)$. In this work, we assume a constant illumination intensity that results in an approximation of the current model in [37]. Specifically, if illumination starts at time $t=t_{\mathrm{on}}$ and remain on, then we assume that current $I(t)$ passing through the membrane at time $t>t_{\mathrm{on}}$ is of the form

[TABLE]

where $\tau_{\mathrm{on}}$ is an optogenetic time constant. Furthermore, if the illumination turns off at time $t_{\mathrm{off}}$, then the current through the membrane decays according to

[TABLE]

where $\tau_{\mathrm{off}}$ is another time constant. Using this current approximation, we find very good agreement with the membrane current dynamics shown in [37, Fig. 9] for short illumination times, i.e., on the order of less than 20 ms, when we set $\tau_{\mathrm{on}}=\tau_{\mathrm{off}}=2\,$ms (which we assume for the remainder of this work). By including $t_{\mathrm{on}}$ and $t_{\mathrm{off}}$, both (1) and (2) also readily model the current if the illumination is repeatedly turned on and off.

II-B Izhikevich Neuron Model

The Izhikevich model uses a two-dimensional system of ordinary differential equations where the variables are the membrane potential $v(t)$ and the membrane recovery variable $u(t)$. $u(t)$, which accounts for the activation of potassium ionic current and the inactivation of sodium ionic currents, provides negative feedback to $v(t)$. The system of equations was obtained via fitting to natural spike initiation dynamics of cortical neurons and is as follows [36, Eqs. (1)–(3)]:

[TABLE]

where (3) and (4) update the rates of change of $v(t)$ and $u(t)$, respectively, and (7) resets $u(t)$ and $v(t)$ after a spike occurs. Time and potential are measured in $\mathrm{ms}$ and $\mathrm{mV}$, respectively. $I(t)$ is the synaptic or input current through the ion channels in the dendrites and it is normalized. The parameters $a$, $b$, $c$, and $d$ are the fitting parameters and they can be tuned for different types of neurons; see nominal values in Table I. We emphasise that parameter values are obtained by fitting to neuron membrane dynamics and thus vary for individual neurons. $a$ sets the time scale of the decay of recovery variable $u(t)$ after a spike occurs. $b$ describes the sensitivity of $u(t)$ to subthreshold fluctuations of $v(t)$, and furthermore it can be used to define the membrane resting potential. $c$ is the reset potential for $v(t)$ after a spike occurs, and $d$ determines the reset of $u(t)$ after a spike occurs.

Results in [36, 33] demonstrate that the Izhikevich model can produce the behaviors of different types of cortical neurons by appropriately tuning the parameters $\{a,b,c,d\}$, even though the model itself is not analytically derived and so is not biophysically meaningful. Each type of neuron is associated with a characteristic firing pattern, where each firing pattern is a sequence of spikes. Nominal model parameters for a selection of neuron types that are suitable for a broad range of neural behavior are listed in Table I.

III Simulating Neuron Spikes

In this section, we present the simulation of spikes in the Izhikevich model when stimulation is provided by the simple optogenetic model. First, we describe how we directly couple the two models. Next, we demonstrate the simulation of a sequence of spikes to motivate the selection of a suitable simulation time step. We also use these simulations to motivate our interest in studying individual spikes. Finally, we assess the impact of the model’s initial conditions and derive the steady-state potentials of the Izhikevich model in the absence of an input current.

III-A Coupling the Izhikevich Model with Optogenetics

We take a direct approach to couple the two physical models for an individual neuron. We assume that the optogenetic stimulation is the membrane’s only external current source at the dendrites and it defines the input current $I(t)$ in (3), such that membrane current is bounded within $0\leq I(t)\leq I_{\mathrm{max}}$. In practice, this is an approximation, since the Izhikevich model was initially developed for natural neurons, where input currents enter via the activation of neurotransmitter receptors at the dendrites. We assume that we can control where the light-gated channels are expressed in the membrane to imitate the conditions for the Izhikevich model. Otherwise, alternative means to describe the membrane dynamics would be required, which can be considered in future work. To use the simplified optogenetic model in (1) and (2), we assume that we can turn the light source on and off as needed, and that it provides constant illumination $\ell(t)$ when the light source is on. Thus, to simulate the complete system, we only need to initialize $\{u(t),v(t),I(t)\}$ and use (3)–(7) in a loop to update $u(t)$ and $v(t)$, where we update $I(t)$ or fire the neuron when required.

III-B Choice of Time Step

We must choose a time step $\Delta t$ to set the resolution with which we evaluate (3)–(7). Specifically, $\Delta t$ is needed to update $u(t)$ and $v(t)$ from $\frac{\mathrm{d}u(t)}{\mathrm{d}t}$ and $\frac{\mathrm{d}v(t)}{\mathrm{d}t}$, respectively, i.e., we update $v(t)$ as

[TABLE]

and correspondingly update $u(t)$. In Fig. 3, we test different values of $\Delta t$ for a regular spiking neuron by setting the (normalized) input current to a constant $I(t)=I_{\mathrm{max}}=10$ (practical values for plateau currents can be on the order of 100 pA or more; see [37, 41]). The default value of $\Delta t$ in [36, 33] is $\Delta t=10^{-3}\,\textnormal{s}$, but we see in Fig. 3a) that this results in an insufficient level of granularity for our analysis, i.e., $\frac{\mathrm{d}u(t)}{\mathrm{d}t}$ and $\frac{\mathrm{d}v(t)}{\mathrm{d}t}$ change too much over the scale of $\Delta t=10^{-3}\,\textnormal{s}$ to accurately update $v(t)$ in (8). Thus, it appears that spikes are occurring before $v(t)$ reaches the spiking potential of 30 mV and furthermore that the spikes are occurring at random potentials. This can be mitigated by decreasing $\Delta t$. However, decreasing $\Delta t$ also increases the computational resources required to simulate the neuron. The timing of the spikes is indistinguishable for $\Delta t=10^{-5}\,\textnormal{s}$ and $\Delta t=10^{-6}\,\textnormal{s}$, and the membrane potential for these cases always peaks at about 30 mV for each spike, but we use $\Delta t=10^{-6}\,\textnormal{s}$ in the remainder of this work to have sufficient resolution for the numerical fits. Unless otherwise stated, we also use $I_{\mathrm{max}}=6$.

III-C Generating Multiple Spikes

From Fig. 3, we also observe that the interspike intervals are not constant, and this is independent of the choice of $\Delta t$. This behavior is expected for regular spiking neurons and other types of neurons as well. However, our objective is to fit expressions to describe a neuron’s behavior and control when it fires. As an early work in this direction, we seek to ignore the effects of interspike interference, so we focus here on predicting the generation and recovery of individual spikes, as shown in Fig. 2. We then use the results as a baseline for sequences of multiple spikes where the neuron is only stimulated while it is charging from rest. Repeated spiking patterns due to on-going input current is a scenario for future work.

III-D Initial Conditions and the Steady State

To maintain accuracy in our numerical analysis, we need to impose consistent conditions on the membrane. To generate a single spike, we will turn the illumination “on” until the neuron fires and then leave the illumination “off”. In the absence of illumination, the membrane current will decrease to 0 and the membrane potential of the neuron should eventually converge to a resting potential (unless it is bistable or inhibition induced; see [33]). We can calculate the resting potential by setting the left hand sides of (3) and (4) to [math], setting the input current $I(t)$ to [math], and then solving the two equations for $u(t)$ and $v(t)$. From (4) we can write $u(t)=bv(t)$, which we can substitute into (3) and re-arrange for $v(t)$ to show that the two possible resting potentials are

[TABLE]

The more negative solution of (9), $v_{\mathrm{rest}}^{-}$, is stable. The more positive $v_{\mathrm{rest}}^{+}$ is unstable and is in fact the spike generation threshold. In other words, if the membrane potential is above $v_{\mathrm{rest}}^{+}$, then $v(t)$ will increase and the neuron will fire even if $I(t)=0$ (though firing within this model could still be avoided with a sufficiently large negative current). Strictly speaking, we do not need to keep illuminating once the membrane potential reaches $v_{\mathrm{rest}}^{+}$, but we assume that it would be easier to detect when the peak membrane potential is reached. Furthermore, automatic firing can still be accelerated by providing an input current, thus providing more precise control over firing times and a margin of error to avoid stopping illumination too early when generating multiple spikes.

If the membrane potential is lower than $v_{\mathrm{rest}}^{+}$ and no input is applied, then the potential will converge to $v_{\mathrm{rest}}^{-}$. Throughout this work, we assume that the potential has converged once it remains within $\epsilon=0.5\%$ of $v_{\mathrm{rest}}^{-}$. We will see that this is a conservative estimate; in practice, we will not need to be so close to the resting potential before we can stimulate again without noticeable interspike interference.

We refer to the time needed for the neuron to fire as the charging time and the time to reach the stable resting potential as the recovery time. We show in Fig. 4, where $v_{\mathrm{rest}}^{-}=-70\,\mathrm{mV}$, that both of these times are sensitive to the initial membrane potential. To facilitate the application of this model to the generation of multiple spikes, we impose for the rest of this work that the initial membrane potential is also the resting potential $v_{\mathrm{rest}}^{-}$, and that the recovery variable $u(0)$ is initially $bv(0)$ (i.e., (4) is 0).

IV Numerical Fitting Results

In this section, we assess whether we can predict the timing behavior, i.e., the charging and recovery times of the Izhikevich neuron model, based on knowledge of the model parameters. Specifically, we seek numerically-derived equations for a neuron’s behavior as a function of $\{a,b,c,d,I_{\mathrm{max}}\}$. We are not predisposed towards any particular class of equations, but we seek results that are sufficiently accurate to use as a guide to control firing times and know how long to wait between firing times (i.e., for the membrane to return to the resting potential before we should start charging it again). Our assumptions limit the usefulness of very high precision; the optogenetic model is a simplifying approximation that produces results over short illumination periods that are consistent with [37], the model parameters $\{a,b,c,d\}$ cannot be directly measured, and we assume that all of the models are deterministic, i.e., there are no physical noise sources. However, the maximum current $I_{\mathrm{max}}$ can be externally controlled to some extent by modifying the illumination intensity (though it may not be constant in practice). We seek to gain intuition about controlling a neuron, and in particular we will estimate and measure the maximum firing frequency that can be achieved without interspike interference.

The remainder of this section is organized as follows. First, we measure the charging time and the recovery time as functions of the individual model parameters (including the maximum input current $I_{\mathrm{max}}$), where the remaining model parameters are fixed. This helps us decide which parameters to focus on in a joint model. For all types of neurons considered (RS, FS, LTS, and IB), the charging time is most sensitive to $b$ and $I_{\mathrm{max}}$ (we note that (3)–(7) show that charging time is independent of $c$ and $d$), and the recovery time is most sensitive to $a$ and $d$. Next, we measure the charging time as a function of both $b$ and $I_{\mathrm{max}}$ and the recovery time as a function of both $a$ and $d$. All fitting functions are found via nonlinear least squares in MATLAB using the fit function with default tolerances. Finally, we consider the stimulation of multiple spikes, where we predict the interference-free firing frequency and measure the deviations from target firing times as a function of the target firing frequency.

IV-A Fitting to a Single Spike

As we are primarily interested in the charging time and recovery time for each neuron type, we use curve fitting to develop accurate models for these properties under various parameter values. This is a challenging task given the five-dimensional parameter space. We first consider fits to individual parameters, keeping other parameters at a “typical” value for a particular neuron type, and then consider fits to multiple parameters. We give a detailed explanation and analysis of our method using a regular spiking (RS) neuron as an example; results for the other types of neurons are summarised in the corresponding tables.

We measure the accuracy of the fitting functions with three methods. $R^{2}$ measures the proportion of the variance in the behavior that is predictable from the model parameters, where $R^{2}\in[0,1]$. The root mean square error (RMSE) measures the standard deviation of the behavior from that predicted by the fitting functions. The maximum error (Max Error) is simply the absolute value of the largest deviation from the fitting function over the parameter range or ranges considered. Consider that we are fitting to a total of $N$ parameter value combinations (where we vary one or more of the parameters $\{a,b,c,d,I_{\mathrm{max}}\}$). We then suppose that $y_{n}$ is the charging time (in ms) for the $n$th combination of the model parameters, $\hat{y}_{n}$ is the corresponding estimated charging time due to some fitting function, and $\overline{y}$ is the average charging time over all $N$ model parameter combinations. A similar description can be made for recovery time. Then, $R^{2}$ for the charging time is measured as

[TABLE]

the RMSE is measured in ms as

[TABLE]

and the maximum error in ms is

[TABLE]

To fit behavior to the individual parameters, we consider polynomial functions up to degree 4 (i.e., from linear to quartic, beyond which minimal improvement was observed), exponential functions with either 1 or 2 terms, and power functions of the form $y=nx^{m}+p$. These fitting functions were the most relevant in MATLAB’s Curve Fitting Toolbox. To fit the behavior to the individual model parameters, we vary one parameter while holding the remaining parameters constant. The chosen range of each parameter is in consideration of the types of neurons listed in Table I. Using the RS neuron as an example, our default parameter values are $\{a,b,c,d\}=\{0.02,0.2,-65,8\}$, which is consistent with RS, and our default maximum current is $I_{\mathrm{max}}=6$. If we vary one of the parameters $\{a,b,c,d,I_{\mathrm{max}}\}$, then the remainder are fixed at the default value. The range of each varied parameter, a selection of fitted equations for their behavior (chosen for quality and space), and the accuracy of each fit are summarized in Table II for charging time and Table III for recovery time (see the Appendix for additional functions that fit the behavior of the RS neuron to the maximum current $I_{\mathrm{max}}$). Our considered parameter ranges varied for different neuron types in order to guarantee that the neuron would fire and would not fire more than once. We note that, as we might expect from Table I, some of the fits for different neuron types are identical because there are common parameter values and ranges. This is particularly the case for charging time because it is only a function of two of the Izhikevich model parameters. For example, since RS, FS, and IB neurons all have the same nominal value of $b$ and the same range for $a$, they also have the same fitting function of $a$ for charging time.

In Fig. 5, we plot the charging and recovery times for a single spike of a nominal RS neuron while varying one individual model parameter. A representative numerical fit accompanies each plot, and is generally chosen to be the simplest fit that results in $R^{2}>0.995$. The results are generally consistent with the other types of neurons that we consider, and are also consistent with what we would expect given (3)-(7).

The charging time in Fig. 5 depends on $\{a,b,I_{\mathrm{max}}\}$. While the charging time is nearly independent of $a$, it noticeably decreases with increasing $b$ or $I_{\mathrm{max}}$. The recovery time depends on all of the model parameters, but is nearly independent of $c$ and $I_{\mathrm{max}}$. It is not surprising for the maximum magnitude of the current to not have a significant impact on the the recovery time, since the illumination is always turned off during recovery and the current decays. However, it might be surprising that parameter $c$, which via (7) dictates the reset potential after the neuron fires, has a negligible impact on the time to recover. This is due to the exponential recovery behavior. The recovery time is most sensitive to $a$ and $d$.

Perhaps with the exception of the maximum current $I_{\mathrm{max}}$, because it is an external and controllable parameter, fitting to multiple model parameters is preferable. So, based on the single-parameter fitting for an RS neuron in Fig. 5, we consider two-parameter fits for an RS neuron. In particular, we fit to the charging time by varying $b$ and $I_{\mathrm{max}}$, and we fit to the recovery time by varying $a$ and $d$. We consider polynomial surfaces up to degree 4, where for simplicity both parameters always have the same degree. We hold the remaining model parameters constant according to the nominal parameter values in Table I. A fitted surface for each type of neuron and the accuracy of each fit are summarized in Table IV for charging time and Table V for recovery time (see the Appendix for additional functions that fit the charging time of the RS neuron to the current $I_{\mathrm{max}}$ and parameter $b$).

In Fig. 6, we plot the charging time as a function of $b$ and $I_{\mathrm{max}}$ and the recovery time as a function of $a$ and $d$ for a nominal RS neuron. We include the third order polynomial surface fit for both the charging time and the recovery time. Both surface fits agree with the numerical data, as indicated in Tables IV and V. We can see that the charging time is sensitive to both $b$ and $I_{\mathrm{max}}$ for the entire range of parameter values considered, whereas the recovery time is relatively more sensitive to $a$ than to $d$.

One might question how reliably we can depend on the particular model parameter values if the Izhikevich model itself was obtained via numerical fitting to experimental data. Since the charging time is generally much faster than the recovery time, we measure the sensitivity of the charging time to random model parameters $a$ and $b$ in Fig. 7, where we predict the charging time as a function of the maximum stimulation current $I_{\mathrm{max}}$ and we assume that the $a$ and $b$ parameters are the nominal values for an RS neuron, i.e., $\{a,b\}=\{0.02,0.2\}$. For each considered value of $I_{\mathrm{max}}$, we generate $10^{3}$ realizations of $a$ and $b$ parameters that are uniformly distributed over the ranges $a\in[0.02,0.036]$, $b\in[0.2,0.21]$, calculate the charging time from rest for each realization by solving (3)–(7), and then plot the distribution of the charging times. Fig. 7 shows that the actual charging times deviate from the predicted value by less than 10% for most of the range of maximum currents $I_{\mathrm{max}}\in[4,12]$.

IV-B Applying Fits to Spike Trains

Thus far, we have focused on the generation of individual spikes in isolation. However, we can apply our results to the generation of spike trains. In particular, we can assess how well we can generate a spike train at a target frequency, where the period is the sum of the charging and recovery times. We turn the illumination on for the expected time to charge the neuron, leave the illumination off for the neuron to recover, and then charge the neuron again. We are interested in measuring the deviation in spike times from the target frequency when we provide insufficient time for the neuron to fully recover. Similar to our work with the IF neuron model in [30, 31], we can use (11) to calculate the RMSE associated with a spike train of $N$ spikes, but where $y_{n}$ is the $n$th target firing time (according to a specified firing frequency), and $\hat{y}_{n}$ is the corresponding observed firing time.

First, we observe deviations visually. Using Tables IV and V, we expect an RS neuron with nominal model parameters and maximum current $I_{\mathrm{max}}=6$ to take $8.13$ ms to charge and $138.2$ ms to recover. Thus the interference-free firing frequency is approximately 6.8 Hz. In Fig. 8, we observe the input current and membrane potential of an RS neuron versus time as we try to generate spikes at 10 Hz and 13 Hz, where in each case we turn on the stimulating current for $8.13$ ms. At 10 Hz, we observe that the spikes can still be generated but that deviations from the target firing time are visually apparent with the third spike (since we expect the neuron to fire as soon as the current is turned off). At 13 Hz, there is a more visible deviation with the second spike and then the third spike is missed entirely. We can achieve faster controlled spiking with a different neuron type. In Fig. 9, we observe the input current and membrane potential of an FS neuron, which we can calculate has an interference-free firing frequency of approximately 33.1 Hz. A spike train at 13 Hz can be generated without a problem, but a 60 Hz spike train misses spikes.

To provide more detailed insight into the generation of spike trains at different frequencies, we measure the RMSE for sequences of 10 spikes (after the distortion-free first spike) as a function of the target firing frequency for RS, FS, LTS, and IB neurons in Fig. 10 where we set the maximum current $I_{\mathrm{max}}=6$ and found the charging time from $I_{\mathrm{max}}$ in Table II. For each type of neuron, the distortion jumps to infinity when we miss a spike. RS neurons are the least accommodating of rapid stimulation, followed by IB neurons, LTS neurons, and then FS neurons. Generally, for each type of neuron, the maximum possible frequency without missing spikes is less than twice that predicted by the interference-free charging and recovery times, e.g., 11 Hz for the RS neuron and 53 Hz for the FS neuron. Thus, the interference-free estimates provide a “rule-of-thumb” to predict achievable firing times in a spike train.

Finally, we measure the distribution of sequence distortions for FS neurons, where we set $I_{\mathrm{max}}=6$ and generate $10^{3}$ realizations of target FS model parameters over the ranges $a\in[0.084,0.1]$, $b\in[0.2,0.21]$, $c\in[-65,-62]$, and $d\in[2,3.2]$. For each realization of target model parameters, we generate actual parameter values that are normally distributed about the target values and with variances that are 1% of the chosen ranges. We use Tables IV and V to determine the target charging and recovery times, and then in Fig. 11 measure the distribution of RMSE distortion as a function of the normalized frequency. The frequencies are normalized to the frequency predicted by the charging and recovery times in Tables IV and V. Even though we are simulating neurons with model parameters that do not match those used to predict the charging and recovery times, the results in Fig. 11 are still consistent with those in Fig. 10, such that FS neurons can be stimulated with RMSE distortion usually below 1 ms if the firing frequency is no more than $50\,\%$ greater than that predicted by the charging and recovery times. This demonstrates the robustness of our methodology to control the firing of individual neurons.

V Conclusion

In this paper, we have considered the use of an optogenetic stimulation model to control the timing of individual neuron spikes. We used the Izhikevich neuron model for the membrane potential dynamics and fitted the neuron charging and recovery times to functions of the model’s parameters and the input current. We have demonstrated that simple functions can help predict lower bounds on the highest firing frequency that can be achieved in regular spiking, fast spiking, low-threshold spiking, and intrinsically bursting neurons with minimal interspike interference. We have also measured deviations due to imperfect knowledge of the neuron model parameters. Future work can consider mismatch between neuron model parameters and experimentally-observed membrane potential dynamics, develop a new model for membrane potential dynamics to align with where light-gated channels are expressed and opened, and study information-theoretic measures for the information that can be embedded in externally-stimulated spike trains.

In Table VI, we list additional equations found for fitting the charging and recovery times of the nominal RS neuron to the maximum input current $I_{\mathrm{max}}$. In Table VII, we list additional equations found for fitting the charging times of the nominal RS neuron to the maximum input current $I_{\mathrm{max}}$ and model parameter $b$.