Homework Module: Evolving Spiking-Neuron Parameters

Purpose: Gain hands-on experience both a) evolving parameters for a single neuron, and b) quantitatively comparing neural spike trains.

1

Assignment

You will use your evolutionary algorithm (EA) to evolve a few critical parameters for a single artificial neuron, which is based on a popular new model by Eugene Izhikevich [1]. These parameters will control the dynamics of the neuron as it runs through many timesteps (i.e. 1000), on each of which you will record the activation level of the neuron to produce an activation time series, often called a spike train for the neuron. The spike train will then be compared to a goal/target spike train (provided to you), and the distance between the two trains will be used to compute a fitness value for the chromosome, all of which is summarized in Figure 1. The two key concepts in this assignment are a) the Izhikevich model, which is a simple, elegant, set of equations that can generate a wide variety of spike trains, depending only upon the setting of a few parameters, and b) metrics for computing the distance between two spike trains. It turns out that judging the similarity between two spike trains, though relatively easy for the human eye, is non-trivial to formalize mathematically in a unanimously-approved manner: there are many different metrics, all of which capture different types of similarities. In this assignment, you will work with a few of the simpler ones.

2

The Basic Izhikevich Spiking-Neuron Model

Izhikevich has received considerable recognition for his spiking-neuron model, whose dynamics can closely match those of a wide variety of neurons via the proper choice of values for a small set of critical parameters. The two key variables of the model are v (the membrane potential or activation level ) and u (a recovery factor), which acts to draw down v. The values of u and v depend upon each other and 5 key parameters: a, b, c, d and k, as shown by the calculations of a) the time derivatives of v and u (denoted v 0 and u0 , respectively), and b) the threshold-driven resets, all of which appear below. The two key variable-update equations are:

v0

=

u0

=

1 (kv 2 + 5v + 140 − u + I) τ a (bv − u) τ 1

(1) (2)

Genome 10.....

a

01.....

00......

b

11.......

c

d

01.......

k

activation

Target Behavior

time

distance Spike Train Distance Metrics

fitness

Figure 1: Illustration of the main elements of the assignment. An EA evolves binary chromosomes that encode the basic parameters for a single neuron (large blue circle). For fitness testing, the neuron is then run for many (e.g. 1000) timesteps to produce a time series of activation values (a.k.a. spike train), which is compared to a desired (a.k.a. target) spike train, with proximity to the target entailing higher fitness. Several different distance metrics will be used to compare spike trains. The upper, orange neuron in the figure is one giving a constant input to the blue neuron, so it has no important parameters other than its output value (which is 10 for this assignment).

2

For this assignment, the effective time constant, τ , is 10, while the external input, I, has a constant value of 10 as well. Remember, these are equations for computing the derivatives, not the actual values, of u and v. The updated values are the old values plus the derivatives. When v exceeds a spike threshold, assumed to be 35 mV for this assignment, the following resets occur:

← c

(3)

u ← u+d

(4)

v

If v exceeds the threshold at time t, then in the activation time series for v, its value at time t should be the threshold or above. But once timestep t is over, you can reset v to c and then use this as the starting value of timestep t+1. The final value for step t+1 will then be the starting value plus the derivative. You are also free to force v to remain at c for the duration of step t+1. The spike patterns in Figure 5 are generated by an algorithm using the former reset strategy. To simulate an Izhikevich neuron, perform the following initializations:

v

←

−60

u ← 0

(5) (6)

Next, run the update (and reset, when necessary) equations for N timesteps (where N is 1000 for this assignment) to generate a time series of N values for v. This is your spike train.

2.1

The 5 Evolved Izhikevich-model parameters

Prior to running the model, each of the 5 parameters: a, b, c, d and k, must be assigned values, which will NOT change throughout the generation of an N-step spike train. For this assignment, the valid ranges for these parameters will be:

• a ∈ [0.001, 0.2] • b ∈ [0.01, 0.3] • c ∈ [−80, −30] • d ∈ [0.1, 10] • k ∈ [0.01, 1.0]

Your evolutionary algorithm’s chromosome will encode values (within these ranges) for these 5 parameters. They will strongly determine the behavior of the u-v update and reset calculations, and thus the shape of the spike train.

3

3

Spike Train Distance Metrics (SDMs)

