Automatic Fitting of Spiking Neuron Models to Electrophysiological Recordings

Spiking models can accurately predict the spike trains produced by cortical neurons in response to somatically injected currents. Since the specific characteristics of the model depend on the neuron, a computational method is required to fit models to electrophysiological recordings. The fitting procedure can be very time consuming both in terms of computer simulations and in terms of code writing. We present algorithms to fit spiking models to electrophysiological data (time-varying input and spike trains) that can run in parallel on graphics processing units (GPUs). The model fitting library is interfaced with Brian, a neural network simulator in Python. If a GPU is present it uses just-in-time compilation to translate model equations into optimized code. Arbitrary models can then be defined at script level and run on the graphics card. This tool can be used to obtain empirically validated spiking models of neurons in various systems. We demonstrate its use on public data from the INCF Quantitative Single-Neuron Modeling 2009 competition by comparing the performance of a number of neuron spiking models.

There are several important diffi culties in this model fi tting problem.Firstly, because of the threshold property, the mapping from a time-varying signal to spike trains is generally discontinuous in spiking models (Brette, 2004).Besides, the fi tness criterion we used (the gamma coincidence factor used in the INCF competition; Jolivet et al., 2008) is discrete, because it is a function of the number of coincidences within a predefi ned temporal window.These facts prevent us from using optimization methods based on gradients (e.g.conjugated gradient descent).Secondly, a single evaluation of the criterion for a given set of parameter values involves the simulation of a neuron model over a very long time.For example, recordings in challenges A and B of the 2009 INCF competition last 60 s, sampled at 10 kHz, totalling 600,000 values of the input signal.Thus, evaluating the fi tness criterion for any spiking model involves several millions of operations.Thirdly, not only the parameter values are unknown, but there are also many candidate models.In particular, it is not unreasonable to think that different neuron types may be best fi t by different phenomenological models.Therefore, the optimization tools should be fl exible enough to allow testing different models.
To address these issues, we developed a model fi tting toolbox for the spiking neural network Brian (Goodman and Brette, 2009).Brian is a simulator written in Python that lets the user defi ne a model by directly providing its equations in mathematical form (including threshold condition and reset operations).Using an interpreted language such as Python comes at a cost, because the interpretation overhead can slow down the simulations, but this problem can be solved by vectorizing all operations when the network model includes many neurons.It turns out that the same strategy applies to programming GPUs, which are most effi cient when all processors execute the same operation (the Single Instruction, Multiple Data (SIMD) model of parallel programming).Therefore, we developed several vectorization techniques for spiking model optimization (see Vectorization

INTRODUCTION
Neurons encode time-varying signals into trains of precisely timed spikes (Mainen and Sejnowski, 1995;Brette and Guigon, 2003), using a diverse set of ionic channels with specifi c characteristics.Recently, it was found that simple phenomenological spiking models, such as integrate-and-fi re models with adaptation, can in fact predict the response of cortical neurons to somatically injected currents with surprising accuracy in spike timing (Jolivet et al., 2004;Brette and Gerstner, 2005;Gerstner and Naud, 2009).This unexpected performance is probably related to the fact that detailed conductance-based models with widely diverse ion channel characteristics can in fact have the same properties at neuron level (Goldman et al., 2001).These encouraging results triggered an interest in quantitative fi tting of neuron models to experimental recordings, as assessed by the recent INCF Quantitative Single-Neuron Modeling competition.The competition has seen several successful submissions, but there is no available method to systematically fi t arbitrary models to experimental data.Such computational tools would be particularly useful for modellers in systems neuroscience, for example, who could use empirically validated models in their studies.We developed a model fi tting library, which works with the Brian simulator (Goodman and Brette, 2009), and allows the fi tting of user-defi ned spiking models to electrophysiological data consisting of spike trains elicited by time-varying signals (for example, intracellularly injected currents).If the user's machine has a graphics processing unit (GPU) − a cheap parallel processor available on many standard machines − the algorithms run in parallel on it.The library is available as a part of Brian, which is distributed under a free open-source license1 .
Spiking models can accurately predict the spike trains produced by cortical neurons in response to somatically injected currents.Since the specifi c characteristics of the model depend on the neuron, a computational method is required to fi t models to electrophysiological recordings.The fi tting procedure can be very time consuming both in terms of computer simulations and in terms of code writing.We present algorithms to fi t spiking models to electrophysiological data (time-varying input and spike trains) that can run in parallel on graphics processing units (GPUs).The model fi tting library is interfaced with Brian, a neural network simulator in Python.If a GPU is present it uses just-in-time compilation to translate model equations into optimized code.Arbitrary models can then be defi ned at script level and run on the graphics card.This tool can be used to obtain empirically validated spiking models of neurons in various systems.We demonstrate its use on public data from the INCF Quantitative Single-Neuron Modeling 2009 competition by comparing the performance of a number of neuron spiking models.
Since the simulations are independent within one iteration, they can be run simultaneously.In the PSO algorithm, each particle accelerates towards a mixture of the location of the best particle and the best previous location of that particle.The state update rule combines deterministic and stochastic terms in order to prevent the particles from getting stuck in locally optimal positions: where X i (t) and V i (t) are the position and speed vectors of particle i at time t, respectively.X t i l ( ) is the best position occupied by particle i before time t (local best position), and X g (t) is the best position occupied by any particle before time t (global best position).ω, c l and c g are three positive constants which are commonly chosen as follows in the literature (Shi and Eberhart, 1998;den Bergh, 2006): ω = 0.9, c l = c g = 1.9.We chose these values for most results shown in this paper.For the in vitro recordings, we chose c l = 0.1 and c g = 1.5 which we empirically found to be more effi cient.Finally, r l and r g are two independent random numbers uniformly resampled between 0.0 and 1.0 at each iteration.
Boundary constraints can be specifi ed by the user so that the particles are forced to stay within physiologically plausible values (e.g.positive time constants).The initial values V i (0) are set at 0, whereas the initial positions of the particles are uniformly sampled within user-specifi ed parameter intervals.The convergence rate of the optimization algorithm decreases when the interval sizes increase, but this effect is less important when the number of particles is very large.

Fitness function.
The computation of the gamma factor can be performed in an offl ine or online fashion.The offl ine method consists in recording the whole spike trains and counting the number of coincidences at the end of the simulation.Since the fi tness function is evaluated simultaneously on a very large number of neurons, this method can be both memory-consuming (a large number of spike trains must be recorded) and time-consuming (the offl ine computation of the gamma factor of the model spike trains − which can consist of several hundred thousand spikes − against their corresponding target trains is performed in series).The online method consists in counting coincidences as the simulation runs.It allows us to avoid recording all the spikes, and is much faster than the offl ine algorithm.Moreover, on the GPU the online algorithm requires a GPU to CPU data transfer of only O(num particles) bytes rather than O(num spikes), and CPU/GPU data transfers are a major bottleneck.Our fi tting library implements the online algorithm.About 10% of the simulation time at each iteration is spent counting coincidences with our algorithm on the CPU.

Vectorization over data sets
The particle swarm algorithm can be easily adapted so that a single neuron model can be fi tted to several target spike trains simultaneously.The algorithm returns as many parameter sets as target spike trains.The gamma factor values for all particles and all spike trains are also computed in a single run, which is much faster than computing the gamma factors in series.For example, in the INCF Quantitative Modeling dataset, all trials (13 in challenge A and 9 in challenge B) can be simultaneously optimized.
Techniques).These techniques apply both to CPU simulations and to parallel GPU simulations.We used just-in-time compilation techniques to keep the same level of fl exibility when models are simulated on the GPU, so that using the GPU is transparent to the user (see GPU Implementation).We demonstrate our algorithms (see Results) by fi tting various spiking models to the INCF competition data (challenges A and B).Consistent with previous studies, we found that adaptive spiking models performed very well.We also show how our tool may be used to reduce complex conductance-based models to simpler phenomenological spiking models.

