Root Mean Square Error of Neural Spike Train Sequence Matching with Optogenetics

TL;DR

This paper analyzes the accuracy of neural spike train matching in optogenetics by deriving and verifying the root mean square error between target and stimulated spike sequences, providing insights into timing distortion.

Contribution

It introduces a mathematical framework for quantifying spike train matching error in optogenetics, including a closed-form expression for expected distortion at low firing rates.

Findings

Derived a closed-form expression for expected RMS error.

Verified theoretical results with simulations.

Provides a tool for assessing neural stimulation accuracy.

Abstract

Optogenetics is an emerging field of neuroscience where neurons are genetically modified to express light-sensitive receptors that enable external control over when the neurons fire. Given the prominence of neuronal signaling within the brain and throughout the body, optogenetics has significant potential to improve the understanding of the nervous system and to develop treatments for neurological diseases. This paper uses a simple optogenetic model to compare the timing distortion between a randomly-generated target spike sequence and an externally-stimulated neuron spike sequence. The distortion is measured by filtering each sequence and finding the root mean square error between the two filter outputs. The expected distortion is derived in closed form when the target sequence generation rate is sufficiently low. Derivations are verified via simulations.

Figures

Equations

$$\tau =n_{\text{min}}\frac{I_{\mathrm{on}}\Delta t}{C}.$$

$$V(t_2)-V(t_1)=\frac{1}{C}{\int }_{t_1}^{t_2}{I}_{\mathrm{on}}dt=\frac{I_{\mathrm{on}}(t_2-t_1)}{C}.$$

$$\Delta V=V(t+\Delta t)-V(t)=\frac{I_{\mathrm{on}}\Delta t}{C}$$

$$d\left(\vec{u},\vec{v}\right)={\left(\sum _n{\left|f\left[n;\vec{u}\right]-f\left[n;\vec{v}\right]\right|}^p\right)}^{1/p},$$

$$f\left[n;\vec{u}\right]=\sum _{i=1}^{M}h\left[n-u_i\right].$$

$$d\left(\vec{u},\vec{v}\right)=$$

$$=$$

$$-2\sum _{i=1}^M\sum _{j=1}^{N}h\left[n-u_i\right]h\left[n-v_j\right]\text{bBigg@}3]\text{bBigg@}3{)}^{\frac{1}{2}}$$

$$=$$

$$+\underset{\text{Target density}}{\underbrace{2\sum _{i=1}^M\sum _{j=i+1}^M\sum _{n}h\left[n-u_i\right]h\left[n-u_j\right]}}+$$

$$+\underset{\text{Generated density}}{\underbrace{2\sum _{i=1}^N\sum _{j=i+1}^N\sum _{n}h\left[n-v_i\right]h\left[n-v_j\right]}}$$

$$-\underset{\text{Overlap measure}}{\underbrace{2\sum _{i=1}^M\sum _{j=1}^N\sum _{n}h\left[n-u_i\right]h\left[n-v_j\right]}}\text{bBigg@}3{)}^{\frac{1}{2}},$$

$$d\left(\vec{u},\vec{v}\right)=$$

$$\displaystyle+2\sum_{i=1}^{M}\sum_{j=i+1}^{M}\Bigg{(}\sum_{n}h\left[n-u_{i}\right]h\left[n-u_{j}\right]$$

$$\displaystyle+\sum_{n}h\left[n-u_{i}-a_{i}\right]h\left[n-u_{j}-a_{j}\right]\Bigg{)}$$

$$-2\sum _{i=1}^M\sum _{j=1}^M\sum _{n}h\left[n-u_i\right]h\left[n-u_j-a_j\right]\text{bBigg@}3{)}^{\frac{1}{2}}.$$

$$d\left(\vec{u},\vec{v}\right)={\left(2M-2\sum _{i=1}^M\sum _{j=1}^M\left(u_i\stackrel{?}{=}{u}_j+a_j\right)\right)}^{\frac{1}{2}},$$

$$d\left(\vec{u},\vec{v}\right)\approx {\left(2M-2\sum _{i=1}^M\left(a_i\stackrel{?}{=}0\right)\right)}^{\frac{1}{2}}.$$

$$\mathrm{Pr}(a_i=0)\approx$$

$$=$$

