Fitting a recurrent dynamical neural network to neural spiking data: Tackling with the sigmoidal gain function issues
Resat Ozgur Doruk


Neural Spiking and System Identification]Fitting a recurrent dynamical neural network to neural spiking data: Tackling with the sigmoidal gain function issues



Department of Electrical and Electronic Engineering, Atılım University, İncek, Gölbaşı, Ankara, TURKEY
resat.doruk@atilim.edu.tr




  July 2018



This is a continuation of a recent study on the modeling of the information coding in sensory system in the brain. The data from a sensory neurons are available as discrete spike timings with no amplitude information. In the simulations, these are generated from a data generator model which has certain differences from the model being estimated. The model under consideration is simpler than the one used as a data generator as it has no sigmoidal gain function parameters. This choice is based on a suggestion from a recent study which states that inclusion of gain functions to the estimation algorithm increases issues such as parameter confounding which leads to performance degradation issues like increased or unpredictable variance changes with different stimuli configurations. To resolve this issue we consider a more generic model that has no sigmoidal gain functions to be estimated. Like that of the first research, we applied a Fourier series stimulus to both data generator and the estimated network. In order to test the performance of the proposed scheme, we have repeated the simulations with different sample sizes, stimulus component counts, amplitude and base frequencies. Mean values of estimation are presented as tables and the statistical analysis results are presented in graphical forms. In addition, since we do not have any true parameter data (for the estimated model) we compare the firing rate responses of both data generator and estimated models to different stimuli. It appeared that the responses of both models are almost the same. 




 

 Keywords: Neural spiking, Maximum likelihood estimation, Sigmoidal gain functions, Parameter confounding













§ INTRODUCTION



 §.§ General Neuron Modeling


Neuron modeling is an important milestone in the theoretical or computational neuroscience fields. The nonlinear biological oscillator Hodgkin-Huxley model <cit.> is known as a sparking study that triggers the derivation of further simplified or complicated models. Some models only describe the dynamical behaviour of the neurons in consideration. <cit.> models can be considered as an example where the bursting behaviors of a biological neuron is being considered. On the contrary, models including <cit.> involve both the dynamical behavior of the membrane potential and some of the internal dynamics due to sodium, potassium and calcium channels and also leakage of chlorine ions.      

In addition there are also some intermediate models such as <cit.>. In those, a linear filter and nonlinear cascade is combined. They focus on the computational aspects of a neuron (dynamical features). There are also so called "blackbox models" which are basically a statistical relationship that defines the statistics of the response against a given stimulus. Related examples are <cit.>. 

Another option in neural modeling is the utilization of modified generic neural networks <cit.>. These models may be either of static feedforward <cit.> or recurrent dynamical type <cit.>. The generic neural networks are a good option if we are interested in the information coding rather than the biophysical components. In addition, the recurrent neural networks are known to have universal approximation properties <cit.>. Because of that, they are eligible to be used in the modeling of different neurological functions <cit.>. 



 §.§ Neural Spiking Issue
 

Most of the biological neurons transmit information through a phenomenon called as neural spiking <cit.>. That is a series of action potentials fired by a biological neuron against a stimulus. As repeated action potentials (or bursts) have the same shape, there should be a mechanism that plays a role in  information coding. In sensory neurons, the temporal locations and rates of the spiking (rate of firing of action potentials) are believed to code the transmitted information <cit.>. Recent studies also showed that, the neural spiking is not a deterministic process due to stochasticity in ion channels (sodium, potassium, clacium etc.) and synapses <cit.>. In addition this stochasticity is found to obey an inhomogeneous Poisson distribution <cit.> with neural firing rate as its event rate <cit.>. There are numerous studies utilizing the findings of that research.  <cit.> develops a static feed-forward network model for the tonotopic auditory cortex in marmoset monkeys. These models are rather computational which are dealing with the information processing side of the neurons. In <cit.>, one can see another model where a continuous time recurrent neural network is considered in the same fashion. In that paper, stimulus is modeled as a Phased Cosine Fourier Series (PCFS). The advantage of last paper over works like <cit.> is that, the time dependent features of the model is taken into consideration. In any case, the main output of the models is the neural firing rate which represented by a set of spikes recorded.       



 §.§ About This Work
 