VECTORIZATION TECHNIQUES
Fitting a spiking neuron model to electrophysiological data is performed by maximizing a fi tness function measuring the adequacy of the model to the data.We used the gamma factor (Jolivet et al., 2008), which is based on the number of coincidences between the model spikes and the experimentally-recorded spikes, defi ned as the number of spikes in the experimental train such that there is at least one spike in the model train within ±δ, where δ is the size of the temporal window (typically a few milliseconds).The gamma factor is defi ned by , where N coinc is the number of coincidences, N exp and N model are the number of spikes in the experimental and model spike trains, respectively, and r exp is the average fi ring rate of the experimental train.The term 2δN exp r exp is the expected number of coincidences with a Poisson process with the same rate as the experimental spike train, so that Γ = 0 means that the model performs no better than chance.The normalization factor is chosen such that Γ ≤ 1, and Γ = 1 corresponds to a perfect match.The gamma factor depends on the temporal window size parameter δ (it increases with it).However, we observed empirically that the parameter values resulting from the fi tting procedure did not seem to depend critically on the choice of δ, as long as it is not so small as to yield very few coincidences.For most results shown in the Section "Results", we chose δ = 4 ms as in the INCF competition.This fi tting problem can be solved with any global optimization algorithm that does not directly use gradient information.These algorithms are usually computationally intensive, because the fi tness function has to be evaluated on a very large number of parameter values.We implemented several vectorization techniques in order to make the fi tting procedure feasible in a reasonable amount of time.Vectorization allows us to use the Brian simulator with maximum effi ciency (it relies on vectorization to minimize the interpretation overhead) and to run the optimization algorithm in parallel.

