NDVis

Separating fast and slow spikers in exp3

2008-08-13T08:05:00.000-07:00

In the last post we talked about various methods for optimizing functions that are very expensive to compute. In this post I'll talk about applying some of those techniques to separating a few neuron model regions.

In this post we use a very simple optimizer for a multivariate function f(x1,...,xn). The optimizer starts with an intial quess g=(g1,...,gn) and an initial steps size d and then for each dimension i of the domain it attempts to improve the guess g by adding or subtracting d from gi and recalculating f on these new points. If either gives a lower bound, it replaces g with the new guess giving a better minimum and repeats the process with 2*d as the step size. Eventually no improvement will be obtained for that dimension and the process continues to the next dimension. If none of these 2n calculations have provided a better guess, it divides d by two and repeats the process. Otherwise, it just repeats the process with the new guess and the same d. This process is iterated some fixed number of times.

The Fast/Slow Spiker boundary
We tried this with the fast slow spiker boundary starting with the initial guess of (0,0,0,0,0,0,0,0,0) and step size 1 running for 100 iterations through all 9 parameters of the error function. The error function we chose is the sum of the relative predition errors. That is the sum of the percentage of Fast spikers that are misclassified by the hyperplane and the the percentage of Slow spikers that are misclassified. We defined fast spikers to be those with periods less than 0.1s and the slow spikers to be those with periods in the range [0.1,0.42].

Here is a link to the scheme code used to run this optimization.

After 100 iterations, the simple optimizer found a hyperplane whose error rate is 3.24% corresponding to 4594 out of 211106 misclassified fast spikers (2.18%) and 87 out of 8261 misclassified slow spikers (1.05%). The resulting hyperplane was


  0.0                   * na
  0.9914505598049774    * cat
  0.09669412841985557   * cas
 -0.08057844034987964   * a
 -0.03223137613995186   * kca
 -0.01074379204665062   * kd
  0.00537189602332531   * h
  0.0)                  * leak
<
  0.25516506110795223

and the hyperplane coefficients (H) on the left side of the inequality have been normalized to have norm 1.

Another approach to finding a separating hyperplane is to look at all of the boundary points (i.e. points with Hamming distance one which belong to different domains) and then the find a best fit hyperplane for those boundary points using SVD. When we used that approach with the fast slow spikers with the original PBM data we found a hyperplane with an error rate of 4.539%.

When we ran our optimization routine using the SVD generated hyperplane we get a much better result with an error of 1.4% (1163 fast and 70 slow spikers misclassified)
compared to 3.2% (4594 fast and 87 slow). That hyperplane is as follows:


 -0.017757873435532652*na  +
  0.9794311939718561  *cat +
  0.1431868622579713  *cas +
 -0.10587125724245507 *a   +
 -0.06542288691577072 *kca +
 -0.06622104704782268 *kd  +
  0.0038126574778711846*h  +
 -0.0028418323284959096*leak
<
  0.010050138208689753

Our next step will be to try one of the smarter optimization routines and see if this outperforms the result obtained by using SVD to get a good initial guess.

The S5-nonS5 boundary
We then ran the same experiment on the S5 boundary (defined as those silent neurons whose membrane potential is less than -42mV). The resulting hyperplane is


 0.16947428370562861 * na +
 0.5028910495715982  * cat +
 0.7643315339676329  * cas +
-0.29368837294981337 * a   +
-0.01206938518971836 * kca +
-0.05230066915544622 * kd  +
 0.03645960109394088 * h   +
-0.20907695385939198    * leak
<
0.16193091796205464

This is satisfied by


  144241 = 97.95% of the 147262 S5 silent neuron models


   3021 = 2.05% of the S5 silent neuron models are misclassified

and it is also satisfied by only


  36714 = 2.4% of the 1532354 other neuron models

The total error rate is 2.4%+2.05% = 4.45%.

At this point we may want to use NDVis to visualize the quality of this hyperplane classification. Perhaps most of the incorrectly classified neurons are in one section of the space and we can use a piecewise linear function to better classify that the space as the union of two linear domains....

We will also want to apply this approach to several of the other boundaries so as to get quantitative rather than just qualitative results.

Multidimensional Optimization to find separating hyperplanes

2008-08-12T11:32:00.000-07:00

The problem we are considering is how to find simple numerical constraints that can be used to determine the type of a neuron in terms of the maximal conductance values. For example, if we let S1 be the set of fast spiking neuron models (with period < 0.1s)
and S2 the slow spikers (with period > 0.1s and < 0.38s). Then we want to find a hyperplane H and an offset C such that for any 8-tuple X of maximal conductance values with each Xi in the range [0,5] then the two regions are well classified by the test


  F_H(X) == (H*X >= C)

that is, if we define


  R1(H) = {x: F(x) = true}
  R2(H) = {x: F(x) = false}