The research is to accompany a recent study where the parameters of a continuous time nonlinear dynamical recurrent neural network is fitted from neural spiking data <cit.>. However there were certain issues which are not addressed. One major peculiarity is the assumption that the parameters associated with the sigmoidal gain functions are known. In reality, these parameters are not known and should be estimated. However, another recent study <cit.> suggests that parameter confounding issues lead to very inefficient estimates. For example in <cit.>, the variance of the estimates are decreasing with the increasing sample sizes (or number of iterations). Contrary to that, the work in <cit.> reveals that inclusion of the gain parameters into the evaluation leads to unrelated results (i.e. no definite pattern versus the sample size). So the problem should be reconsidered so that one has a useful generic model with no presumed parameters. In addition, if only the firing rate and the associated spike timings are topics of interest one can try a generic continuous time recurrent neural network without any tunable gain function parameters <cit.> in the identification procedures. The parameter estimation procedure is similar to that of <cit.> and based on joint maximum likelihood estimation <cit.>. 



§ MATERIALS AND METHODS



 §.§ Neural Model in Consideration

 
Like in <cit.>, the model will also be a recurrent dynamical neural network <cit.> with excitatory (e) and inhibitory (i) components:

    ẋ_e   =-β_e [  x_e+w_eeg(x_e)-w_eig(x_i)+c_eu ]  
    ẋ_i   =-β_i [  x_i+w_ieg(x_e)-w_iig(x_i)+c_iu ]
 
with g(x) being a generic sigmoid function as shown below:

    g(x)=1/1+e^-α x
     
where α is a small scaling constant. It can be as small as α=0.001. In this model setting, x_e and x_i are dimensionless variables representing the dynamics of the excitatory and inhibitory components of the model in (<ref>) respectively. This is a small difference from <cit.> as the state variables have units in mVs. The input u, represents a generic stimulus which can be modeled as a function of time (such as Fourier Series).  

Like in <cit.> the output of the model will be the excitatory firing rate r_e. It can be related to the dynamics of the excitatory neuron by a sigmoidal relationship as shown:


    r_e=F_e/1+exp(-α x_e)


The above is similar to (<ref>) with a small difference which is the maximum firing rate parameter F_e. 
As a difference from the model setting in <cit.>, the firing b The parameters to be estimated are θ=[β_e,β_i,w_e,w_i,w_ee,w_ei,w_ie,w_ii,F_e] and their definitions are given in Table <ref> below:



The estimation procedure will need a data generator that will generate the required neural spiking data sets to be used in computation. The data should represent a realistic neural spiking event set. As this will be a theoretical study, one can achieve this goal by simulating the same or another neuron model and obtain a set of useful spiking data. Obtaining the data from (<ref>) with true values of parameters θ will not be considered here as a similar situation already exists in <cit.>. In addition, the true values of parameters θ will actually be not known. Thus, one can generate the neural spiking data from a different model and train the network in (<ref>). The details will be discussed in Section <ref>.



 §.§ Neural Spiking and Poisson Processes

 


  §.§.§ Inhomogeneous Poisson Processes
  
It is stated in Section <ref> that neural spiking in sensory neurons largely obeys an Inhomogeneous Poisson Distribution <cit.> which can be expressed as shown below:

    [N(t+Δ t)-N(t)=k]=e^-λλ^k/k!

where

    λ=∫_t^t+Δ tr_e(τ)dτ

is the average number of spikes based on the firing rate r_e(t) of the excitatory unit and N(τ) indicates the cumulative total number of spikes up to time τ, so that  N(t+Δ t)-N(t) is the number of spikes within the time interval [t,t+Δ t). 

Another interpretation is that, the probability of observing k spikes in the time interval (t,t+Δ t) will be given by (<ref>). 

Suppose that one has K spikes (t_1,t_2,…,t_K) occured in the interval (0, T) (here 0≤ t_1 ≤ t_2 ≤…≤ t_K ≤ T so t and Δ t become t=0 and Δ t=T). 