$$\overline{d}\left(\vec{u},\vec{v}\right)=$$

$$=$$

$$\times {\left(1-p_{\mathrm{T}}\right)}^{x(n_{\text{min}}-1)}{\left(1-{\left(1-p_{\mathrm{T}}\right)}^{n_{\text{min}}-1}\right)}^{M-1-x},$$

$$d\left(\vec{u},\vec{v}\right)=$$

$$\displaystyle-2\sum_{i=1}^{M}\sum_{j=1}^{M}\Big{[}(h_{0}^{2}+h_{1}^{2})(u_{i}\stackrel{{\scriptstyle?}}{{=}}u_{j}+a_{j})$$

$$\displaystyle+h_{0}h_{1}(u_{i}\stackrel{{\scriptstyle?}}{{=}}u_{j}+a_{j}\pm 1)\Big{]}\Bigg{)}^{\frac{1}{2}},$$

$$d\left(\vec{u},\vec{v}\right)\approx$$

$$\displaystyle-2(h_{0}^{2}+h_{1}^{2})\left(1+\sum_{i=2}^{M}(u_{i}-u_{i-1}\stackrel{{\scriptstyle?}}{{\geq}}{n_{\textrm{min}}})\right)\Bigg{)}^{\frac{1}{2}}.$$

$$p\left(x_1,x_2\right)=p\left(x_2|x_1\right)p\left(x_1\right),$$

$$\overline{d}\left(\vec{u},\vec{v}\right)\approx$$

$$\displaystyle+2h_{0}h_{1}x_{1}\big{]}^{\frac{1}{2}}p\left(x_{1},x_{2}\right)$$

$$=$$

$$\times \sum _{x_2=0}^{M-1-x_1}{\left[(2M-2-2x_2)(h_0^2+h_1^2)+2h_0{h}_1{x}_1\right]}^{\frac{1}{2}}$$

$$\times \left(\genfrac{}{}{0ex}{}{M-1-x_1}{x_2}\right)(1-p_{\mathrm{T}}{)}^{x_2(n_{\text{min}}-2)}$$

$$\times {\left(1-(1-p_{\mathrm{T}}{)}^{n_{\text{min}}-2}\right)}^{M-1-x_1-x_2},$$

Full text

I Introduction

The nervous system is the most complex system of the human body and understanding this system is considered to be one of the biggest challenges in all of biology; see [1, Ch. 45]. The neural network, with up to $10^{14}$ connections within the brain, also controls bodily functions such as muscle contraction. The transfer of information is not entirely internal; sensory neurons, such as those in the retina, generate and propagate signals in response to external stimuli.

There is significant interest in developing methods to control the external excitation of neurons to improve our understanding of the nervous system and develop treatments for neurological diseases. One prominent example is the emerging field of optogenetics; see [2]. Optogenetics uses a relatively simple genetic modification to induce a neuron to express light-sensitive receptors on its membrane. These light-gated receptors can then be used to adjust the ion current across the membrane, which enables one to alter its electrical potential and control when it fires. Experiments in [3, 4] identified opsin-based receptors such as Channelrhodopsin (ChR) to be particularly suitable for optogenetic studies, due to its simplicity and its compatibility for implanting in living animals, as first demonstrated in the worm C. elegans in [5].

From a communication perspective, nanoscale stimulators were proposed in [6, 7] to control neurons and interface with a neural network. In [8], it was proposed that such stimulators could be implemented using optogenetics and be implanted for long term use. More generally, the notion of precise neural control raises questions about the amount of information that can be carried using neurons and the reliability of that information. Information-theoretic analysis of a single ChR receptor in [9] showed that it has a remarkably high capability of receiving information. However, the information propagated by neurons is typically observed via the pulses that fire and not the behavior of individual receptors. There is no one universal method for neurons to encode information, but researchers typically measure the number and timing of fired pulses or “spikes”; see [10]. Examples of the importance of timing include [11], where neural spike timing patterns in songbirds were manipulated with millisecond-scale variations to control respiratory behavior.

