Root Mean Square Error of Neural Spike Train Sequence Matching with Optogenetics
RRoot Mean Square Error of Neural Spike TrainSequence Matching with Optogenetics
Adam Noel, Dimitrios Makrakis
School of EECSUniversity of Ottawa, Ottawa, Ontario, CanadaEmails: [email protected], [email protected]
Andrew W. Eckford
Dept. of EECSYork University, Toronto, Ontario, CanadaEmail: [email protected]
Abstract —Optogenetics is an emerging field of neurosciencewhere neurons are genetically modified to express light-sensitivereceptors that enable external control over when the neuronsfire. Given the prominence of neuronal signaling within the brainand throughout the body, optogenetics has significant potential toimprove the understanding of the nervous system and to developtreatments for neurological diseases. This paper uses a simpleoptogenetic model to compare the timing distortion between arandomly-generated target spike sequence and an externally-stimulated neuron spike sequence. The distortion is measuredby filtering each sequence and finding the root mean squareerror between the two filter outputs. The expected distortion isderived in closed form when the target sequence generation rateis sufficiently low. Derivations are verified via simulations.
I. I
NTRODUCTION
The nervous system is the most complex system of thehuman body and understanding this system is considered tobe one of the biggest challenges in all of biology; see [1,Ch. 45]. The neural network, with up to connectionswithin the brain, also controls bodily functions such as musclecontraction. The transfer of information is not entirely internal;sensory neurons, such as those in the retina, generate andpropagate signals in response to external stimuli.There is significant interest in developing methods to controlthe external excitation of neurons to improve our under-standing of the nervous system and develop treatments forneurological diseases. One prominent example is the emergingfield of optogenetics; see [2]. Optogenetics uses a relativelysimple genetic modification to induce a neuron to expresslight-sensitive receptors on its membrane. These light-gatedreceptors can then be used to adjust the ion current across themembrane, which enables one to alter its electrical potentialand control when it fires. Experiments in [3], [4] identifiedopsin-based receptors such as Channelrhodopsin (ChR) to beparticularly suitable for optogenetic studies, due to its simplic-ity and its compatibility for implanting in living animals, asfirst demonstrated in the worm C. elegans in [5].From a communication perspective, nanoscale stimulatorswere proposed in [6], [7] to control neurons and interface witha neural network. In [8], it was proposed that such stimulatorscould be implemented using optogenetics and be implanted forlong term use. More generally, the notion of precise neural
This work was supported in part by the Natural Sciences and EngineeringResearch Council of Canada (NSERC). control raises questions about the amount of informationthat can be carried using neurons and the reliability of thatinformation. Information-theoretic analysis of a single ChRreceptor in [9] showed that it has a remarkably high capabilityof receiving information. However, the information propagatedby neurons is typically observed via the pulses that fire and notthe behavior of individual receptors. There is no one universalmethod for neurons to encode information, but researcherstypically measure the number and timing of fired pulses or“spikes”; see [10]. Examples of the importance of timinginclude [11], where neural spike timing patterns in songbirdswere manipulated with millisecond-scale variations to controlrespiratory behavior.In this paper, we measure how effectively we could useoptogenetics to externally stimulate a spike train to match a“target” spike train. We model an ideal neuron that is chargedby a light source. Inspired by the metrics-based approach re-viewed in [10] for natural responses to stimuli, we compare thetarget and generated sequences by measuring the “distance”between them. From our perspective, the distance between thesequences is a distortion between the generated train and thetarget train. In practice, there should be a threshold distortion below which the pertinent information in the spike train can berecovered. We apply simple filter-based metrics and measurethe timing distortion as the root mean square error (RMSE)between the filtered output of the two sequences. We derivethe distortion and approximate the expected distortion whenmatching randomly-generated target sequences.The rest of the paper is organized as follows. Section IIdescribes the neuron firing model. We derive the averagesequence distortion in Section III. We verify our derivationswith simulations in Section IV and conclude in Section V.II. S
YSTEM M ODEL
We consider a system that has a light source and a singleneuron with multiple light-sensitive receptors on its surface.The light source can illuminate the neuron, which opens thereceptor ion channels on its surface and increases its internalpotential until it fires. For the sake of analysis, we will makeseveral (mostly realistic) assumptions about this process:1) The light source is binary, i.e., either on or off. Thisassumption can be appropriate for lasers or LEDs, whichare both common in optogenetics; see [2]. a r X i v : . [ q - b i o . N C ] A ug n on on onoff offNeuron fires and resetson Fig. 1. Illustration of the neuron model with integrate-and-fire. In thisexample, n min = 3 , i.e., spikes must be separated by at least t . In eachinterval when the light source is on, the voltage increases by ∆ V . Once thethreshold τ = 3∆ V is reached, the neuron fires and resets to V ( t ) = 0 .
2) The process of receptors opening is stochastic (e.g.,see [4]), but we assume that the current is equal to itsexpected value I on when the light is on. This assumptionis appropriate if the number of receptors on the neuronis sufficiently large.3) The neuron uses the integrate-and-fire model from [12]with capacitance C and threshold τ . The integrate-and-fire model provides analytical simplicity at the expenseof some fidelity; generalization of this model is possible.4) Time is discretized into slots of ∆ t that are shorter thanthe time necessary to charge the neuron. Specifically,there exists an integer n min such that the integrate-and-fire threshold satisfies τ = n min I on ∆ tC . (1)In continuous time, the minimum firing time is t min = n min ∆ t , but we will focus on the discrete time model.In the following, we give some interpretive statements aboutour assumptions. Let V ( t ) represent the neuron potential asa function of time. From assumptions 1-3, when the light ison, the neuron behaves as an ideal capacitive circuit with acurrent I on . Thus, if the light is on from time t to t , thenthe change in potential V ( t ) − V ( t ) is V ( t ) − V ( t ) = 1 C (cid:90) t t I on dt = I on ( t − t ) C . (2)If the light is off, the current is zero, so V ( t ) − V ( t ) = 0 .From assumption 3, once V ( t ) exceeds the threshold τ , theneuron fires and V ( t ) is immediately reset to zero.From assumption 4, we will say that the light is synchro-nized with the discrete-time clock, and is either on or off foran entire interval ∆ t . Then from (2) we can define ∆ V as ∆ V = V ( t + ∆ t ) − V ( t ) = I on ∆ tC (3)when the light source is on. Finally, since τ = n min ∆ V from(1), the light must be on for n min slots in order for the neuronto fire. This is depicted in Fig. 1.In this work, we consider using the light source to generatea train of spikes to match a target sequence, where we areconstrained by the time it takes to charge and fire the neuron. Target SequenceGenerated Sequence
Fig. 2. Example of target sequence matching in discrete time, where thegenerated sequence is constrained by a charging time of n min = 3 slots. Slotsare labeled chronologically and colored when there is a spike at the start ofthe slot. The target sequence has 4 pulses (shown in blue). The generatedsequence can match the first 2 pulses (green), but the final 2 pulses (yellow)each have a delay of one slot. This matching problem is demonstrated in Fig. 2. We supposethat there is a target spike train (cid:126)u that we define by the timingof its individual spikes, i.e., (cid:126)u = { u , u , . . . , u M } , where u i is the time slot when the i th spike fires. We assume that (cid:126)u isknown a priori . We may not be able to generate (cid:126)u perfectly,e.g., it may be the superposition of spike trains from multipleneurons, but instead we use the light source to generate thesequence (cid:126)v = { v , v , . . . , v N } . The only constraint on neuronfiring times in (cid:126)v that we consider is the time that it takes tocharge up the neuron to the threshold voltage τ . Since (cid:126)u isknown a priori , we can turn on the light source and begincharging the neuron before the corresponding target firingtime. As long as a given target spike is at least n min slotssince the previous spike, then we can generate a correspondingspike at the precise target time. This is a simplified andideal generation model but it facilitates tractable analysis. Aninteresting variation for future work is to refrain from matchinga spike if that enables us to match future spikes, e.g., omitgenerating a spike to match the target in the 7th time slot inFig. 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 (cid:126)v from the sequence (cid:126)u , subject to adistance or distortion measure d ( (cid:126)u, (cid:126)v ) . Obviously, if all spikesin (cid:126)u are separated by at least n min slots, then in our model wecan generate (cid:126)v = (cid:126)u and we should have d ( (cid:126)u, (cid:126)v ) = 0 .III. RMSE IN G ENERATED S PIKE T RAINS
In this section, we present the filter-based metric model from[10] for comparing spike trains and apply it to measure the rootmean square error (RMSE) of the timing distortion betweenthe target sequence (cid:126)u and the generated sequence (cid:126)v . We thensimplify the RMSE when the filter has a length of 1 or 2 timeslots and derive the expected distortion.
A. Filter-Based Model from [10]
There is no one universal method to encode information ina spike train, but in this work we focus on a filter-based metricthat enables some discretion in how to measure the distortion.We begin by mapping the spike trains (cid:126)u and (cid:126)v onto the vectorspace of functions; see [10] for a more general discussion anddditional examples. A discrete time model of the (cid:96) p normdistortion in [10, Eq. (17)] between sequences (cid:126)u and (cid:126)v is d ( (cid:126)u, (cid:126)v ) = (cid:32)(cid:88) n | f [ n ; (cid:126)u ] − f [ n ; (cid:126)v ] | p (cid:33) /p , (4)where n is the time index and f [ n ; · ] is the mapping functionthat maps a sequence to a vector space. We use a filter functionwith a kernel h [ n ] , such that the sequence (cid:126)u maps as f [ n ; (cid:126)u ] = M (cid:88) i =1 h [ n − u i ] . (5)For ease of analysis, we are interested in finite-lengthkernels. Other kernels considered in [10] include the Gaussianfilter and the exponential filter, but they are outside the scopeof this work. Notably, the exponential filter associates themapping function with a neuron’s post-synaptic conductance. B. Filter-Based Metric with RMSE
We now focus on the (cid:96) norm, i.e., the Euclidean distanceor RMSE between the two sequences in vector space, whichwe will see is sensitive to the timing of the individual spikesin (cid:126)u and (cid:126)v . From (4) and (5), the distortion can be written as d ( (cid:126)u, (cid:126)v ) = (cid:88) n (cid:32) M (cid:88) i =1 h [ n − u i ] − N (cid:88) i =1 h [ n − v i ] (cid:33) = (cid:88) n (cid:32) M (cid:88) i =1 h [ n − u i ] (cid:33) + (cid:32) N (cid:88) i =1 h [ n − v i ] (cid:33) − M (cid:88) i =1 N (cid:88) j =1 h [ n − u i ] h [ n − v j ] = M (cid:88) i =1 (cid:88) n h [ n − u i ] (cid:124) (cid:123)(cid:122) (cid:125) Target energy + N (cid:88) i =1 (cid:88) n h [ n − v i ] (cid:124) (cid:123)(cid:122) (cid:125) Generated energy + 2 M (cid:88) i =1 M (cid:88) j = i +1 (cid:88) n h [ n − u i ] h [ n − u j ] (cid:124) (cid:123)(cid:122) (cid:125) Target density ++ 2 N (cid:88) i =1 N (cid:88) j = i +1 (cid:88) n h [ n − v i ] h [ n − v j ] (cid:124) (cid:123)(cid:122) (cid:125) Generated density − M (cid:88) i =1 N (cid:88) j =1 (cid:88) n h [ n − u i ] h [ n − v j ] (cid:124) (cid:123)(cid:122) (cid:125) Overlap measure , (6)where we label the terms in (6). The two energy termsdescribe the energy of the two filtered sequences. The density terms describe the proximity of the individual spikes in eachsequence to the other spikes in the same sequence. The overlap measure describes the proximity of the individual spikes in (cid:126)v tothe spikes in (cid:126)u . It can be shown, as expected, that the distortionis minimized to d ( (cid:126)u, (cid:126)v ) = 0 when the overlap measure ismaximized, i.e., when (cid:126)v = (cid:126)u .We have not yet placed any constraints on the form of thekernel h [ n ] or the two sequences. To simplify the distortionmeasure, we now impose that we generate a sequence of thesame length as the target sequence, i.e., N = M , such that wecan write the timing of each spike in (cid:126)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 d ( (cid:126)u, (cid:126)v ) = M (cid:88) i =1 (cid:88) n h [ n − u i ]+ 2 M (cid:88) i =1 M (cid:88) j = i +1 (cid:32) (cid:88) n h [ n − u i ] h [ n − u j ]+ (cid:88) n h [ n − u i − a i ] h [ n − u j − a j ] (cid:33) − M (cid:88) i =1 M (cid:88) j =1 (cid:88) n h [ n − u i ] h [ n − u j − a j ] . (7)Our notion of “proximity” between spikes when measuringthe overlap or density of sequences is particularly sensitive tothe length of the kernel h [ n ] . To explore this further, we nextconsider kernels of length ∆ W = 1 and ∆ W = 2 discretetime slots. Furthermore, we assume for simplification that thespike times in the target sequence (cid:126)u are all unique (i.e., thereis no more than one spike in a given slot), and without lossof generality that they are sorted in increasing order.
1) Distortion with Kernel of Length 1: If ∆ W = 1 , then thedensity of the target sequence must be 0, i.e., there is no partialoverlap between spikes in the same sequence. Every spike inthe generated sequence either perfectly matches or misses atarget spike; partial overlap between the two sequences is notpossible. Relative to other kernel lengths, the overlap term isminimized. From the perspective of matching (cid:126)u with (cid:126)v , thisdistortion measure discards every generated spike that does notalign with any target spike and the severity of the misalignmentdoes not matter (unlike, for example, a distortion metric thatmeasures the delay). In other words, the timing of the spikesmust be perfectly synchronous to have d ( (cid:126)u, (cid:126)v ) = 0 .Applying a Kronecker delta kernel, where h [ n ] = δ [ n ] , thedistortion in (7) becomes d ( (cid:126)u, (cid:126)v ) = M − M (cid:88) i =1 M (cid:88) j =1 (cid:16) u i ? = u j + a j (cid:17) , (8)where (cid:16) u i ? = u j + a j (cid:17) is an indicator function with value 1if the equality is true. Eq. (8) is in a form that we can readilyevaluate from the target sequence (cid:126)u and the offset sequence (cid:126)a = { a , a , . . . , a M } . We note that (8) does not impose that generated spike has to align with its corresponding targetin order to prevent distortion due to that spike, nor does itmake any assumptions about the individual offset values a j (i.e., they could be any integer). However, since we cannotgenerate multiple spikes simultaneously (because we must re-charge the neuron after every spike generation), we know thatthe indicator function can only be true for at most one valueof j for every value of i (and vice versa).Next, we consider the expected distortion measure d ( (cid:126)u, (cid:126)v ) .In order to do so, we must make additional assumptions aboutthe target sequence and the offset sequence. We impose thateach offset a j must be a delay , such that a j ≥ . We alsoassume that the target sequence is sparse, such that it hasno more than 2 spikes within any interval of n min slots,where n min is the minimum number of slots between spikesin (cid:126)v . Thus, each a j will only depend on the values of thecorresponding u j and u j − , such that we are never waitingto generate more than 1 spike at a time. By imposing theseassumptions, a generated spike that occurs at the same timeas a target spike must be intended for that target, and we canapproximate the distortion in (8) as d ( (cid:126)u, (cid:126)v ) ≈ (cid:32) M − M (cid:88) i =1 (cid:16) a i ? = 0 (cid:17)(cid:33) . (9)Let us consider whether the assumptions for (9) preventus from satisfying (cid:16) u i ? = u j + a j (cid:17) in (8) when a i (cid:54) = 0 . Wecan prove by contradiction that this is true. If a i (cid:54) = 0 , then (cid:16) u i ? = u j + a j (cid:17) could only be true for some i (cid:54) = j . We’veimposed that a i must be non-negative, so we can only consider j < i . The most recent case to satisfy the indicator functionwould be j = i − , such that u i = u i − + a i − . However,if a i − (cid:54) = 0 , then u i − − u i − < n min , and the timing of the ( i − th, ( i − th, and i th spikes violates our assumption thatthere can be no more than 2 target spikes within any intervalof n 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 reducedfor each spike that is generated without delay.We are able to determine the expectation of the approximatedistortion in (9). First, we need the probability that a i = 0 .Since (cid:126)u is increasing, the first offset a = 0 . For i > , weknow that a i = 0 if there is sufficient separation between thecurrent and previous target spikes, i.e., if u i − u i − ≥ n min .To maximize the entropy of the target sequence, we assumethat the number of slots separating consecutive target spikesfollows a geometric distribution with probability p T ; see [13].Thus, we can estimate the probability that a i = 0 as Pr( a i = 0) ≈ Pr( u i − u i − ≥ n min )= (1 − p T ) n min − , (10)and so ( a i ? = 0) , i > , is a Bernoulli random variable withsuccess probability (1 − p T ) n min − . Furthermore, the summa-tion X = (cid:80) Mi =2 ( a i ? = 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 − trials and value x . From the properties of functions of random variables (see[14]), we can then write the expected distortion as follows: d ( (cid:126)u, (cid:126)v ) = E [ d ( (cid:126)u, (cid:126)v )] ≈ M − (cid:88) x =0 (2 M − − x ) p ( x )= M − (cid:88) x =0 (2 M − − x ) (cid:18) M − x (cid:19) × (1 − p T ) x ( n min − (cid:16) − (1 − p T ) n min − (cid:17) M − − x , (11)where p ( x ) is the probability mass function (PMF) of therandom variable X .
2) Distortion with Kernel of Length 2: If ∆ W > , thenpartial overlap between spikes is possible, and the degree ofdistortion is less sensitive to the precise alignment of thespikes. For example, let us consider the case ∆ W = 2 , suchthat h [ n ] = h δ [ n ] + h δ [ n − . The degree of overlapwithin a sequence is still limited and can only occur betweenconsecutive spikes, i.e., the i th spike can overlap the ( i + 1) thspike but not the ( i + 2) th. By applying this constraint, it canbe shown that the distortion in (7) simplifies to d ( (cid:126)u, (cid:126)v ) = (cid:32) M ( h + h ) + 2 M − (cid:88) i =1 h h ( u i +1 ? = u i + 1) − M (cid:88) i =1 M (cid:88) j =1 (cid:104) ( h + h )( u i ? = u j + a j )+ h h ( u i ? = u j + a j ± (cid:105)(cid:33) , (12)which is exact. To simplify (12) and approximate the meandistortion, we assume that every offset must be a delay andthat the target sequence has no more than 2 spikes withinany interval of n min + 1) slots. Then, as in the case where ∆ W = 1 , the only slots that could have overlap between (cid:126)u and (cid:126)v is when a spike can be generated at the same time asits corresponding target with no delay, i.e., a i = 0 . We alsoassume that the number of slots separating consecutive targetspikes follows a geometric distribution with probability p T .Applying these assumptions to (12) leads to d ( (cid:126)u, (cid:126)v ) ≈ (cid:32) M ( h + h ) + 2 h h M − (cid:88) i =1 ( u i +1 ? = u i + 1) − h + h ) (cid:32) M (cid:88) i =2 ( u i − u i − ≥ n min ) (cid:33) (cid:33) . (13)The summation X = (cid:80) M − i =1 ( u i +1 ? = u i +1) is a Binomialrandom variable with M − trials, success probability p T , andvalue x . Analogously to the case where ∆ W = 1 , the sum-mation X = (cid:80) Mi =2 ( u i − u i − ≥ n min ) is a Binomial randomvariable with M − trials, success probability (1 − p T ) n min − ,nd value x . Qualitatively, X is the number of target spikesthat are in the slot immediately following the previous targetspike, and X is the number of target spikes that are at least n min slots after the previous target spike.Next, we approximate the expected distortion. The differ-ence between this case and ∆ W = 1 is that we must determinethe joint PMF p ( x , x ) of two dependent Binomial randomvariables X and X . We derive the joint PMF using themultiplicative rule for joint probabilities, i.e., p ( x , x ) = p ( x | x ) p ( x ) , (14)where p ( x ) is the Binomial PMF with M − trials andsuccess probability p T . Given knowledge of x , there are fewertrials for the occurrence of x (reduced to M − − x ) but ahigher success probability (1 − p T ) n min − . Thus, we can writethe expected distortion as d ( (cid:126)u, (cid:126)v ) ≈ M − (cid:88) x =0 M − − x (cid:88) x =0 (cid:2) (2 M − − x )( h + h )+ 2 h h x (cid:3) p ( x , x )= M − (cid:88) x =0 (cid:18) M − x (cid:19) p x T (1 − p T ) M − − x × M − − x (cid:88) x =0 (cid:2) (2 M − − x )( h + h ) + 2 h h x (cid:3) × (cid:18) M − − x x (cid:19) (1 − p T ) x ( n min − × (cid:0) − (1 − p T ) n min − (cid:1) M − − x − x , (15)which is exact for a low target spike density. However, itsverbosity makes it less intuitive than (11), i.e., when ∆ W = 1 .
3) Consideration of Longer Kernels:
A comparison be-tween (8) and (12) demonstrates the increase in complexitywhen we need to account for partial overlap between fil-tered spikes in the same or different sequences. However,the approximations that we applied to simplify the distortionand derive the expected distortion suggest that a tractableexpression to approximate the distortion for any arbitrary filterlength is feasible. We leave such consideration for future work.IV. N
UMERICAL R ESULTS
In this section, we evaluate the performance of the RMSEtiming distortion metric that we derived in Section III. Weassume throughout this section that the target sequence (cid:126)u has M = 10 spikes that are generated with success probability p T ∈ [10 − , in every discrete time slot. We generate targetsequences by simulating the geometric distribution, and allsimulation results are generated by averaging over at least sequences for each p T , such that every plotted simulationpoint has an insignificant confidence interval. We considera minimum charging time t min = 2 ms , which is consistentwith the typical neuron recovery time; see [1, Ch. 45]. Weuse discrete time slots of length ∆ t = 0 . , so there is aminimum of n min = 4 slots between generated spikes. Target Generation Probability − − − A v e r ag e R M S E w i t h ∆ W = Fig. 3. Average RMSE d ( (cid:126)u, (cid:126)v ) as a function of the target firing probability p T in each time slot. The expected analytical curve is evaluated from (11).The actual observed distortion is calculated using (8), and the approximateobserved distortion is calculated using (9). For the filter of length ∆ W = 1 , we consider the Kro-necker delta kernel so that we can apply the results fromSection III-B1. For the filter of length ∆ W = 2 , we considerfilter coefficients { h , h } = {√ . , √ . } , so that eachcoefficient is weighted equally and the sum of the squaresof the coefficients is equal to that of the Kronecker delta.In Fig. 3, we measure the average RMSE as a function of thetarget spike generation probability p T for the filter of length ∆ W = 1 . We use two methods to calculate the simulateddistortion. The true simulated distortion is measured using(8). We approximate the simulated distortion using (9), wherewe assume that p T is sufficiently low, i.e., that (cid:126)u is sparse.The expected analytical curve is plotted using (11) and alsoassumes that (cid:126)u is sparse.We observe in Fig. 3 that all three curves agree wellwhen p T < , such that the timing of a given generatedspike in (cid:126)v primarily depends on the timing of only oneprevious spike in (cid:126)u . For p T ≥ , we often have two ormore spikes within an interval of n min slots, i.e., spikesoccur sufficiently often that multiple previous spikes in (cid:126)u affect the timing of spikes in (cid:126)v . This generally leads to theexpected distortion acting as a lower bound. However, for very high spike generation probabilities, i.e., p T ≥
40 % , thetrue distortion becomes lower than predicted by the expectedcurve. This is because there are so many target spikes in (cid:126)u that delayed spikes in (cid:126)v are likely to occur at the same timeas future spikes in (cid:126)u , i.e., u i = u j + a j for some j < i .Examples of this are shown in Fig. 4. Such occurrences arenot accounted for in the derivations of the approximations,where spikes in (cid:126)v must align with the corresponding spikes in (cid:126)u , but from (8) these asynchronous overlaps lead to a smallerdistortion. Thus, the expected distortion is a lower bound inthe “low density” regime but an upper bound in the “highdensity” regime. We also note that the approximate simulated Target SequenceGenerated Sequence u u u u u u v v v v Fig. 4. Example of target sequence matching where the target sequence (cid:126)u has a high spike generation probability and the generated sequence (cid:126)v isconstrained by a charging time of n min = 3 slots. (cid:126)u has 6 pulses (at start ofslots shown in blue). (cid:126)v can match the first pulse (green), but the spikes in the th and th slots match future spikes that are not the corresponding targetspikes (yellow). The spike in the th slot does not match any spike (red). Thespikes in (cid:126)v to match u and u occur after the th slot and are not shown. Target Generation Probability − − − A v e r ag e R M S E w i t h ∆ W = Fig. 5. Average RMSE d ( (cid:126)u, (cid:126)v ) as a function of the target firing probability p T in each time slot. The metric kernel has length ∆ W = 2 and coefficients { h , h } = {√ . , √ . } . The expected analytical curve is evaluatedfrom (15). The actual observed distortion is calculated using (12), and theapproximate observed distortion is calculated using (13). distortion converges to the expected distortion as p T → .This is because every spike generated after the initial onehas no overlap with its corresponding target spike, so bothdistortion measures are maximized to the same value, i.e., d ( (cid:126)u, (cid:126)v ) = √ M − √ .In Fig. 5, we measure the average RMSE as a functionof the target spike generation probability p T for the filter oflength ∆ W = 2 . Analogously to the filter of length 1, we usetwo methods to calculate the simulated distortion. The truesimulated distortion is measured using (12). We approximatethe simulated distortion using (13), where we assume that p T is sufficiently low. The expected analytical curve is plottedusing (15) and also assumes that p T is sufficiently low.As in Fig. 3 when ∆ W = 1 , Fig. 5 shows that the expecteddistortion is a lower bound on the approximate simulateddistortion. This bound is accurate for low p T (here when p T < ) and then converges when p T → . However,unlike the ∆ W = 1 case, Fig. 5 demonstrates that the expected distortion is an upper bound on the true distortion for all p T .This is a side effect of the longer filter; non-zero overlapoccurs between a target spike in (cid:126)u and the correspondingspike generated in (cid:126)v when the latter is generated with a one-slot delay. Such “imperfect” overlap reduces the measure ofdistortion calculated using (12) but is not accounted for ineither the approximate or expected distortion.V. C ONCLUSIONS
In this paper, we used a simple optogenetic model toexternally stimulate a neuron and generate a spike train. Weconstrained the neuron’s charging time and measured thedistortion in the spike train as the filtered train’s RMSE froma filtered target sequence. We showed that the expected dis-tortion can be accurately predicted when the spike generationrate in the target sequence is sufficiently low.Ultimately, this is a preliminary work to understand theinformation that can be carried in a sequence of neuron pulses.Future work includes a more complete statistical description ofthe distortion and consideration of other filters and distortionmetrics. We are also interested in comparing our analysiswith experimental neuron firing data, and considering thepropagation of stimulated pulses between connected neurons.R
EFERENCES[1] D. E. Sadava, D. M. Hillis, H. C. Heller, and M. Berenbaum,
Life: TheScience of Biology , 10th ed. Sinauer Associates, 2014.[2] L. Fenno, O. Yizhar, and K. Deisseroth, “The development and applica-tion of optogenetics,”
Annu. Rev. Neurosci. , vol. 34, no. 1, pp. 389–412,Jul. 2011.[3] G. Nagel, D. Ollig, M. Fuhrmann, S. Kateriya, A. M. Musti, E. Bamberg,and P. Hegemann, “Channelrhodopsin-1: A light-gated proton channelin green algae,”
Science , vol. 296, no. 5577, pp. 2395–2398, Jun. 2002.[4] G. Nagel, T. Szellas, W. Huhn, S. Kateriya, N. Adeishvili, P. Berthold,D. Ollig, P. Hegemann, and E. Bamberg, “Channelrhodopsin-2, a directlylight-gated cation-selective membrane channel.”
Proc. Nat. Acad. Sci. ,vol. 100, no. 24, pp. 13 940–5, Nov. 2003.[5] G. Nagel, M. Brauner, J. F. Liewald, N. Adeishvili, E. Bamberg, andA. Gottschalk, “Light activation of channelrhodopsin-2 in excitable cellsof Caenorhabditis elegans triggers rapid behavioral responses,”
Curr.Biol. , vol. 15, no. 24, pp. 2279–2284, Dec. 2005.[6] S. Balasubramaniam, N. T. Boyle, A. Della-Chiesa, F. Walsh,A. Mardinoglu, D. Botvich, and A. Prina-Mello, “Development of arti-ficial neuronal networks for molecular communication,”
Nano Commun.Net. , vol. 2, no. 2-3, pp. 150–160, Jun. 2011.[7] F. Mesiti and I. Balasingham, “Nanomachine-to-neuron communicationinterfaces for neuronal stimulation at nanoscale,”
IEEE J. Sel. AreasCommun. , vol. 31, no. 12, pp. 695–704, Dec. 2013.[8] S. A. Wirdatmadja, S. Balasubramaniam, Y. Koucheryavy, and J. M.Jornet, “Wireless optogenetic neural dust for deep brain stimulation,” in
Proc. IEEE Healthcom , Sep. 2016, pp. 1–6.[9] A. W. Eckford, K. A. Loparo, and P. J. Thomas, “Finite-state channelmodels for signal transduction in neural systems,” in
Proc. IEEEICASSP , Mar. 2016, pp. 6300–6304.[10] C. Houghton and J. D. Victor, “Measuring representational distances:The spike-train metrics approach,” in
Vis. Popul. Codes , N. Kriegeskorteand G. Kreiman, Eds. MIT Press, 2011, ch. 8, pp. 213–244.[11] K. H. Srivastava, C. M. Holmes, M. Vellema, A. R. Pack, C. P. H.Elemans, I. Nemenman, and S. J. Sober, “Motor control by preciselytimed spike patterns.”
Proc. Natl. Acad. Sci. , vol. 114, no. 5, pp. 1171–1176, Jan. 2017.[12] L. F. Abbott, “Lapicque’s introduction of the integrate-and-fire modelneuron (1907).”
Brain Res. Bull. , vol. 50, no. 5-6, pp. 303–4, 1999.[13] S. M. Ross,
Introduction to Probability and Statistics for Engineers andScientists , 4th ed. Academic Press, 2009.[14] J. G. Proakis,