The probability density function for a given spiking train (t_1,t_2,…,t_K) can be derived from the inhomogeneous Poisson process <cit.> as shown in the following:


    p(t_1,t_2,…,t_K)
    =exp(-∫_0^Tr_e(t, u,θ)dt)∏_k=1^K r_e(t_k,u,θ)


The above formulation describes how likely a particular spike train (t_1,t_2,…,t_K) is generated by an inhomogeneous Poisson process with event rate r_e(t, u,θ). It should be noted that this firing rate depends on the parameters θ as well as stimulus u. 

The function in (<ref>) is derived from Local Bernoulli approximation which will also be  utilized in the neural spike data generator mechanism. The discussion can be found in <cit.> and the simulation approach is beriefly presented in Section <ref>.



  §.§.§ Maximum Likelihood Methods and Parameter Estimation
  
The maximum-likelihood estimation of the network parameters θ in Table <ref> is based on the likelihood function given by (<ref>), 
which takes the individual spike timings into account. 

We know from estimation theory that maximum likelihood estimation is asymptotically efficient, i.e., reaching the Cramér-Rao bound when the data sample size becomes larger. 

In order to improve the statistical content of the obtained data, one will need to collect multiple data elicited from multiple stimuli. 

Suppose that we have, N_it independent stimuli that yielded N_it different spike trains and suppose also that  m^th stimulus (m=1,…, N_it) elicits a spike train with a total of K_m spikes in the time window [0,T], and the spike timings are given by
S_m=(t_1^(m),t_2^(m),…,t_K_m^(m)). Using (<ref>) one can evaluate the likelihood for the spike train S_m as:


    p(S_m|θ)
    =exp(-∫_0^Tr_e^(m)(t)dt)∏_k=1^K_m r_e^(m)(t_k^(m))
  

In the above, r_e^(m) is the firing rate response to the m^th stimulus. This function depends implicitly on the model parameters θ and on the stimulus. The left-hand side of (<ref>) emphasizes the dependence on network parameters θ, which is convenient for parameter estimation. If the responses to each individual stimulus m are independent (this will be a reasonable assumption especially when the duration between consecutive stimuli are sufficiently large). Under all those, the joint likelihood function for N_it spike trains can be written as shown below:


    L(S_1,S_2,…,S_M|θ)
    =∏_m=1^Mp(S_m |θ)


To ease the computation, one is recommended to evaluate the natural logarithm of (<ref>) and derive the log-likelihood function as:


    l (S_1,S_2,…,S_M|θ)
    =-∑_m=1^M∫_0^Tr_e^(m)(t)dt 
    +∑_m=1^M∑_k=1^K_mlnr_e^(m)(t_k^(m))


The maximization of the above with respect to the parameter vector θ will yield the maximum likelihood estimate. That is:


    θ̂_ML=max_θ[ l  (S_1,S_2,…,S_M|θ)]
 

The above operation can be performed by MATLAB optimization routines such as . A zero lower bound should be assigned for each parameter in order to ensure that positive valued results are obtained. 



 §.§ Data Generating Model
 

In this research, we will generate the neural spiking data by simulating an inhomogeneous Poisson process of which event rate (firing rate) is computed by simulating a different neuron model with predetermined parameters. In this work, we will achieve this goal by simulating the model in <cit.> with its true parameters. In this section we will present a summary for the sake of completeness. 



  §.§.§ Data Generator CTRNN


The data generator neural network is also a continuous time recurrent neural network in the following form:

    v̇_e   =-β_e [  v_e+w_eeg_e(v_e)-w_eig_i(v_i)+c_eu ]  
    v̇_i   =-β_i [  v_i+w_ieg_e(v_e)-w_iig_i(v_i)+c_iu ]
  

where the v_e and v_i represent dynamics of the excitatory and inhibitory unit respectively. They can be thought like a membrane potential. The critical difference from (<ref>) is related to the sigmoidal gain functions g_e(v_e) and g_i(v_i). These have a different form than (<ref>) as shown below:


    g_j(v_j)=Γ_j/1+exp( -a_j(v_j-h_j))

with j being either e for excitatory units or i for inhibitory units. The firing rate from excitatory unit r_e is evaluated as:


    r_e=g_e(v_e)=Γ_e/1+exp(-a_e(v_e-h_e))