Vectorization over parameters
Particle swarm optimization algorithm.We chose the particle swarm optimization (PSO) algorithm (Kennedy and Eberhart, 1995), which involves defi ning a set of particles (corresponding to parameter values) and letting them evolve in parameter space towards optimal values (see Figures 1 and 7).Evaluating the fi tness of a particle requires us to simulate a spiking neuron model with a given set of parameter values and to calculate the gamma factor.

Vectorization over time
The gamma factor is computed by simulating the neurons over the duration of the electrophysiological recordings.The recordings can be as long as several tens of seconds, so this may be a bottleneck for long recordings.We propose vectorizing the simulations over time by dividing the recording into equally long slices and simulating each neuron simultaneously in all time slices.Spike coincidences are counted independently over time slices, then added up at the end of the simulation when computing the gamma factor.The problem with this technique is that the initial value of the model at the start of a slice is unknown, except for the fi rst slice, because it depends on the previous stimulation.To solve this problem, we allow the time windows to overlap by a few hundreds of milliseconds and we discard the initial segment in each slice (except the fi rst one).The initial value in the initial segment is set at the rest value.Because spike timing is reliable in spiking models with fl uctuating inputs (Brette and Guigon, 2003), as in cortical neurons in vitro (Mainen and Sejnowski, 1995), spike trains are correct after the initial window, that is, they match the spike trains obtained by simulating the model in a single pass, as shown in Figure 2.
The duration of one slice is overlap + recording duration/number of slices, so that the extra simulation time is overlap × number of slices, which is, relative to the total simulation time, overlap/slice duration.Thus, there is a trade off between the overhead of simulating overlapping windows and the gain due to vectorisation.In our simulations, we used slices that were a few seconds long.The duration of the required overlap is related to the largest time constant in the model.