In this paper, we measure how effectively we could use optogenetics to externally stimulate a spike train to match a “target” spike train. We model an ideal neuron that is charged by a light source. Inspired by the metrics-based approach reviewed in [10] for natural responses to stimuli, we compare the target and generated sequences by measuring the “distance” between them. From our perspective, the distance between the sequences is a distortion between the generated train and the target train. In practice, there should be a threshold distortion below which the pertinent information in the spike train can be recovered. We apply simple filter-based metrics and measure the timing distortion as the root mean square error (RMSE) between the filtered output of the two sequences. We derive the distortion and approximate the expected distortion when matching randomly-generated target sequences.

The rest of the paper is organized as follows. Section II describes the neuron firing model. We derive the average sequence distortion in Section III. We verify our derivations with simulations in Section IV and conclude in Section V.

II System Model

We consider a system that has a light source and a single neuron with multiple light-sensitive receptors on its surface. The light source can illuminate the neuron, which opens the receptor ion channels on its surface and increases its internal potential until it fires. For the sake of analysis, we will make several (mostly realistic) assumptions about this process:

The light source is binary, i.e., either on or off. This assumption can be appropriate for lasers or LEDs, which are both common in optogenetics; see [2]. 2. 2.

The process of receptors opening is stochastic (e.g., see [4]), but we assume that the current is equal to its expected value $I_{\mathrm{on}}$ when the light is on. This assumption is appropriate if the number of receptors on the neuron is sufficiently large. 3. 3.

The neuron uses the integrate-and-fire model from [12] with capacitance $C$ and threshold $\tau$. The integrate-and-fire model provides analytical simplicity at the expense of some fidelity; generalization of this model is possible. 4. 4.

Time is discretized into slots of $\Delta t$ that are shorter than the time necessary to charge the neuron. Specifically, there exists an integer ${n_{\textrm{min}}}$ such that the integrate-and-fire threshold satisfies

[TABLE]

In continuous time, the minimum firing time is $t_{\textrm{min}}={n_{\textrm{min}}}\Delta t$, but we will focus on the discrete time model.

In the following, we give some interpretive statements about our assumptions. Let $V(t)$ represent the neuron potential as a function of time. From assumptions 1-3, when the light is on, the neuron behaves as an ideal capacitive circuit with a current $I_{\mathrm{on}}$. Thus, if the light is on from time $t_{1}$ to $t_{2}$, then the change in potential $V(t_{2})-V(t_{1})$ is

[TABLE]

If the light is off, the current is zero, so $V(t_{2})-V(t_{1})=0$.

From assumption 3, once $V(t)$ exceeds the threshold $\tau$, the neuron fires and $V(t)$ is immediately reset to zero.

From assumption 4, we will say that the light is synchronized with the discrete-time clock, and is either on or off for an entire interval $\Delta t$. Then from (2) we can define $\Delta V$ as

[TABLE]

when the light source is on. Finally, since $\tau={n_{\textrm{min}}}\Delta V$ from (1), the light must be on for ${n_{\textrm{min}}}$ slots in order for the neuron to fire. This is depicted in Fig. 1.

In this work, we consider using the light source to generate a train of spikes to match a target sequence, where we are constrained by the time it takes to charge and fire the neuron. This matching problem is demonstrated in Fig. 2. We suppose that there is a target spike train $\vec{u}$ that we define by the timing of its individual spikes, i.e., $\vec{u}=\{u_{1},u_{2},\ldots,u_{M}\}$, where $u_{i}$ is the time slot when the $i$th spike fires. We assume that $\vec{u}$ is known a priori. We may not be able to generate $\vec{u}$ perfectly, e.g., it may be the superposition of spike trains from multiple neurons, but instead we use the light source to generate the sequence $\vec{v}=\{v_{1},v_{2},\ldots,v_{N}\}$. The only constraint on neuron firing times in $\vec{v}$ that we consider is the time that it takes to charge up the neuron to the threshold voltage $\tau$. Since $\vec{u}$ is known a priori, we can turn on the light source and begin charging the neuron before the corresponding target firing time. As long as a given target spike is at least ${n_{\textrm{min}}}$ slots since the previous spike, then we can generate a corresponding spike at the precise target time. This is a simplified and ideal generation model but it facilitates tractable analysis. An interesting variation for future work is to refrain from matching a spike if that enables us to match future spikes, e.g., omit generating a spike to match the target in the 7th time slot in Fig. 2 so that we can match the target in the 10th time slot.