The model in (<ref>) and (<ref>) are simulated with the true values of parameters from <cit.> repeated in Tables <ref> and <ref> for convenience. 





The firing rate r_e will be solved from the equation r_e=g_e(v_e). The approach for generation of the spikes are presented in the next section.



  §.§.§ Spike Generation for Data Collection


In order to have a set of neural spiking data that is to be a representative of a collected data from an in vivo experiment we will need to have a reliable algorithm for spike generation. 

In Section <ref>, we stated that sensory neural spiking largely obeys an Inhomogeneous Poisson Process. So one can implement a Poisson process simulation algorithm with the firing rate r_e in (<ref>) being the Poisson event rate. 

There are numerous methodologies to generate the Poisson events given the event rate r_e(t). These ranging from discrete simulation <cit.> to thinning <cit.>.

Discrete simulation may be beneficial when one solves the dynamical models by fixed step solvers such as Euler Integration or Runge-Kutta methods. The only disadvantage of this approach is that, it confines the spikes into discrete time bins. However, if one has a sufficiently small discrete time bin such as Δ t=1 ms, the statistical distribution of the spikes should approach to that of an Inhomogeneous Poisson Process <cit.>.  
Discrete simulation of neural spiking can be summarized as shown below <cit.>:


	
  * Given the firing rate of any neuron as r(t)
	
  * Find the probability of firing at time t_i by evaluating p_i=r(t_i)Δ t where Δ t is the integration interval. It should be as small as 1 ms. 
	
  * Compute a random variable by drawing a sample from a distribution which is uniform between 0 and 1. Define this as x_rand=U[0,1] where U stands for uniform distribution. 
	
  * If p_i>x_rand fire a spike at t=t_i, else do nothing.
	
  * Collect spikes as S=[t_1,…,t_N_s] where N_s will be the number of spikes obtained at a single run of simulation. 




 §.§ Stimulus Model


The stimulus in this research will be of Phased Cosine Fourier Series Type as shown below (same as that of <cit.>):

    I=∑_n=1^N_UA_n cos(ω_nt+ϕ_n)

where A_n is the amplitude, ω_n=2π f_0 n is the frequency of the n^th Fourier component in rad/ and ϕ_n is the phase of the component. Here the amplitude A_n and the base frequency f_0 (in Hz) are fixed but the phase ϕ_n will be a randomly chosen from a uniform distribution between [-π, π] radians. The amplitude parameter A_n is fixed for all mode n as A_n=A_max.

One can see a typical variation of stimulus in (<ref>) for different component sizes N_U=5… 50, A_max=100 and f_0=3.3333 Hz in Figure <ref>.



The typical firing rate response of the data generator model in (<ref>), (<ref>) and (<ref>) with the parameters in Table <ref> and <ref> to the stimuli in Figure <ref> can be seen in Figure <ref>. In Figure <ref> one can see the neural spiking responses of the excitatory unit against the same stimuli. 







§ APPLICATION



 §.§ Step-by-step description of the Problem and Simulation

        
In this section we will present a step-by-step introduction of the example problem. One can refer to the following on this manner:


	
  * A single run of simulation will last for T_f=3 seconds. 
	
  * The data generator neuron model in (<ref>) and (<ref>) will be simulated at the parametric values given in Table <ref> and Table <ref> and firing rate of the excitatory unit is stored as r_e^m(t) where m is the current number of simulation. 
	
  * Firing rate data r_e^m(t) is used to generate neural spikes S_m in the m^th run using the methodology defined in Section <ref>. This data will be used to compute the likelihood. The number of spikes will be K_m at the m^th run.  
	
  *  Repeat the simulation N_it times to obtain enough number of spikes. 
	
  * The spiking data needed by (<ref>) will be obtained at the <ref>th step. However, the firing rate component of (<ref>) should be computed at the current iteration of the optimization. 
	
  * Run an optimization algorithm on the joint likelihood function (<ref>) of which objective computes the firing rate at the current iterated value of the parameters but the spikes from Step <ref>. This procedure will estimate the parameters of the main model in (<ref>) and (<ref>) using maximum likelihood estimation as defined in (<ref>).       




 §.§ Optimization Algorithm