GPU IMPLEMENTATION
A graphics processing unit (GPU) is a type of chip available on modern graphics cards.These are inexpensive units designed originally and primarily for computer games, which are increasingly being used for non-graphical parallel computing (Owens et al., 2007).
The chips contain multiple processor cores (240 in the current state of the art designs, and 512 in the next generation which will be available in 2010) with a limited ability to communicate between each other.This makes them ideal for the particle swarm algorithm where many independent simulations need to be run for each iteration of the algorithm, especially with the vectorization techniques that we presented in the Section "Vectorization Techniques".Programming for a GPU is rather specialised.Each processor core on a GPU is much simpler than a typical CPU, and this places considerable limitations on what programs can be written for them.Moreover, although recent versions of these chips allow more of the functionality of a full CPU (such as conditional branching), algorithms that do not take into account the architecture of the GPU will not use it effi ciently.In particular, the GPU places constraints on memory access patterns.Consequently, although 240 cores may be present in the GPU, it is unrealistic in most cases to expect a 240× speed increase over a CPU.However, speed improvements of tens to hundreds of times are often possible (see the showcase on the NVIDIA CUDA Zone)2 .In our case, we have achieved a roughly 50−80× speed improvement, thanks largely to the fact that the model fi tting algorithm is "embarassingly parallel", that is, that it features a number of independent processes which do not interact with each other.This means that memory can be allocated in the topologically continuous fashion that is optimal for the GPU, and the code does not need to introduce synchronization points between the processes which would slow it down.
We use the GPU in the following way.We have a template program written in the subset of C++ that runs on the GPU.Code that is specifi c to the neuron model is generated and inserted into the template when the model fi tting program is launched.The code is generated automatically from the differential equations defi ning the neuron model, compiled and run automatically on the GPU (just-in-time compilation) using the PyCUDA package (Klöckner et al., 2009).
The code generation works as follows.We start from a basic template that handles standard operations shared amongst all models, including loading the input currents to the model, computing the coincidence count and so forth.The template has several slots for specifying code specifi c to the model, including the numerical integration step and the thresholding and reset behavior.The differential equations defi ning the neuron model are stored as strings in Brian.For example, if the differential equation was 'dV/dt=-V/tau' then we store the substring '-V/tau'.An Euler solver template for this might look like: for var, expr in diffeq: (var)__tmp = (expr); for var in diffeq: (var) += dt*(var)__tmp; Vm Time (s) A leaky integrate-and fi re neuron receiving a dynamic input current is simulated over 2.4 s (top row, black).The same model is simultaneously simulated in three time slices (blue) with overlapping windows (dashed lines).Vertical dotted lines show the slice boundaries.The model is normalized so that threshold is 1 and reset is 0. The overlap is used to solve the problem of choosing the initial condition and is discarded when the slices are merged.
In the case above this would become (after substitutions): This latter code snippet is valid C code that can be run on the GPU.A similar technique is used for thresholding and reset behaviour.
These techniques mean that the huge speed improvements of the GPU are available to users without them having to know anything about GPU architectures and programming.Usage is transparent: on a system with a GPU present and the PyCUDA package installed, the model fi tting code will automatically be run on the GPU, otherwise it will run on the CPU.

DISTRIBUTED COMPUTING
As mentioned in the previous section, the model fi tting algorithm presented here is close to "embarassingly parallel", meaning that it can relatively easily be distributed across multiple processors or computers.In fact, the algorithm consists of multiple iterations, each of which is "embarassingly parallel" followed by a very simple computation which is not.More precisely, the simulation of the neuron model, the computation of the coincidence count, and most of the particle swarm algorithm can be distributed across several processors without needing communication between them.The only information that needs to be communicated across processors is the set of best param-eters found so far.With N particles and M processors, each processor can independently work on N/M particles, communicating its own best set of parameters to all the other processors only at the end of each iteration.This minimal exchange of information means that the work can be effi ciently distributed across several processors in a single machine, or across multiple machines connected over a local network or even over the internet.We use the standard Python multiprocessing package to achieve this.See the section titled "Benchmark" and Figure 6 for details on performance.

MODEL FITTING
We checked our algorithms on synthetic data and applied them to intracellular in vitro recordings.

Synthetic data
We fi rst check that the optimization algorithms can fi t a spiking neuron model to spike trains generated by the same model (Figure 3).We simulated a leaky integrate-and-fi re neuron model with an adaptive threshold (equations in Table 1) responding to a fl uctuating current injection over 500 ms (Ornstein-Uhlenbeck process).As expected, the optimization algorithm converged to a perfect fi t (Γ = 1 with precision δ = 0.1 ms) after just a few iterations.We chose this small value of δ to show that the fi t can be nearly perfect, as expected in from the optimization algorithm match the original parameters within 15% error and the resulting trace (blue; threshold in green; the model is normalized so that V m has no unit) is close to the original trace (black), even though only spike timings were used for fi tting.The gamma factor was 1.0 at precision δ = 0.1 ms.
this situation.The corresponding parameter values were very close to the original ones (±15%), so that the original and optimized traces match, even though only the spike trains were used for the fi tting procedure (a perfect match of the traces can be obtained by using more particles than we used for the fi gure, giving parameters values within ±3%).However, this is not likely to be a general result: for example, if the neuron model included after-spike effects (such as refractoriness) and the interspike intervals were longer than these effects (that is, at low fi ring rates), then the method could not recover these properties since they would not be visible in the spike trains.