Our goal in the remainder of this paper is to measure the “distance” of the sequence $\vec{v}$ from the sequence $\vec{u}$, subject to a distance or distortion measure $d\left(\vec{u},\vec{v}\right)$. Obviously, if all spikes in $\vec{u}$ are separated by at least ${n_{\textrm{min}}}$ slots, then in our model we can generate $\vec{v}=\vec{u}$ and we should have $d(\vec{u},\vec{v})=0$.

III RMSE in Generated Spike Trains

In this section, we present the filter-based metric model from [10] for comparing spike trains and apply it to measure the root mean square error (RMSE) of the timing distortion between the target sequence $\vec{u}$ and the generated sequence $\vec{v}$. We then simplify the RMSE when the filter has a length of 1 or 2 time slots and derive the expected distortion.

III-A Filter-Based Model from [10]

There is no one universal method to encode information in a spike train, but in this work we focus on a filter-based metric that enables some discretion in how to measure the distortion. We begin by mapping the spike trains $\vec{u}$ and $\vec{v}$ onto the vector space of functions; see [10] for a more general discussion and additional examples. A discrete time model of the $\ell^{\mathrm{p}}$ norm distortion in [10, Eq. (17)] between sequences $\vec{u}$ and $\vec{v}$ is

[TABLE]

where $n$ is the time index and $f\left[n;\cdot\right]$ is the mapping function that maps a sequence to a vector space. We use a filter function with a kernel $h\left[n\right]$, such that the sequence $\vec{u}$ maps as

[TABLE]

For ease of analysis, we are interested in finite-length kernels. Other kernels considered in [10] include the Gaussian filter and the exponential filter, but they are outside the scope of this work. Notably, the exponential filter associates the mapping function with a neuron’s post-synaptic conductance.

III-B Filter-Based Metric with RMSE

We now focus on the $\ell^{2}$ norm, i.e., the Euclidean distance or RMSE between the two sequences in vector space, which we will see is sensitive to the timing of the individual spikes in $\vec{u}$ and $\vec{v}$. From (4) and (5), the distortion can be written as

[TABLE]

where we label the terms in (6). The two energy terms describe the energy of the two filtered sequences. The density terms describe the proximity of the individual spikes in each sequence to the other spikes in the same sequence. The overlap measure describes the proximity of the individual spikes in $\vec{v}$ to the spikes in $\vec{u}$. It can be shown, as expected, that the distortion is minimized to $d(\vec{u},\vec{v})=0$ when the overlap measure is maximized, i.e., when $\vec{v}=\vec{u}$.

We have not yet placed any constraints on the form of the kernel $h\left[n\right]$ or the two sequences. To simplify the distortion measure, we now impose that we generate a sequence of the same length as the target sequence, i.e., $N=M$, such that we can write the timing of each spike in $\vec{v}$ as $v_{i}=u_{i}+a_{i}$, where $a_{i}$ is the offset of the $i$th generated spike from the target time. Thus, we can write the distortion as

[TABLE]

Our notion of “proximity” between spikes when measuring the overlap or density of sequences is particularly sensitive to the length of the kernel $h\left[n\right]$. To explore this further, we next consider kernels of length $\Delta W=1$ and $\Delta W=2$ discrete time slots. Furthermore, we assume for simplification that the spike times in the target sequence $\vec{u}$ are all unique (i.e., there is no more than one spike in a given slot), and without loss of generality that they are sorted in increasing order.

III-B1 Distortion with Kernel of Length 1

If $\Delta W=1$, then the density of the target sequence must be 0, i.e., there is no partial overlap between spikes in the same sequence. Every spike in the generated sequence either perfectly matches or misses a target spike; partial overlap between the two sequences is not possible. Relative to other kernel lengths, the overlap term is minimized. From the perspective of matching $\vec{u}$ with $\vec{v}$, this distortion measure discards every generated spike that does not align with any target spike and the severity of the misalignment does not matter (unlike, for example, a distortion metric that measures the delay). In other words, the timing of the spikes must be perfectly synchronous to have $d(\vec{u},\vec{v})=0$.

Applying a Kronecker delta kernel, where $h\left[n\right]=\delta[n]$, the distortion in (7) becomes