Then most of S1 will be in R1 and most of S2 will be in R2. We can formally measure
the fitness of the hyperplane H by summing the relative errors of this separation


  E(H) = |(R1(H) * S2|/|S2| + |R2(H)*S1|/|S1|

where A*B denotes the intersection of sets A and B and |A| is the size of the set A.

Here is a nice article on optimization techniques for this type of problem ...
Optimization by direct search: New perspectives on some classical and modern methods (2003) SIAM Review
by Tamara G. Kolda, Robert Michael Lewis, Virginia Torczon
We will try a few of these methods for the fast/slow spiker case.

The Kolda,Lewis,Torczon paper mentions that the simplest general optimization method is the one used in the 50's by Enrico Fermi and Nicholas Metropolis in which they do a search for the local minimum starting at point x and proceeding to nearby grid points x+D*ei or x-D*ei where e1,e2,...,en is the basis of the underlying vector space and D is a constant, initially 1. If they reach a point which is a local minimum in that grid, then they replace D with D/2 and continue until D is sufficiently small. This is a very simple method but will often converge slowly or not at all in practice and is a kind of simulated annealing approach.

A better method, by Powell (1964) is to start with the same initial basis e1,...,en but to then try each basis vector in turn and move in that direction so as to minimize the objective function. The method proceeds in steps starting with a point x0 and generating new points x1, x2, .... It also changes the basis at each step to better reflect the geometry of the optimization function. The jth step proceeds from point pj as follows. Find new vectors


 pi = xj + a1*e1+a2*e2+...+ai*ei

where the ai is chosen by trying several points around 0 (perhaps on a grid or using binary search etc.) Thus, we find an optimization point in the e1 direction, and then move in the e2 direction to a better point, etc. Finally, we search in direction pn to find the local minimum


   x(j+1) = xj + b*(pj-xj)

for some appropriate b. They then replace ej with pj and proceed to j+1 (moding by n as appropriate).

The Powell method is often effective because it modifies the basis as it goes along to reflect the properties of the optimization function. It is also very fast on quadric optimization problems.

Another interesting class of optimizer works by moving a simplex through the space. The idea is to generate n+1 points x0,x1,...,xn and then to select the point xi at which the function is worst (e.g. takes the maximum value). The algorithm then reflects that point through the hyperplane consisting of the other points and evaluates there. If the value is lower it then continues. Otherwise, it takes the vertex with the next worst value and continues. If all vertices have been tried and none of the reflections provide a better value, then the algorithm halts or one can try shortening some of the edges. For example, for each vertex xi, try searching along the line from the projection of xi into the plane formed by the other vertices and look for a better value.

The next step is to try some of these optimization techniques to separate different regions. We'll report on that in a future post.

Exploring the use of QHull and NDVis

2008-08-10T07:53:00.000-07:00

In this post I will document an attempt to study the boundary of a region in the PBM database using the convex hull approach. The region we consider is the set of slow spiking neuron models (i.e. whose period is between 0.1 and 0.38 seconds) which have low values of KCa and CaT (both in the range [0,0.5]).

Our first step is to use the database exp3 described above and pull out those neuron models using the following query:


select na,cat,cas,a,kca,kd,h,leak from exp3 
    where cat=0 and kca=0 and ncat=1 
      and val >0.1 and val < 0.38;

This query generates 2489 rows.

Next we add two lines to the top of this file giving the dimension (8) and the number of points (2489) and then we run the qhull program


[macbookpro:NDVIS/experiments/exp3] tim% time ~/Desktop/qhull-2003.1/src/qhull s n A0.97 FA Fa < slowsp00.txt > slowsp00.out

Convex hull of 2489 points in 8-d:

  Number of vertices: 1222
  Number of facets: 2502
  Number of non-simplicial facets: 2387

Statistics for:  | /Users/tim/Desktop/qhull-2003.1/src/.libs/qhull s n A0.97 FA Fa

  Number of points processed: 1449
  Number of hyperplanes created: 12910511
  Number of distance tests for qhull: 340125400
  Number of distance tests for merging: 247209498
  Number of distance tests for checking: 2112583
  Number of merged facets: 2928370
  CPU seconds to compute hull (after input): 1930
  Approximate facet area:   4.7217202
  Approximate volume:       0.98670223
  Maximum distance of vertex below facet: -1.7 (0.4x)

1884.274u 47.013s 1:24:14.80 38.2% 0+0k 0+10io 2pf+0w
[macbookpro:NDVIS/experiments/exp3] tim%

The flags we used for this call are
s - to create a summary
n - to show the normals and offsets to the facets
A0.97 - merge any neighboring facets where the angle alpha between them is relatively small, that is where cos(alpha) > 0.97
FA - show the total area of the facets in the summary
Fa - print a table showing the area of each individual facet

This generates an output file with 2502 rows. The output consists of the dimension (9)
and the number of rows (2502), each on its own line, followed by the rows. Each row consists of 9 numbers. The first 8 are the coordinates of the normal vector n for the face and the ninth is the constant c. The vector n is normalized to have length 1. The hyperplane for that face consists of the set of all vectors x such that


 x.n=c

The next step is to load this into a database and visualize it using NDVis.
We could also use a k-mean clustering program or some other clustering program to try to generate a smaller set of normals to the faces of the region.

To do this we store the 2502 rows in a file called slowsp00facets.txt and then load that file into mysql as follows:


mysql> use neuro;
mysql> create table expslowbdr(
    na double, cat double, cas double, a double, 
    kca double, kd double, h double, leak double, 
    offset double);
mysql> load data infile 
  "/Users/tim/Desktop/CURRENT/NDVIS/experiments/exp3/slowsp00facets.txt" 
  into table expslowbdr fields terminated by ' ';

Now we fire up NDVis and create a color map to visualize the density of the normals inside the space [-1,1]^8. For this example we divide the intervals into 5 regions
centered about zero using the following function to convert doubles x in the range [a,b] to integers in the range [0,4]
floor((x-a)/(b-a)*5)
For the initial range [-1,1] this results in five intervals of size 0.4:


[-1,-0.6] [-0.6,-0.2] [-0.2,0.2] [0.2,0.6] [0.6,1.0]

We'll start with [-1,1] and then we may zoom in if needed.

Here's the query we use


select
     greatest(0,least(4,floor((na  +1)/2*5))) as xna,
     greatest(0,least(4,floor((cat +1)/2*5))) as xcat,
     greatest(0,least(4,floor((cas +1)/2*5))) as xcas,
     greatest(0,least(4,floor((a   +1)/2*5))) as xa,
     greatest(0,least(4,floor((kca +1)/2*5))) as xkca,
     greatest(0,least(4,floor((kd  +1)/2*5))) as xkd,
     greatest(0,least(4,floor((h   +1)/2*5))) as xh,
     greatest(0,least(4,floor((leak+1)/2*5))) as xleak,
     log(log(1+count(*))) as r from expslowbdr group by xna,xcat,xcas,xa,xkca,xkd,xh,xleak

Here is the image we get of this query with the projection x=(a,h,kca,cat)
and y=(kd,cas,na,leak). You may need to click on the image to see a larger version.

Note that we have created the gray background by using a NDVis projection with 6 values for each dimension, and then using a query that specifies that any point with a 5 in any parameter should be gray. This makes it easy to read the NDVis coordinates of any non-gray point ...

We see from this plot that the QHull is not that effective, there are some patterns with the normals (e.g. no normals with kd < -0.2, and very few with a<-0.2 or cas> 0.2, but the set of normals is quite large. This could be because the surface is broken up into many small facets and is very bumpy (due to the use of the random perturbation of the grid). Also, if there are many facets where the surface of the domain turns a corner, this could result in normals in all directions between the two major faces.

A better approach might be to use the area of each facet to provide a weighting scheme and the to only show those pixels that correspond to normals with a high enough weight. The lesson for this experiment is that using QHull and NDVis directly does not provide enough information to get an informative quantitative description of the domain.

On Quantifying NDVis patterns

2008-08-08T12:19:00.001-07:00

In several of the previous posts we have seen the NDVis allows one to visually detect interesting relationships involving several simultaneous variables. In this post I want to talk about a few ways in which those visually observed patterns translated into numerical constraints on the parameters.

There are two basic approaches. The first isto find a formula for the boundary between two neighboring regions. The second is to find formulas for the boundaries between a region and all other neighboring regions. In this note I will just talk about a few different algorithms that could be used in these situations.

Method 1: calculating and curve fitting the boundary explicitly
One of the simplest methods for finding a hyperplane that separates two neighboring regions is to approximate the boundary and the to find a hyperplane that best fits those points using Singular Value Decomposition (SVD) for example. We have tried this for a few different boundaries. For example, the silent neurons form two or three regions visually and the largest region (which we called S5) displays a strong linear pattern.

One way to approximate the boundary between two neighboring regions in high dimensional space is to look for pairs of points P,Q that have Hamming distance 1 (i.e. which differ by 1 in exactly one dimension), and which lie in different regions. For each such neighboring pairs, the average (P+Q)/2 is a point that lies approximately on the boundary (indeed the error should be no more than 0.5). In the case where one uses a randomly perturbed grid one can look at pairs whose "parent" grid points have Hamming distance 1 and again use their midpoint as an approximation of the boundary.

We can automate this process by iterating through all grid points and all n dimensions and then testing whether or not the 2*n neighbors of the point are in the same region. If they are not, then we add the midpoint to the set of boundary points.

The calculation of the approximating hyperplane can be done using the SVD procedure in R or octave.

Method 2: Simulated Annealing
Another approach is to start with an initial hyperplane and then use simulated annealing to modify its parameters in an attempt to get a better fit. One error measure that might be appropriate would be the weighted sum of the errors, e.g. if R1 and R2 are the regions and H is the affine hyperplane candidate, let s1 be the number of points if R1 that are on the wrong side of the hyperplane and s2 the number of points of R2 that are similarly misclassified, then for an error measure we can use the weighted average of s1 and s2,
E(H) = s1(H)/|R1| + s2(H)/|R2|
There are many kinds of machine learning algorithms that could be applied to this problem. We haven't yet tried any of these approaches.

Method 3: Computing Convex Hulls
Another approach is to compute an approximate convex hull for each region and to use the (n-1)-dimensional facets of that convex hull as the hyperplanes that separate it from neighboring regions. The are several algorithms that one could use, the qhull is probably the best known http://www.qhull.org/

Skiena's book, The Algorithm Design Manual has a nice section on convex hull algorithms and also mentions the "Wrapping Paper" algorithm for higher dimensional polytopes.

What we really want is an algorithm for partitioning a subset of a hypercube into convex domains that closely approximate some feature-based decomposition. Moreover, it would be great to allow non-linear boundaries (e.g. quadrics rather than hyperplanes). We could try to handle the non-linear case by mapping the points into a higher dimensional space (e.g. x=(x1,...,xn) --> (x1^2, x1*x2,...,x1*xn, x2^2,...., xn^2) )
and then applying the convex hull algorithm.

I've tried the qhull algorithm on the slow spikers, but the problem is that it generates a huge polytome with many millions of sides when ideally I want something with only a few sides. The goal is to find a few hyperplanes that well approximate the domain and are simple enough to capture fundamental properties of the domain.

Viewing Single Spike neuron models by period

2008-08-01T10:38:00.001-07:00

In this post we look at those neuron models whose steady state voltage plot has a single spike. In PBM they observed that plotting spiking neurons by frequency yielded a bimodal distribution with two essentially separate distributions. Also, in PBM they separate tonic spikers from 1 spike bursters using an area calculation, in exp3 we did not calculate that area feature so we cannot separate tonic spikere from one spike bursters. For this post we will call these single spike neurons they can be tonic spikers or bursters.

Given that understanding, we can plot the periods of all single spike neurons as shown below and we see that there are several peaks

and we can use these peaks to define several domains of single spike neuron models, classified by their periods:

D1: period in 0.05-0.10 sec with peak at 0.016 sec
D2: period in 0.10-0.38 sec with peak at 0.30 sec
D3: period in 0.38-0.78 sec with peak at 0.50 sec
D4: period in 0.78-3.00 sec with peak at 0.92 sec

The following colormap colors neurons in these domains using red, green, yellow, and purple respectively:

and here is an ndvis view of this colormap with projection kca-na-cas-h for x and cat-kd-a-leak for y

In this plot the four domains D1,D2,D3,D4 all seem to map to spatial separate regions which we can call R1,R2,R3,R4 in this post. You may need to click on the image to get a larger version where it is easier to see these four regions.

Observe that D1 is concentrated in the region roughly given by
KCa<2 Kd>=1, Na>=1, and as KCa increases it tends to thin out where Na+Kd is low
Moreover there is a tendency for D1 to be in the lower right corners of the Na*Kd cells.

The D2 region appears only when CaT=0 and where Kd>Na (i.e the upper left corner of the level 2 cells, that is the Na*Kd cells)

The D3 and D4 regions satisfy
KCa>=1, (Na<1 or Kd<1 or Na+Kd<2)
to get a better way of separating these regions we need to look at another projection.

Below is the image for the projection with x= na-cas-kca-h and y=kd-cat-a-leak

Observe that the D4 region mainly satisfies CaS<2, while the D3 region has CaS>2.

Although we haven't tried to verify it yet, the correspondence between these domains and the PBM domains seems to be roughly the following:

D1: fast spikers of PBM
D2: slow spikers of PBM
D3 and D4: one spike bursters of PBM

In a later post we will investigate methods for finding hyperplane equations for the boundaries of these regions.

Creating a Color Map for Silent Neurons

2008-08-01T06:22:00.000-07:00

We have seen in the previous post that the distribution of membrane potentials among silent neurons sampled using a random perturbation around regular grid is a multimodal distribution, shown below:

There are two main peaks: the tall peak centered around -50mV and the broad bimodal region from about -25mV to -10mV. There is also a small bump, only visible in the log plot (in red) from about -40mV to -25mV. Looking more carefully, we can refine these boundary potentials to -18mv, -31mV, and -42mV and use this to define a color map.

Here is the color map we use (move the slider to view the colors at the far right):

This example also shows how we use convert perturbed grid data into ndvis coordinates, a better approach would be to precompute the integer coordinates (using round) and use those to create a model number which we could index over. This would allow for more rapid access to the data since we wouldn't have to call "round" eight times for each pixel.

Here is the ndvis image for this color map with the "standard" projection na-cat-cas-a for x and kca-kd-h-leak for y:

Some Observations
One thing we see immediately from this picture is that the four regions which were defined by their membrane potential ranges all lie in different regions of parameters space. Lets let D1,D2,D3,D4 represent the four different kinds of silent neurons we have identified according to their membrane potential:

D1: -18mV to 2mV, orange

D2: -31mV to -18mV, blue

D3: -42mV to -31mV, aqua

D4: -63mV to -42mV, green

These domains are defined by neuron model properties. What we have observed is that these domains also map nicely to four regions:

The orange domain D1 (-18mV to 2mV) is concentrated in the region R1 which is the bottom row of each cell in the bottom row of the image, that is, in the region KCa=Kd=0 and it seems a little more dense in the subregions where CaT is large.

The blue domain (-31mV to -18mV) is in the region R2 which is mainly in the lower left cell, where Na=KCa=0, but it seems to extend with linear drop off into NA=1,KCa=0 and NA=1,KCa=1 but concentrated in all cases in the subregion where CaT >=3.

Observe that the blue and orange domains seem to intermix in the intersection of their corresponding regions (R1 and R2) - the two cells on the bottom left with the blue encroaching on the orange region. These two domains are also not cleanly separated in their voltage distributions as they are the two halves of a bimodal distribution. This would lead one to hypothesize that the perhaps the domains ought to defined by a combination of neuron properties and parameter properties. Perhaps there is something fundamentally different about the silent neurons in regions R1 and R2.

The small aqua region is concentrated in the region R3 which is in the lower left corner (Na=KCa=0) and in that
grid cell it is mostly in the upper right (Kd large, CaT small).

The green region can be better visualized with the following dimension order: cas-a-na-h for x and cat-leak-kd-kca for y, as seen in the following image:

This clearly has a "linear" appearance to is and in a later post we will derive a formula for the boundary as a simple hyperplane.

As a final note, let us compare the image we have just analyzed with the image we get when we apply a similar color map to the data from exp1 which used a non-perturbed regularly spaced grid. The non-perturbed data has more peaks in its distribution and these are mostly represented by the red pixels which don't appear in the image from exp3.

Observe that we have roughly the same regions (the green, blue, and orange) which are roughly in the same locations but the perturbed sampling produces an image with a wider range of densities. The non-perturbed has larger solid color bars. We can interpret the densities in the perturbed sampling image as indicating how many neurons have the property of interest in a particular interval represented by that grid cell. If the density is low then that might mean that the region is concentrated in only one small part of the paramter space represented by that grid.

In future posts we will try to find more precise formulas defining the regions R1,R2,R3,R4 of parameter space that correspond to the domains D1,D2,D3,D4 of silent neurson classified by ranges of membrane potentials. We will also look at the behavior of active neurons near the boundaries of these regions.

The Silent/Active Boundary for the PBM neuron model

2008-07-29T13:40:00.000-07:00

Lets now take a look at the boundary between the silent and active neurons in the PBM model using the two data sets discussed earlier: exp1 with a straight grid sampling method and exp3 with a randomly perturbed grid sampling method.

The silent neurons are characterized by one paramter, their resting membrane potential, and this ranges between -100mV and +100mV. Lets look at the distribution plots we get from the two experiments.

First the plot from the fixed grid experiment exp1 is below:

This plot shows a multi-modal distribution with four large regions and two smaller ones at the left. The plot was constructed by grouping neuron models in to bins of size 1mV based on their steadystate membrane potential. The green plot shows the number of neuron models with that potential. The red plot is the log of the distribution function (multiplied by 1000 to fit on the same scale).

Next lets look at the distribution obtained from the random perturbation of parameters around the grid points:

This was created in the same way as the exp1 plot but using the exp3 data. Note that it is a much simpler distribution with only two major modes. The leftmost peak corresponds to a similar peak in the exp1 data, while the rightmost peak corresponds roughly to two exp1 peaks.

The exp3 distribution represents a randomly chosen model but with a distribution that places double the probability on neurons that are within 0.5 of zero in a parameter.

In our book Chapter on silent neurons using the PBM data we performed a detailed analysis of the silent neuron regions based on the strict grid sampling and we identified six distinct regions of the parameter space that corresponded the six main peaks in the distribution. We will not go into that data in detail except to note that all but two of those regions disappear in this new sampling.

This raises the question of whether neuron models with certain parameters exactly equal to zero are interesting when their behavior just slightly away from zero is much different. In other words, is it realistic to assume that some parameters will be exactly zero and not slightly above zero? This is clearly a neuroscience question and may have different answers depending on which parameter is being considered and in which model. For now, we content ourselves with pointing out that the two sampling methods - grid based and perturbed grid based - yield quite different distributions at least for the silent neurons.

We will look at some of the visualizations of silent neurons in the exp3 case in the next post.

Randomly Perturbed Grid-based simulation

2008-07-29T08:25:00.000-07:00

The next major experiment we ran on the cluster was exp3. We ran a little one (exp2) which sampled a 6D grid just to see how long it would take, I may talk about that at some later point...

The goal in exp3 was to explore the consequences of replacing our straight grid based sampling with a different approach in which we randomly perturbed each parameter using a uniform random offset in the range -0.5 .. 0.5 before simulating. Thus the simulation for the grid point
(5.000 3.000 0.000 0.000 0.000 0.000 0.000 2.000)
was
(4.929 3.067 0.132 0.289 0.236 0.328 0.074 2.040)
and note that we take the absolute value after perturbing the zero parameter so all values are positive. The rest of the experiment was identical to exp1. As we will see the results of the exp1 and exp3 runs are very similar for values away from zero, giving us much the same visualizations, but the visualizations near zero can be quite different since zero is a special value in the ODE simulation parameters.

Lets first look at the "basic view" for the randomly perturbed grid model (exp1)

and compare it to the standard grid model from exp1:

Lets make some observations about these two images:

Antialiasing
The randomly perturbed model generates less sharp boundaries especially in the rows and columns corresponding to zero values. Notice for example the second column where na=1, the triangular spiking regions (in yellow) have fuzzier edges in the randomized version which indicates that the the lower right corner (where CaT=5 and Kd=0) is fully inside the single spike region, but in the boundary region where Kd/CaT in in the range [0.5,1] there is a mix of spikers and bursters and when the ration is over 1, there are very few spikers. Thus, this visualization provides a more nuanced look at the relationship between parameter values and simulated behavior.

Behavior near zero easier to understand
If we look at the lower left corner (where Na=KCa=0) we see that the randomized simulation gives a much different picture of what happens in that region. For the strict grid simulation the region was almost entirely silent and there were no bursters at all. In the randomized version we see that there are indeed some bursters (blue) close with Na and KCa in the region [0,0.5] and these tend to occur when Kd>CaT. Morevoer, there are some spikers (yellow) when CaT near 0 and when Kd near 0.

This pattern is repeated in the rest of the first row as well where we see quite a bit of bursting (blue) for all values of Na while the strict grid only showed bursting when Na=1,2 and in a fairly constrained region. Those nice boundaries disappear in the randomized version suggesting that they are an artifact of the behavior of the model for values exactly equal to zero.

In the next post we will look more closely at the boundaries between silent and active neuron models in these two visualization and will argue that the randomized perturbation model gives more useful information, or at least a different kind of information.

Repeating the Prinz,Billimoria,Marder computation

2008-07-28T12:34:00.000-07:00

The first experiment we ran on the cluster was to regenerate part of the Prinz,Billimoria,Marder data.

Programming
We had obtained the C++ code for running a single simulation from Dr. Prinz and so we wrote a main.cpp program that accepted two inputs from the user (a1, a2) in the range 0..5 It would then call the PBM simulator on all possible inputs for the 8 ion channels starting with a1 for Na and a2 for CaT and then using i3,i4,i5,i6,i7,i8 all in the range 0..5 for the remaining channels CaS, A, KCa, Kd, H, Leak which would require simulating 6^6 = 46,656 neurons for each of the 36 (a1,a2) pairs. The simulator was run for 10 simulated seconds (which is 200,000 simulation steps since we simulated with a time step of 1/20,000 second. On our Dell Poweredge 1900 dual quad cores, this took about 1.3 seconds on a single core. We broke the simulation over the entire space of 6^8=1,679,616 models into 36 jobs simulating 6^6 neurons. Each ran for about 10 hours. We ran these jobs on a 2 node subcluster with 16 cores and everything was completed in about 30 hours (corresponding to about 360 compute hours!)

Classification
For this experiment we did not follow the PBM classification criteria exactly but decided to try something simpler though less accurate. We assumed that the neuron would enter steady state aft er 4 simulated seconds and so then started looking for patterns in the last 6 seconds (or 120,000 simulation steps). The simple classification we considered had three components:

silent - if the minimum and maximum voltage over the six simulated seconds are within 5 mV of each other (which ignores small oscillations ....)
tonic spiker - if all of the minima are all within 1 mV of the previous minima, and likewise for the maxima, and likewise for the time interval between successive minima and maxima, that is, if each spike looks very much like the previous spike over all 6 simulated seconds, we way it is a tonic spiker. This lumps together the 1 spike bursters with the PBM tonic spikers, but it is easy to calculate and we wanted to compare the difference. In this case we also stored the length of the last period (i.e. time between the last two maxima).
other - here we assume it is some sort of burster and we return 12/E where E is the number of extrema (maxima and minima) in the last 6 seconds. This is an approximation for the number of spikes per second for the "burster"

Thus, the classification is much less precise than the PBM classification but it has the advantage of being very fast to compute and it generates quite similar results, at least for the "big picture"
For example, we get 289,394 silent neurons where as PBM had 286400. Our classification has many more "spikers" than the PBM because it lumps together the 1 spike bursters with the tonic spikers. We will tease that apart later.

A First Visualization
Here is an example of the standard NDVis view of the data with the order Na, CaT, CaS, A for the x axis and KCa, Kd, H, Leak for the y-axis. We have color coded silent neurons in red, "spikers" in yellow, and "bursters" in blue.

The leftmost column is where Na=0 and we see that this is heavily populated by "spikers" (in yellow) except in the first row which is where Kca=0. That region (Na=KCa=0) is mostly silent.

For those that aren't used to reading NDVis images for this domain. The image consists of a 6x6 grid from which you can read off the value of Na (the column number) and KCa (the row number). Likewise, each of these 36 grid cells is itself a 6x6 grid of smaller cells which correspond to CaT for the x-axis (the columns) and Kd for the y-axis (the rows). Each of these small grid cells is again a 6x6 grid of yet smaller cells from which one can read off CaS on the x and H on the y. Finally the lowest level cell consists of 36 pixels arranged by A on the x-axis and Leak on the y-axis. Changing the dimension order dramatically changes the view and emphasizes different properties of the dataset.

Continuing with our analysis, in the bottom row, where KCa=0, the neurons are mostly spikers except for the bottom line of silents (red) which is where kd=0.

Finally, we observe that almost all of the 25 grid cells where Na>0 and KCa>0 there is a row of spikers at the bottom (where Kd=0) (with KCa=1 having only some spikers).

Thus a rough analysis of this image would lead to the following rough rules for determining whether a neuron on one of these grid points is silent, spiking, or a burster:

KCa=0 and Na=0 => silent 97% (spiker 3%)
KCa=0 and Kd=0 => silent 97% (spiker 3%)

Otherwise,

Na=0 and KCa>0 => spiker 73% (silent 19%)
KCa=0 and Kd>0 and Na>0=> spiker 79% (silent 12%, burster 9%)
KCa>0 and Kd=0 => spiker 66% (silent 26%, burster 8%)

Otherwise,

Na>0 and KCa>0 and Kd>0 => burster 74% (spiker 14%, silent 12%)

The rules for silents are pretty good in this case being 97% accurate, the rules for spikers and bursters range from 66% to 79% accurate, and so indicate strong trends but clearly miss some important behavior.

In future posts we will try to improve the characterizations of these domains, but for now we have demonstrated that the simplified PBM approach gives similar results to the original PBM approach, as far as this first simplest classification goes...

This is admittedly imprecise but it does capture many of the patterns, but we will first move to a slightly different sampling technique (random perturbation of a grid) and examine how that changes the visualization and our conclusions..

Visualizing Grid-based Simulation Results

2008-07-26T08:55:00.001-07:00

Let me start by giving an overview of what I consider the most interesting problems in the NDVis/Neuro project.

The general problem
The general situation we are considering is when someone has developed a mathematical model (usually of some physical or biological process) where the model consists of a collection of coupled non-linear ODEs or more generally some dynamical system which is completely characterized by a vector v at time t and whose behavior at time t+dt can be well approximated by some computational function f(v,t,dt). This computational function approximates a real function F(t) with some initial conditions F(0)=v0.

Moreover, we assume that this system has several parameters p = (p1,p2,...,pk) which naturally occur in some well understood intervals I=(I1,...,Ik). Thus, the computational function should really be indexed by p giving us f_p(v,t,dt) which represents an analytic function F_p(t) with initial conditions F_p(0)=v0.

The general problem is to seek to explain how the behavior of F depends on the parameters p and the initial condition v0. We are particularly interested in models in which F rapidly approaches a periodic function in most cases, (assuming a common starting condition) which is the situation with many neuron models. In this case, there is the additional problem of automatically classifying the types of periodic behavior B(p) that arise for all values of the parameters p in their natural intervals I, and determining the relationship between p and B(p).

The Prinz, Billimoria, Marder approach
For the Prinz, Billimoria, Marder model. Astrid Prinz developed a classification consisting of five types of behavior: silent (i.e. constant voltage), tonically spiking (periodic with a simple sharp spike), bursting (periodic with multiple spikes per period or with a single "fat" spike), and non-periodic within the bounds of the simulation time period. By statistically analyzing the distribution of period times for spiking neurons, Astrid discovered that the period distribution for tonic spikers was bimodal with a clear separation between the peaks for periods of around 0.091 seconds. This allowed her to subdivide the tonic spikers into fast spikers (with periods less that 0.091 seconds) and slow spikers (those with larger periods).

The Prinz,Billimoria,Marder paper studied this problem by running an exhaustive parameter search with 6 values for each of the 8 parameters (Na, CaT, CaS, A, KCa, Kd, H, Leak) and for each possible combination running several simulations. For example, running until a period is found and recording information about the period, then simulating various current injections, and observing the resulting behavior and recording it, etc.

Using their data (which they graciously made publically available), we were able to perform an analysis of the boundaries between the silent and non-silent neurons which you can find in

Visualization and Interactive Exploration of High Dimensional Datasets
by John T. Langton, Elizabeth Gifford, and Timothy Hickey
in "Applications of Computational Intelligence in Biology: Current Trends and Open Problems,", Smolinski, Milanova, Hassanien (eds), Springer, 2008, pp 231-255.

We have published several other papers analyzing that data which you can find at our group research site:

http://neuron.cs.brandeis.edu

This summer we are taking a different tack and using the new high performance cluster at Brandeis (with over 400 cores and over a Terabyte of memory) to generate new simulations of the model and generate new databases.

I will discuss some of our results in the coming blog entries, but for now let me just give an overview of the kinds of studies we're doing.

Random Perturbations of the Sampling Grid
One of the most successful and interesting variations has been to apply a random perturbation in the range [-0.5,0.5] to each parameter before simulating (we also need to take the absolute value of the result to avoid having physically impossible negative parameters for this model). We have found that this approach produces images that are very similar to the straight grid approach but that they are somewhat more informative. In particular they allow one to visually estimate the distance a region boundary is from the fixed parameter values (0,1,2,3,4,5) by looking at the density of points from that region . This is an exciting development and I would recommend that all future grid based studies try this random perturbation approach as well.

The perturbation technique can also allow one to make predictions about the behavior of the model throughout the region rather than just on the simulated subgrid. We have used this approach to study the boundary between silent and non-silent neurons and as we will show the rigid subgrid approach indicates that there are seven types of silent neurons (identified by their membrane potential) and we were able to do a detailed analysis of the structure of those regions in the grid. When we move to a random perturbation model however it turns out that there are only two significant regions and their boundaries are much simpler to describe than in the rigid case. We will devote some of our future postings to these results as well.

The Simulate/Visualize/Resimulate loop
Another variation has been to use the NDVis tool to identify regions of interest and then to resimulate centered around that region and again use the NDVis tool to gain understanding of the model in that region. We have done this with the boundary between fast and slow spikers.

Yet another approach has been to simulate on smaller grids with fewer dimensions. This gives a lower resolution view of the behavior but can be simulate much more quickly, sometimes in minutes rather than hours, allowing for a more interactive exploration of the space.

Other problems we are currently studying
There are several problems we are looking to tackle this summer. One of them is to develop a simple set of formulas to predict the min/max,period of tonic spikers in terms of their parameters. We would like a formula which has an error of less than say 5% in almost all cases, say more than 90% of the time. Another problem is to automatically classify the period behavior we find in the database. We are currently looking at using cubic splines to connect max, min, and inflection points of the membrane potential graph. These can be classified by 5 numbers (the voltage and slope at the endpoints and the distance between the endpoints). Thus a voltage plot can be viewed as a cycle of points in 5 dimensional space and we can try to classify these. There are lots of interesting ideas here but that will have to wait for another post. Finally, we are looking for new ways to find periods or approximately periodic behavior. One problem is that using too strict a criterion can result in a single period that actually consists of two or more almost (but not quite) identical subperiods. Also, for some bursters that have hundreds of spikes per burst, its clear that you have a periodic spiking behavior but the number of spikes may differ slightly for each period and the system might classify it as non-periodic.

The future
In future posts I'll start presenting some of the results mentioned above and discuss more fully some of the interesting questions we're looking into. Although we are currently focusing on the PBM model, the tools we are building will work for any non-linear ODE and it is to develop tools for this more general class of applications that we are studying the PBM model so closely. We hope to expand into some other models in the coming year.

Neuron Simulation

2008-07-17T14:07:00.000-07:00

I will use this blog to describe my experience using John Langton's NDVis tool for Exploratory Data Analysis in several domains. The two domains I am currently actively exploring are a neuron simulation domain and a Flow Cytometry domain.

The Neuron Model domain was first explored by Dr. Astrid Prinz (Emory), Dr. Cyrus Billimoria (BU), and Dr Eve Marder when the first two were in Dr. Marder's lab at Brandeis. We have several publications on the use of the NDVis tool to discover interesting properties of this domain.

The Flow Cytometry domain is a more recent interest and is based on data provided by Dr. Ryan Brinkman at BCCRA http://www.bccrc.ca/tfl/people_rbrinkman.html

This blog is meant to be read mainly by my collaborators and students in my lab, so outside readers may need to read some of the papers on my website http://www.cs.brandeis.edu/~tim to make sense of what we are doing here...