In vitro recordings
We then applied our fi tting procedure to in vitro intracellular recordings from the 2009 Quantitative Single-Neuron Modeling competition (challenges A and B).In these recordings, fl uctuating currents were injected into the soma of cortical neurons (L5 regular spiking pyramidal cell for challenge A and L5 fast spiking cell for challenge B). Figure 4 shows the result of fi tting an integrate-and-fi re model with adaptive threshold (equations in Table 1) to an intracellular recording of an L5 fast spiking cell responding to in-vivo-like current injection (Challenge B).For the 1 s sample shown in the fi gure, the optimization algorithm converged in a few iterations towards a very good fi tness value (Γ = 0.9 at precision δ = 2 ms).Tables 1 and 2 report the results of fi tting four different spiking models to the recordings in challenges A (regular spiking cell) and B (fast spiking cell).Each challenge includes the recordings of several identical trials (same input current; 13 trials in A and 9 trials in B), and we report the average and standard deviation for all quantities (gamma factor and parameter values).The recordings were divided into a training period (fi rst 10 s), which was used to optimize the models, and a test period (last 10 s), which was used to calculate the gamma factors reported in the tables.We used the following values for the optimization algorithm: ω = 0.9, c l = 0.1, c g = 1.5.We used a local constant Table 1 | Optimization results for Challenge A. Four neuron models were fi tted to the electrophysiological recordings of a regular spiking L5 pyramidal cell responding to in-vivo-like current injection.There were 13 trials with the same input.The models were optimized on the part of the data between 17.5 seconds and 28 seconds and the gamma factors were calculated between 28 and 38 seconds (δ = 4 ms).The value relative to the intrinsic reliability Γ in is reported in brackets.Parameter D is a time shift for output spikes (recorded spikes typically occur slightly after model spikes because they are reported at the time when the membrane potential crosses 0 mV).We also reported rescaled versions of the parameter values (in brackets) so that they correspond to electrophysiological quantities.

ADAPTIVE THRESHOLD IF
τdV /dt = RI − V 0.54 ± 0.05 (69%) R = 3.85 ± 1.13 /nA (56.5 ± 2.4 MΩ) value much smaller than the global one so that the evolution of the algorithm is not dominated by the local term, which would make optimization slower.Since only the spike trains were used for fi tting, units were arbitrary (e.g.reset is 0 and threshold is 1 for the adaptive IF model).To interpret the parameter values, we also report scaled versions of the parameters obtained by changing the voltage units of the model in a such a way that the average and standard deviation of the model trace agree with those of the intracellular recording.The mean µ d and the standard deviation σ d of the experimental voltage trace were computed over the test period after the action potentials were cut.Then, the mean µ m and the standard deviation σ m of the voltage traces of the fi tted model were also computed over the same test period.Finally, each parameter value X of the model was affi ne-transformed in either of the following two ways: Interestingly, for both challenges, best performance was achieved by a simple model, the integrate-and-fi re model with adaptive threshold: Γ = 0.54 ± 0.05 for challenge A and Γ = 0.76 ± 0.07 for challenge B. The intrinsic reliability of the neurons can be defi ned as the average gamma factor between trials (Γ in = 0.78 for challenge A and Γ in = 0.74 for challenge B).Relative to the intrinsic reliability, the performance was 69% and 102% respectively.The fact that the performance on challenge B is greater than 100% probably refl ects a drift in parameter values with successive trials, most likely because the cell was damaged.Indeed, the performance was tested on data that was not used for parameter fi tting, so it could not refl ect overfi tting to noisy recordings but rather a systematic change in neuron properties.Consistently with this explanation, in both challenges the fi ring rate increased over successive trials (whereas the input was exactly identical).This change in neuronal properties was likely caused by dialysis of the cell by the patch electrode (the intracellular medium is slowly replaced by the electrode solution) or by cell damage.
A related model with an adaptive current instead of an adaptative threshold performed similarly.The performance of the adaptive exponential integrate-and-fi re model (Brette and Gerstner, 2005) (AdEx), which includes a more realistic spike initiation current (Fourcaud-Trocme et al., 2003), was not better.This surprising result is explained by the fact that the optimized slope factor parameter (Δ T ) was very small, in fact almost 0 mV, meaning that spike initiation was as sharp as in a standard integrate-and-fi re model.The Izhikevich (2003)