[TABLE]

where $\left(u_{i}\stackrel{{\scriptstyle?}}{{=}}u_{j}+a_{j}\right)$ is an indicator function with value 1 if the equality is true. Eq. (8) is in a form that we can readily evaluate from the target sequence $\vec{u}$ and the offset sequence $\vec{a}=\{a_{1},a_{2},\ldots,a_{M}\}$. We note that (8) does not impose that a generated spike has to align with its corresponding target in order to prevent distortion due to that spike, nor does it make any assumptions about the individual offset values $a_{j}$ (i.e., they could be any integer). However, since we cannot generate multiple spikes simultaneously (because we must re-charge the neuron after every spike generation), we know that the indicator function can only be true for at most one value of $j$ for every value of $i$ (and vice versa).

Next, we consider the expected distortion measure $d\left(\vec{u},\vec{v}\right)$. In order to do so, we must make additional assumptions about the target sequence and the offset sequence. We impose that each offset $a_{j}$ must be a delay, such that $a_{j}\geq 0$. We also assume that the target sequence is sparse, such that it has no more than 2 spikes within any interval of $2{n_{\textrm{min}}}$ slots, where ${n_{\textrm{min}}}$ is the minimum number of slots between spikes in $\vec{v}$. Thus, each $a_{j}$ will only depend on the values of the corresponding $u_{j}$ and $u_{j-1}$, such that we are never waiting to generate more than 1 spike at a time. By imposing these assumptions, a generated spike that occurs at the same time as a target spike must be intended for that target, and we can approximate the distortion in (8) as

[TABLE]

Let us consider whether the assumptions for (9) prevent us from satisfying $\left(u_{i}\stackrel{{\scriptstyle?}}{{=}}u_{j}+a_{j}\right)$ in (8) when $a_{i}\neq 0$. We can prove by contradiction that this is true. If $a_{i}\neq 0$, then $\left(u_{i}\stackrel{{\scriptstyle?}}{{=}}u_{j}+a_{j}\right)$ could only be true for some $i\neq j$. We’ve imposed that $a_{i}$ must be non-negative, so we can only consider $j<i$. The most recent case to satisfy the indicator function would be $j=i-1$, such that $u_{i}=u_{i-1}+a_{i-1}$. However, if $a_{i-1}\neq 0$, then $u_{i-1}-u_{i-2}<{n_{\textrm{min}}}$, and the timing of the $(i-2)$th, $(i-1)$th, and $i$th spikes violates our assumption that there can be no more than 2 target spikes within any interval of $2{n_{\textrm{min}}}$ slots. So, we only need to look for cases of $a_{i}=0$, i.e., when spikes are generated with no delay, which leads to (9). From (9), we see that the distortion measure is reduced for each spike that is generated without delay.

We are able to determine the expectation of the approximate distortion in (9). First, we need the probability that $a_{i}=0$. Since $\vec{u}$ is increasing, the first offset $a_{1}=0$. For $i>1$, we know that $a_{i}=0$ if there is sufficient separation between the current and previous target spikes, i.e., if $u_{i}-u_{i-1}\geq{n_{\textrm{min}}}$. To maximize the entropy of the target sequence, we assume that the number of slots separating consecutive target spikes follows a geometric distribution with probability $p_{\mathrm{T}}$; see [13]. Thus, we can estimate the probability that $a_{i}=0$ as

[TABLE]

and so $(a_{i}\stackrel{{\scriptstyle?}}{{=}}0),i>1$, is a Bernoulli random variable with success probability $\left(1-p_{\mathrm{T}}\right)^{{n_{\textrm{min}}}-1}$. Furthermore, the summation $X=\sum_{i=2}^{M}(a_{i}\stackrel{{\scriptstyle?}}{{=}}0)$, which is the number of target spikes (after the initial spike) that we can generate with no delay, is a Binomial random variable with $M-1$ trials and value $x$. From the properties of functions of random variables (see [14]), we can then write the expected distortion as follows:

[TABLE]

where $p\left(x\right)$ is the probability mass function (PMF) of the random variable $X$.

III-B2 Distortion with Kernel of Length 2