Deciphering the neural code is clearly one of nature’s most daunting challenges. Despite extensive knowledge of the underlying physicochemical processes that govern neural signal transmission, the functional significance of any particular neural firing pattern is rarely transparent. In fact, merely detecting a pattern becomes an ominous task in the face of multiple neurons with complex temporal spiking patterns. In the simple case, a set of neurons that consistently fire together is a strong candidate for a functionally significant cluster, and detection of that neural assembly from spike-train data is straightforward. However, subtle differences in spike trains can obscure cluster recognition. The scientific task is therefore to formalize relevant spike-train similarity metrics. In experimental neuroscience, these metrics assist in single-neuron cases by determining whether a neuron’s responses to two different stimuli are indeed functionality distinct or merely superficially-distinct isomorphs with the same downstream effect. Also, in multi-neuron experiments, a single tetrode may detect bits and pieces of different spike trains during a time window without revealing the actual number of neurons being probed. Similarity metrics can help group spike trains into equivalence classes to give an indicator of this cardinality. Of course, metrics are also essential for the behavior-based comparisons and classifications of neurons within and between brain regions. For example, [2] uses a novel metric to assess whether hippocampal CA1 cells within a particular standard class, such as spiking pyramidals, have similar or distinct spiking patterns. A wide variety of similarity metrics (a.k.a. distance metrics) are used by experimental and computational neuroscientists for comparing time series of neuronal action potentials. One factor that commonly differentiates these metrics is the degree of focus on the depolarization spikes. Whereas some metrics make direct comparisons of actual spike times, others focus on the length of the time interval between spikes; still others consider aggregates of spike times. At the other extreme, some metrics give no special status to spikes and simply compare adjacent points across two activation-level time series. In this project, we will consider three of the more basic SDMs, those based on a) exact spike times, b) intervals between spikes, and c) general waveforms (with no focus on the spike times), each of which is very straightforward to compute, given the spike trains of two neurons. The basic calculations for each of these is explained below. For the spike-time and spike-interval metrics, the spike train must be preprocessed to find the time points of actual spikes. This can normally be done by the following process:

1. Define an activation threshold, T, above which the neuron is considered to be spiking. For the Izhikevich model described above (and others that closely mimic real neuron behavior) , a typical value of T is 0 mV (millivolts), since most neurons spend a majority of their time with negative activation levels. 2. Move a k-timestep window (e.g. k = 5) along the spike train, and any activation value that is a) in the exact middle of the time window, b) above T, and c) the maximum of all activations in the time window, is considered a spike. 3. By using this window, you avoid the double-counting of spikes, which will often consist of a few time points above the threshold.

Using this (or a similar) algorithm, you will compute a list of spike times to accompany your spike train, where the former simply lists the time points at which spikes occur, while the latter lists the activation levels for every time point of the simulation. The three basic SDMs are now described in detail. 4

3.1

Spike Time Distance Metric

This metric simply compares the times at which corresponding spikes occur, giving reduced similarity (increase distance) with increasing gaps. Two spikes correspond when they have identical indices in their respective spike-time lists. The spike-time distance between two spike trains, T a and T b is defined as:

1 dst (T , T ) = N a

b

"N −1 X

# p1 |tai

−

tbi |p

(7)

i=0

where tai and tbi are the times of the ith spikes in trains a and b, respectively, and N is the minimum of the total spike counts in the two trains. The parameter p typically has a value of 2. Figure 2 shows the basis for this metric. In essence, this metric finds the average length of the thin green rectangles in the figure, with larger values indicating reduced spike-train similarity, i.e. greater spike-train distance.

Figure 2: Illustration of the basic spike-time distance metric, which sums the differences between corresponding spike times (as shown by the thin rectangles on the time axis). Note that the lower spike train is inverted for ease of presentation.

3.2

Spike Interval Distance Metric

Spike intervals are simply the gaps between successive spike times. This metric compares the lengths of corresponding intervals, giving greater similarity for smaller differences. The spike-interval distance between two trains is defined as:

1 dsi (T a , T b ) = N −1

"N −1 X

# p1 |(tai − tai−1 ) − (tbi − tbi−1 )|p

i=1

5

(8)

where T a , T b , tai , tbi , N and p are as defined for equation 7. Figure 3 shows the basis for this metric. In the figure, note that each thin green rectangle (which denotes a spike interval) is bounded by two thicker gray rectangles. It is these gray rectangles that denote the differences between two spike intervals; the metric essentially computes the average combined length of pairs of these gray rectangles. As with the spike-time distance metric, this one can be computed directly from the spike-time list, without requiring values from the original spike train.

Figure 3: Illustration of the basic spike-interval distance metric, which sums the differences between the lengths of the time intervals between corresponding spike times. (as shown by the thick rectangles on the ends of the thinner rectangles). Note that the bottom spike train is inverted for ease of illustration.

3.3

Waveform Distance Metric

The waveform metric is the simplest of all. It compares cotemporaneous activation levels across the two spike trains, with no special treatment for spikes. Hence, it does not use the spike-time list at all. The waveform distance metric is defined as:

1 dwf (T , T ) = M a

b

"M −1 X

# p1 |(via

−

vib |p

(9)

i=0

where M is the number of sampled time points and via and vib are the voltages (i.e. activation levels) at time point i along spike trains a and b, respectively. Figure 4 portrays this comparison of voltage curves. Each thick, green, vertical line represents the difference between corresponding activation levels of the two spike trains; and the waveform metric essentially computes the average height of these lines.

6

Figure 4: Illustration of the basic waveform distance metric, which sums the differences between corresponding points on two activation time series (a.k.a. spike trains). The two time series are displayed as solid blue and dotted red curves, respectively, with corresponding-point distances shown as solid, green, vertical lines. Only some distances are shown.