MODEL REDUCTION
Our model fi tting tools can also be used to reduce a complex conductance-based model to a simpler phenomenological one, by fi tting the simple model to the spike train generated by the complex model in response to a fl uctuating input.We show an example of this technique in Figure 5 where a complex conductance-based model described in benchmark 3 of Brette et al. (2007) is reduced to an adaptive exponential integrate-and-fi re model (Brette and Gerstner, 2005).The complex model consists of a membrane equation and three Hodgkin-Huxley-type differential equations describing the dynamics of the spike generating sodium channel (m and h) and of the potassium rectifi er channel (n).
In this example, the gamma factor was 0.79 at precision 0.5 ms.This technique can in principle be applied to any pair of neuron models.

BENCHMARK
To test the scaling performance across multiple processors, we used a three machine cluster connected over a Windows network.Each computer consisted of a quad-core 64 bit Intel i7 920 processor at 2.6 GHz, 6 GB RAM, and an NVIDIA GeForce GTX 295 graphics card (which have two GPUs).The cluster as a whole then had 12 cores and 6 GPUs.
Performance scaled approximately linearly with the number of processors, either with CPUs (Figure 6A) or GPUs (Figure 6B).In the case of CPUs, performance was close to ideal (that is N processors performing N times faster than a single processor), and was slightly lower than ideal with GPUs.With our relatively inexpensive three machine cluster, we could achieve performance approximately 300 times faster than with a single CPU, allowing us to fi t models in hours which would previously have taken weeks.A single GPU performed approximately 65× faster than a single CPU (or around 16× faster than the four cores available on a single machine in our cluster), and a dual-GPU GTX 295 card performed around 108× faster than a single CPU (or around 27× faster than the four CPUs alone).

DISCUSSION
We presented vectorized algorithms for fi tting arbitrary spiking neuron models to electrophysiological data.These algorithms can run in parallel on a graphics processing unit (GPU) or on multiple cores.It appeared that the speed improved by a factor of 50−80 times when the GPU was used for model simulations.With three dual-GPU cards, the performance was about 300 times faster than with one CPU, which makes it a cheap alternative to clusters of PCs.The algorithms are included as a model fi tting library for the Brian simulator (Goodman and Brette, 2009), which is distributed under a free open-source license3 .This computational tool can be used by modellers in systems neuroscience, for example, to obtain empirically validated models for their studies.
We chose to use the particle swarm algorithm for optimization, but it can be easily replaced by any other global optimization algorithm that uses simultaneous evaluations of different parameter sets, such as genetic algorithms.Indeed, the optimization procedure is defi ned at script level (in Python) and runs on the main processor.The error criterion could also be modifi ed, for example to include an error on the intracellular voltage trace, so that the model can predict both spike times and voltage.
Other model fitting techniques have been previously described by several authors, most of them based on maximum likelihood (Paninski et al., 2004(Paninski et al., , 2007)).Our initial motivation for choosing a more direct approach based on general global optimization methods was that it applies to arbitrary models, including nonlinear ones, whereas maximum likelihood optimization is generally model-specific.Current maximum likelihood A fi tting task was performed over 1-s-long electrophysiological data.(A) Speedup in the simulation of 10 iterations with 400,000 particles as a function of the number of cores, relative to the performance of a single CPU.We used different cluster confi gurations using 1−12 cores spread among one up to three different quad-core computers connected over a local Windows network.(B) Speedup in the simulation of 10 iterations with 2,000,000 particles as a function of the number of GPUs, relative to the performance of a single CPU (using the same number of particles).GPUs were spread on up to three PCs with dual-GPU cards (GTX 295).
techniques apply essentially to linear threshold models, which constitute a large class of models but do not include, for example, the AdEx model and Izhikevich model, which we could evaluate with our algorithm.In Huys et al. (2006), a sophisticated maximum likelihood method is presented to estimate parameters of complex biophysical models, but many aspects of the model must be known in advance, such as time constants and channel properties (besides, the voltage trace is also used).Another important difference is in the application of these techniques.The motivation for maximum likelihood approaches is to determine parameter values when there is substantial variability in the neural response to repeated presentations of the Our results on challenges A and B of the INCF Quantitative Modeling competition confi rm that integrate-and-fi re models with adaptation give good results in terms of prediction of cortical spike trains.Interestingly, the adaptive exponential integrate-andfi re model (Brette and Gerstner, 2005) did not give better results although spike initiation is more realistic (Badel et al., 2008).It appears that the fi tting procedure yields a very small value for the slope factor parameter Δ T , consistent with the fact that spikes are sharp at the soma (Naundorf et al., 2006;McCormick et al., 2007).The Izhikevich model also did not appear to fi t the data very well.This could be because spike initiation is not sharp enough in this model (Fourcaud-Trocme et al., 2003) or because it is based on the quadratic model, which approximates the response of conductancebased models to constant currents near threshold, while the recorded neurons were driven by current fl uctuations.
Our technique can also be used to obtain simplifi ed phenomenological models from complex conductance-based ones.Although it primarily applies to neural responses to intracellular current injection, it could in principle be applied also to extracellularly recorded responses to time-varying stimuli (e.g.auditory stimuli), if spike timing is reproducible enough.An interesting extension, which could apply to the study of phase locking properties of auditory neurons, would be to predict the distribution of spike timing using stochastic spiking models.