If $\Delta W>1$, then partial overlap between spikes is possible, and the degree of distortion is less sensitive to the precise alignment of the spikes. For example, let us consider the case $\Delta W=2$, such that $h\left[n\right]=h_{0}\delta[n]+h_{1}\delta[n-1]$. The degree of overlap within a sequence is still limited and can only occur between consecutive spikes, i.e., the $i$th spike can overlap the $(i+1)$th spike but not the $(i+2)$th. By applying this constraint, it can be shown that the distortion in (7) simplifies to

[TABLE]

which is exact. To simplify (12) and approximate the mean distortion, we assume that every offset must be a delay and that the target sequence has no more than 2 spikes within any interval of $2({n_{\textrm{min}}}+1)$ slots. Then, as in the case where $\Delta W=1$, the only slots that could have overlap between $\vec{u}$ and $\vec{v}$ is when a spike can be generated at the same time as its corresponding target with no delay, i.e., $a_{i}=0$. We also assume that the number of slots separating consecutive target spikes follows a geometric distribution with probability $p_{\mathrm{T}}$. Applying these assumptions to (12) leads to

[TABLE]

The summation $X_{1}=\sum_{i=1}^{M-1}(u_{i+1}\stackrel{{\scriptstyle?}}{{=}}u_{i}+1)$ is a Binomial random variable with $M-1$ trials, success probability $p_{\mathrm{T}}$, and value $x_{1}$. Analogously to the case where $\Delta W=1$, the summation $X_{2}=\sum_{i=2}^{M}(u_{i}-u_{i-1}\stackrel{{\scriptstyle?}}{{\geq}}{n_{\textrm{min}}})$ is a Binomial random variable with $M-1$ trials, success probability $\left(1-p_{\mathrm{T}}\right)^{{n_{\textrm{min}}}-1}$, and value $x_{2}$. Qualitatively, $X_{1}$ is the number of target spikes that are in the slot immediately following the previous target spike, and $X_{2}$ is the number of target spikes that are at least ${n_{\textrm{min}}}$ slots after the previous target spike.

Next, we approximate the expected distortion. The difference between this case and $\Delta W=1$ is that we must determine the joint PMF $p\left(x_{1},x_{2}\right)$ of two dependent Binomial random variables $X_{1}$ and $X_{2}$. We derive the joint PMF using the multiplicative rule for joint probabilities, i.e.,

[TABLE]

where $p\left(x_{1}\right)$ is the Binomial PMF with $M-1$ trials and success probability $p_{\mathrm{T}}$. Given knowledge of $x_{1}$, there are fewer trials for the occurrence of $x_{2}$ (reduced to $M-1-x_{1}$) but a higher success probability $(1-p_{\mathrm{T}})^{{n_{\textrm{min}}}-2}$. Thus, we can write the expected distortion as

[TABLE]

which is exact for a low target spike density. However, its verbosity makes it less intuitive than (11), i.e., when $\Delta W=1$.

III-B3 Consideration of Longer Kernels

A comparison between (8) and (12) demonstrates the increase in complexity when we need to account for partial overlap between filtered spikes in the same or different sequences. However, the approximations that we applied to simplify the distortion and derive the expected distortion suggest that a tractable expression to approximate the distortion for any arbitrary filter length is feasible. We leave such consideration for future work.

IV Numerical Results

In this section, we evaluate the performance of the RMSE timing distortion metric that we derived in Section III. We assume throughout this section that the target sequence $\vec{u}$ has $M=10$ spikes that are generated with success probability $p_{\mathrm{T}}\in[10^{-3},1]$ in every discrete time slot. We generate target sequences by simulating the geometric distribution, and all simulation results are generated by averaging over at least $10^{4}$ sequences for each $p_{\mathrm{T}}$, such that every plotted simulation point has an insignificant confidence interval. We consider a minimum charging time $t_{\textrm{min}}=2\,\mathrm{ms}$, which is consistent with the typical neuron recovery time; see [1, Ch. 45]. We use discrete time slots of length $\Delta t=0.5\,\mathrm{ms}$, so there is a minimum of ${n_{\textrm{min}}}=4$ slots between generated spikes.

For the filter of length $\Delta W=1$, we consider the Kronecker delta kernel so that we can apply the results from Section III-B1. For the filter of length $\Delta W=2$, we consider filter coefficients $\{h_{0},h_{1}\}=\{\sqrt{0.5},\sqrt{0.5}\}$, so that each coefficient is weighted equally and the sum of the squares of the coefficients is equal to that of the Kronecker delta.