3.4

A Spike-Count Difference Penalty

Most formal SDMs are designed to compare realistic spike trains, often those based on real brain data. Hence, although the spike trains being compared may differ significantly, they will both have numerous spikes. Conversely, in this assignment, you are generating spike trains from a parameterized model, and since evolution is generating those parameters, it may produce very odd spike trains, some with extremely many spikes, others with none at all. The formal SDMs are not necessarily designed to handle these odd spike trains. If one train has 100 spikes and the other has zero, some of the SDMs may not register much of a difference - since they find no corresponding spikes in the two trains from which to begin accumulating differences. To account for this factor, you will probably want to include a basic distance (i.e. penalty) based on the spike-count difference. With most of the metrics, if spike train 1 has M spikes, while train 2 has N (with N > M ), then the remaining N-M spikes will not be accounted for. One simple way to incorporate these extra spikes into the distance is to make the following assumptions:

• The extra spikes are evenly distributed about the entire spike train. L L , where M is the • The average distance from one of these spikes to the nearest spike of train 1 is 2M average temporal difference between two spikes in train 1 (i.e. the period of the wave that is punctuated by a spike). L is the length of the entire spike train.

Thus, the penalty is simply

L 2M

per extra spike, or

(N −M )L . 2M

In looking at equations 7 and 8, notice that the total distance is being averaged across all the corresponding spikes. Consequently, the penalty should probably be added to the summation prior to averaging; otherwise it could totally dominate the result. Other, more complex penalties, can surely be devised, but this one provides one simple way to supplement any SDM that relies on spike times to line up and compare spike trains. No extra penalty is necessary for the waveform SDM, since it gives no special status to spike times in its calculations. 7

Spike-Train

40

Spike-Train

40

20

Activation-Level(mV)

Activation-Level(mV)

20 0

Act

20

40

60

800

200

400

Time(ms)

600

800

200

400

600

Time(ms)

800

1000

Spike-Train

20

Activation-Level(mV)

Activation-Level(mV)

40

40

20 0 20

Act

40 60

0

Act

20 40 60

80 1000

Act 20

600

1000

Spike-Train

40

0

200

400

600

Time(ms)

800

800

1000

200

400

600

Time(ms)

800

1000

Figure 5: The 4 target spike trains used for this assignment, with correspondence to the following spike-train files: izzy-train1(top left), izzy-train2 (top right), izzy-train3 (bottom left), izzy-train4 (bottom right).

4

Target Spike Trains

The target spike trains for this assignment are provided in data files downloadable from the same web page as this assignment. Each file consists of a single (long) line of 1001 activation levels (for time points 0 to 1000). Plots of each of these spike trains appear in Figure 5.

5

Test Cases

The assignment involves 12 test cases: the cross product of the 4 target spike trains and the 3 SDMs. So the first test case involves spike-train 1 as the target and the spike-time metric as the SDM, case 2 is spike-train 1 using the spike-interval metric, etc.

6

Deliverables 1. An overview description of your system in text and diagrams. Include a description of your genotype representation and fitness function. 5 points 2. For each of the 12 test cases, show (graphically) the closest-to-target spike train that your EA finds, along with its fitness value. Use diagrams similar to those in Figure 5. Also provide a plot of the

8

fitness progression for that run, and list the 5 parameter values (a,b,c,d and k) encoded by the bestof-run genotype. For each case, list the EA parameter values that were most effective in finding good solutions. These include the mutation and crossover rates, population size and parent selection mechanism. These may be different for different test cases, so indicate them for each case and discuss any cases where unusual combinations of these parameters seemed to work best. Include a brief description (in Norwegian or English) of your results for each of the 12 cases. 6 points 3. Based on what you’ve learned about the different representational forms in this course, how would you classify the genotype-phenotype mapping in this exercise? Explain. 1 point 4. Discuss the practical implications of the tool that you have built. How might a computational neuroscientist use it? 1 point 5. Can you imagine other problem domains where a more general version of this same tool could be used? Briefly discuss the prerequisites for usage of the tool in one such domain, and sketch the overall system (in a manner similar to Figure 1 above). 2 points

The report should be between 5 and 10 pages long, including diagrams.

6.1

Warnings

There are a few different published versions of Izhikevich’s model. You are free to use any of them, but be aware that the target spike trains in Figure 5 were all generated using the version described in this document. Depending upon the actual course and semester - your instructor uses this module in various courses - you may or may not be required to do any or all of the following:

1. Demonstrate this module to the instructor or a teaching assistant. 2. Upload the report and/or code for this module to a particular site such as It’s Learning.

Consult your course web pages for the requirements that apply.

References [1] E. Izhikevich, Simple model of spiking neurons, IEEE Transactions on Neural Networks, 14 (2003), pp. 1569–1572. [2] L. G. Jose Ambros-Ingerson and W. R. Holmes, A classification method to distinguish cell-specific responses elicited by current pulses in hippocampal ca1 pyramidal cells, Neural Computation, 20 (2008), pp. 1512–1536.

9