Theoretically, any optimization algorithm ranging from gradient descent to derivative free pattern search, simulated annealing or genetic algorithms can be utilized to achieve the computational goals of this research. In the literature and computational software several different implementations are provided to perform the mentioned operations. In this research, we will implement all computations related to optimization and other general operations by MATLAB. Regardless of the type of algorithm, all of the methods converge to a local optimum and thus will require an initial parametric guess. When one has a convex problem all of the initial guesses are expected to converge to the same solution. In the case of problems involving dynamical models, this may or may not be the case. In any case, the main criteria on the choice of the algorithms is the speed of convergence. This is critical regardless of the available computational facilities as we have a huge data to process. The computations are performed using a Hyper Computing Facility (HPC).  Some initial evaluations suggested that local optimizer routines provided by MATLAB's  should be preferred concerning speed and computational resource considerations. MATLAB also provides derivative free optimization techniques such as  (pattern search technique),  (simulated annealing technique) or  (genetic algorithm technique) but these may require longer durations to obtain a solution. The  algorithm needed gradient information but it can be provided by itself through finite difference approximations In addition,  and all other optimizers mentioned here are local optimizers and thus multiple initial guesses will be required. When the evaluations are performed on an HPC system, each optimization with a different initial guess can be assigned to a specific processor core. So, a faster parallel computing is achieved and result can be obtained from a single pass. MATLAB Parallel Computing Toolbox will help on this manner. The initial guesses are generated randomly from a uniform distribution. 



 §.§ Simulation Scenarios


The data associated with a nominal simulation scenario is presented in Table <ref>. 

      

It will be convenient to analyze the influences of change in stimulus component size N_U, amplitude level A_max and number of samples N_it we will need to repeat the problem for a set of different values of those parameters. The details can be seen in Table <ref>.



The initial conditions of the excitatory and inhibitory units are not accurately known thus we will assume them as zero i.e. x_e(0)=0 and x_i(0)=0. We will repeat the simulation 10 times for each case, so that we will have statistically sufficient number of results for analyses. 



§ RESULTS


We will present the numerical results of the maximum likelihood estimation of the parameters of our neuron model in (<ref>) in Table <ref> by maximizing the joint likelihood function in (<ref>). The optimization is performed using the gradient based interior-point method provided by MATLAB's  algorithm. All the cases in Table <ref> are examined under the conditions in Table <ref>. The results will be presented as:


	
  * The mean estimated values of θ against varying sample size N_it, amplitude level A_max, stimulus component size N_U and base frequency f_0.  
	
  * The standard deviations of the estimated parameters against varying sample size N_it, amplitude level A_max, stimulus component size N_U and base frequency f_0.
	
  * The variations of the mean square error between excitatory firing rate r_e obtained from data generator model and estimated model. This is required to verify the estimation's accuracy as we do not have any true parameters for the estimated model.     




 §.§ Mean Estimated Values


In this section, we will present the estimated mean values of all parameters in Table <ref> against varying scenario parameters in Table <ref> and <ref>. The results are presented in forms of tables starting from Table <ref>. Explanations are given in the table definitions. 

 

 







 §.§ Change in the Standard Deviation of the Estimates


One can see the variation of the standard deviation of the estimates against the sample size N_it, stimulus component size N_U, maximum amplitude parameter A_max and stimulus base frequency f_0 in Figures <ref>, <ref>, <ref>  and <ref> respectively. 

 









 §.§ Comparing the Outputs of Original and Estimated Models


In order to validate the success of our work, we will need to check whether the firing rate outputs of two models are tracking each other. This can be performed by applying a stimulus profile to both original or data generator model in (<ref>) and our estimated model in (<ref>) and compare their firing rate outputs. It will be better to test the models using at least two different stimuli. In this research, we will do this by driving our models by a square wave stimulus such as the one shown in Figure <ref> and a Fourier series one such as the one shown in Figure <ref>. 

 

 

The response obtained from (<ref>) with the estimated parameters corresponding to the case N_it=400, N_U=5, A_max=100 and f_0=3.3333 (the last case in Table <ref>) to the stimuli in Figures <ref> and <ref> can be seen in Figures <ref> and <ref> respectively. 

 

 