In Fig. 3, we measure the average RMSE as a function of the target spike generation probability $p_{\mathrm{T}}$ for the filter of length $\Delta W=1$. We use two methods to calculate the simulated distortion. The true simulated distortion is measured using (8). We approximate the simulated distortion using (9), where we assume that $p_{\mathrm{T}}$ is sufficiently low, i.e., that $\vec{u}$ is sparse. The expected analytical curve is plotted using (11) and also assumes that $\vec{u}$ is sparse.

We observe in Fig. 3 that all three curves agree well when $p_{\mathrm{T}}<4\,\%$, such that the timing of a given generated spike in $\vec{v}$ primarily depends on the timing of only one previous spike in $\vec{u}$. For $p_{\mathrm{T}}\geq 4\,\%$, we often have two or more spikes within an interval of $2{n_{\textrm{min}}}$ slots, i.e., spikes occur sufficiently often that multiple previous spikes in $\vec{u}$ affect the timing of spikes in $\vec{v}$. This generally leads to the expected distortion acting as a lower bound. However, for very high spike generation probabilities, i.e., $p_{\mathrm{T}}\geq 40\,\%$, the true distortion becomes lower than predicted by the expected curve. This is because there are so many target spikes in $\vec{u}$ that delayed spikes in $\vec{v}$ are likely to occur at the same time as future spikes in $\vec{u}$, i.e., $u_{i}=u_{j}+a_{j}$ for some $j<i$. Examples of this are shown in Fig. 4. Such occurrences are not accounted for in the derivations of the approximations, where spikes in $\vec{v}$ must align with the corresponding spikes in $\vec{u}$, but from (8) these asynchronous overlaps lead to a smaller distortion. Thus, the expected distortion is a lower bound in the “low density” regime but an upper bound in the “high density” regime. We also note that the approximate simulated distortion converges to the expected distortion as $p_{\mathrm{T}}\to 1$. This is because every spike generated after the initial one has no overlap with its corresponding target spike, so both distortion measures are maximized to the same value, i.e., $\overline{d}\left(\vec{u},\vec{v}\right)=\sqrt{2M-2}=\sqrt{18}$.

In Fig. 5, we measure the average RMSE as a function of the target spike generation probability $p_{\mathrm{T}}$ for the filter of length $\Delta W=2$. Analogously to the filter of length 1, we use two methods to calculate the simulated distortion. The true simulated distortion is measured using (12). We approximate the simulated distortion using (13), where we assume that $p_{\mathrm{T}}$ is sufficiently low. The expected analytical curve is plotted using (15) and also assumes that $p_{\mathrm{T}}$ is sufficiently low.

As in Fig. 3 when $\Delta W=1$, Fig. 5 shows that the expected distortion is a lower bound on the approximate simulated distortion. This bound is accurate for low $p_{\mathrm{T}}$ (here when $p_{\mathrm{T}}<5\,\%$) and then converges when $p_{\mathrm{T}}\to 1$. However, unlike the $\Delta W=1$ case, Fig. 5 demonstrates that the expected distortion is an upper bound on the true distortion for all $p_{\mathrm{T}}$. This is a side effect of the longer filter; non-zero overlap occurs between a target spike in $\vec{u}$ and the corresponding spike generated in $\vec{v}$ when the latter is generated with a one-slot delay. Such “imperfect” overlap reduces the measure of distortion calculated using (12) but is not accounted for in either the approximate or expected distortion.

V Conclusions

In this paper, we used a simple optogenetic model to externally stimulate a neuron and generate a spike train. We constrained the neuron’s charging time and measured the distortion in the spike train as the filtered train’s RMSE from a filtered target sequence. We showed that the expected distortion can be accurately predicted when the spike generation rate in the target sequence is sufficiently low.

Ultimately, this is a preliminary work to understand the information that can be carried in a sequence of neuron pulses. Future work includes a more complete statistical description of the distortion and consideration of other filters and distortion metrics. We are also interested in comparing our analysis with experimental neuron firing data, and considering the propagation of stimulated pulses between connected neurons.