ACKNOWLEDGMENTS
This work was partially supported by the European Research Council (ERC StG 240132).

FIGURE 1 |
FIGURE 1 | Fitting spiking neuron models with particle swarm optimization (PSO).The target spike train was generated by injecting a noisy current (white noise, 500 ms long) into a leaky integrate-and-fi re model (R = 10 10 A −1 and τ = 40 ms), and PSO was run to fi nd the correct parameter values.The boundary conditions used for the optimization are the ranges of the axes.The optimization algorithm makes the particles (sets of parameter values) evolve towards the area with high fi tness values (an oblique line in this example).(A) Positions of the particles in parameter space (R, τ) at the start of the algorithm and their evolution at the next iteration (arrows).The colored background represents the value of the gamma factor for all parameter values.(B) Voltage trace of the red particle shown in A compared to the original one (blue), in response to the injected current (top).The y-axis is unitless.(C) Positions of the particles at iteration 4. (D) Voltage trace of the red particle (same as in (A) and (B)) at iteration 4. The gamma factor increased from 0.09 to 0.55 between (B) and (D).

FIGURE 4 |
FIGURE 4 | Fitting a leaky integrate-and-fi re neuron with adaptive threshold to electrophysiological recordings.(A) A fl uctuating in-vivo-like current was injected intracellularly in the soma of a L5 fast spiking cell.(B) The

FIGURE 5 |FIGURE 6 |
FIGURE 5 | Model reduction from a Hodgkin−Huxley model to an adaptive exponential integrate-and-fi re model (AdEx).(A) A 500-ms fl uctuating current is injected into a conductance-based model.(B) The output spike train is used to fi t an AdEx model.(C) Voltage trace of the optimized AdEx model.The gamma factor is 0.79 at precision 0.5 ms and the parameter values are τ = 12 ms, τ w = 25 ms, R = 7.0 × 10 9 A −1 , V R = − 0.78, Δ T = 1.2, a = 0.079, α = 0 (normalized voltage units).

FIGURE 7 |
FIGURE 7 | Evolution of the fi tness values of particles during optimization.The color-coded background shows the distribution of the fi tness values (Γ) of the 10,000 particles in each of the 500 iterations of the fi tting procedure for trial 3 of challenge B (darker means higher density).The dotted line shows the gamma factor of the best particle.The neuron model was the adaptive threshold IF model.The best gamma factor reaches a plateau after 200 iterations but the other particles continue to evolve.