§ DISCUSSION




 §.§ Summary


In this research, we considered a theoretical study of model fitting to a noisy stimulus/response data obtained from a generic sensory neural network. In general the purpose is similar to that that of <cit.>. The main differences can be summarized as shown below:


	
  * We desire a more generic model that relates stimulus to firing rate response of the neural system under consideration. Such a model has smaller number of parameters. 
	
  * In <cit.>, it is found that inclusion of sigmoidal gain function parameters to the estimation problem leads to a degradation in estimation performance. For example, increasing the sample size N_it, did not improve the variance of estimates.  


Thus, we preferred a model which does not require the estimation of the sigmoidal gain function parameters. That will decrease the level of parameter confounding so that we expect a reduced level of estimation errors. In addition data is taken from a data generator model which is different in form (the model used in <cit.>). This means that we do not have any true parameters and thus we will need to compare the firing rate responses of the data generator (original model in other words) and the estimated model. Upon completion of the study we noticed the results presented in the next section. 



 §.§ Discussion of the Results


In Tables <ref>, <ref>, <ref> and <ref> one can see the mean estimated values of the parameters of our model in (<ref>) and (<ref>). The values are remaining in a narrow range against changing simulation parameters sample size N_it, stimulus component count N_U, maximum amplitude A_max and base frequency f_0 parameters. One caveat appears when A_max=25, the parameter w_ei and w_ie deviates from their values correspoding to the amplitudes larger than A_max=25. This result verifies that, the nominal selection of A_max=100, N_U=5 and f_0=10/3 is a good compromise. 

The results associated with the variation of standard deviations against increasing simulation parameters are available in Figures <ref>, <ref>, <ref> and <ref>. 

These results are as expected. The standard deviations of estimates of each parameter decreases with increasing sample size N_it. This occurrence is similar to that of <cit.>. 

There isn't any definite pattern associated with the changing stimulus component size N_U. Thus, there is no need to choose large N_U values that brings unnecessary computational complexity.    

Concerning the stimulus amplitude A_max, one can say that lower levels should be avoided. Figure <ref> suggests that, choosing a value like A_max=25 will negatively affect the estimation performance. Same graphical results suggest that a selection of 100 ≤ A_max≤ 200 should be a good compromise. 

The changes in base frequency is also effective on the estimate variance. Figure <ref> suggests that mid-range frequencies (7/3≤ f_0 ≤10/3 Hz) should be preferred. 

Since we are estimating a model of which true parameters are not known, we needed to develop a verification mechanism. This can be done by comparing the firing rate responses of the original or data generator model and the estimated model. In this research, we did that by stimulating both models by a square wave (Figure <ref>) and a Fourier series (Figure <ref>) stimuli. The corresponding responses were given in Figure <ref> against square wave and Figure <ref> against Fourier series stimuli respectively. One can say that under a smooth stimulus like Fourier series both original (data generator) and estimated models generate almost same firing rate response (Figure <ref>). When a discontinuous stimulus like Figure <ref> is applied, the responses elicited by both models have a small deviation at the discontinuities Figure <ref>. 

Based on all those findings one can state that, the estimation procedure is completed successfully.








§ ACKNOWLEDGEMENT

The computational facilities needed by this research are provided by TR-GRID/TRUBA framework (www.truba.gov.tr) operated by National Academic Center of Computing (ULAKBIM) of Turkish Scientific and Technological Research Institution (TUBITAK).



§ REFERENCES


99
	
	hodgkin1952quantitative
	Hodgkin A L and Huxley A F 1952 The Journal of physiology 117 500
	
	fitzhugh1961impulses
	FitzHugh R 1961 Biophysical journal 1 445–466
	
	hindmarsh1984model
	Hindmarsh J L and Rose R 1984 Proc. R. Soc. Lond. B 221 87–102
	
	morris1981voltage
	Morris C and Lecar H 1981 Biophysical journal 35 193–213
	
	ermentrout2001effects
	Ermentrout B, Pascal M and Gutkin B 2001 Neural Computation 13
	1285–1310
	
	booth1997compartmental
	Booth V, Rinzel J and Kiehn O 1997 Journal of Neurophysiology 78
	3371–3385
	
	rust2005spatiotemporal
	Rust N C, Schwartz O, Movshon J A and Simoncelli E P 2005 Neuron 46 945–956
	
	hosoya2005dynamic
	Hosoya T, Baccus S A and Meister M 2005 Nature 436 71
	
	mante2005independence
	Mante V, Frazor R A, Bonin V, Geisler W S and Carandini M 2005 Nature
		neuroscience 8 1690
	
	borst1999information
	Borst A and Theunissen F E 1999 Nature neuroscience 2 947
	
	barlow1959possible
	Barlow H 1959 Possible principles underlying the transformation of sensory
	messages. in. w. rosenblith (ed.) sensory communication (pp. 217-234)
	
	fairhall2001efficiency
	Fairhall A L, Lewen G D, Bialek W and van Steveninck R R d R 2001 Nature 412 787
	
	hassenstein1956systemtheoretische
	Hassenstein B and Reichardt W 1956 Zeitschrift für Naturforschung B
	11 513–524
	
	haykin2004comprehensive
	Haykin S and Network N 2004 Neural networks 2 41
	
	dimattina2011active
	DiMattina C and Zhang K 2011 Neural computation 23 2242–2288
	
	dimattina2013adaptive
	DiMattina C and Zhang K 2013 Frontiers in neural circuits 7
	
	funahashi1993approximation
	Funahashi K i and Nakamura Y 1993 Neural networks 6 801–806
	
	ledoux2011dynamics
	Ledoux E and Brunel N 2011 Frontiers in computational neuroscience 5
	
	doruk2018
	Doruk R O and Zhang K 2018 Journal of Biological Physics ISSN 1573-0689
	<https://doi.org/10.1007/s10867-018-9501-z>
	
	schafer2006recurrent
	Schäfer A M and Zimmermann H G 2006 Recurrent neural networks are universal
	approximators International Conference on Artificial Neural Networks
	(Springer) pp 632–640
	
	phdtamw
	Tam W 2012 Adaptive modeling of marmoset inferior colliculus neurons in
		vivo Ph.D. thesis Department of Biomedical Engineering, Johns Hopkins
	University, Baltimore, MD, USA
	
	hancock1997modeling
	Hancock K E, Davis K A and Voigt H F 1997 Biological cybernetics 76 419–428
	
	hancock1999wideband
	Hancock K E and Voigt H F 1999 Annals of biomedical engineering 27 73–87
	
	de2008linking
	de la Rocha J, Marchetti C, Schiff M and Reyes A D 2008 The Journal of
		Neuroscience 28 9151–9163
	
	truccolo2005point
	Truccolo W, Eden U T, Fellows M R, Donoghue J P and Brown E N 2005 Journal
		of neurophysiology 93 1074–1089
	
	0295-5075-22-2-013
	Rieke F, Warland D and Bialek W 1993 EPL (Europhysics Letters) 22
	151 <http://stacks.iop.org/0295-5075/22/i=2/a=013>
	
	herz2006modeling
	Herz A V, Gollisch T, Machens C K and Jaeger D 2006 Science 314
	80–85
	
	shadlen1994noise
	Shadlen M N and Newsome W T 1994 Current opinion in neurobiology 4 569–579
	
	uteden2008pointprocesses
	Eden U T 2008 Neural Signal Processing: Quantitative Analysis of Neural
		Activity  45–51
	
	10.7287/peerj.preprints.27015v1
	Doruk O and Zhang K 2018 PeerJ Preprints 6 e27015v1 ISSN
	2167-9843 <https://doi.org/10.7287/peerj.preprints.27015v1>
	
	myung2003tutorial
	Myung I J 2003 Journal of Mathematical Psychology 47 90–100
	
	beer1995dynamics
	Beer R D 1995 Adaptive Behavior 3 469–509
	
	lewis1979simulation
	Lewis P A and Shedler G S 1979 Naval Research Logistics Quarterly 26 403–413
	
	brown2002time
	Brown E N, Barbieri R, Ventura V, Kass R E and Frank L M 2002 Neural
		computation 14 325–346
	








