Stochastic neural dynamics mediate between the underlying biophysical and physiological properties of a neural system and its computational and cognitive properties ( e . g . 1–4 ) ., Hence , from a computational perspective , we are often interested in recovering the neural network dynamics of a given brain region or neural system from experimental measurements ., Yet , experimentally , we commonly have access only to noisy recordings from a relatively small proportion of neurons ( compared to the size of the brain area of interest ) , or to lumped surface signals like local field potentials or the EEG ., Inferring from these the computationally relevant dynamics is therefore not trivial , especially since both the recorded signals ( e . g . , spike sorting errors; 5 ) as well as the neural system dynamics itself ( e . g . , stochastic synaptic release; 6 ) come with a good deal of noise ., The stochastic nature of neural dynamics has , in fact , been deemed crucial for perceptual inference and decision making 7–9 , and potentially helps to avoid local minima in task learning or problem solving 10 ., Speaking in statistical terms , model-free techniques which combine delay embedding methods with nonlinear basis expansions and kernel techniques have been one approach to the problem 11; 12 ., These techniques provide informative lower-dimensional visualizations of population trajectories and ( local ) approximations to the neural flow field , but they may highlight only certain , salient aspects of the dynamics ( but see 13 ) and , in any case , do not directly return distribution generating equations or underlying computations ., Alternatively , state space models , a statistical framework particularly popular in engineering and ecology ( e . g . 14 ) , have been adapted to extract lower-dimensional , probabilistic neural trajectory flows from higher-dimensional recordings 15–25 ., State space models link a process model of the unobserved ( latent ) underlying dynamics to the experimentally observed time series via observation equations , and differentiate between stochasticity in the process and observation noise ( e . g . 26 ) ., So far , with few exceptions ( e . g . 23; 27 ) , these models assumed linear latent dynamics , however ., Although this may often be sufficient to yield lower-dimensional smoothed trajectories , it implies that the recovered dynamical model may be less apt for capturing highly nonlinear dynamical phenomena in the observations , and will by itself not be powerful enough to reproduce a range of important dynamical and computational processes in the nervous system , among them multi-stability which has been proposed to underlie neural activity during working memory 28–32 , limit cycles ( stable oscillations ) , or chaos ( e . g . 33 ) ., Here we derive a new state space algorithm based on piecewise-linear ( PL ) recurrent neural networks ( RNN ) ., It has been shown that RNNs with nonlinear activation functions can , in principle , approximate any dynamical systems trajectory or , in fact , dynamical system itself ( given some general conditions; 34–36 ) ., Thus , in theory , they are powerful enough to recover whatever dynamical system is underlying the experimentally observed time series ., Piecewise linear activation functions , in particular , are by now the most popular choice in deep learning algorithms 37–39 , and considerably simplify some of the derivations within the state space framework ( as shown later ) ., They may also be more apt for producing working memory-type activity with longer delays if for some units the transfer function happens to coincide with the bisectrix ( cf . 40 ) , and ease the analysis of fixed points and stability ., We then apply this newly derived algorithm to multiple single-unit recordings from the rat prefrontal cortex obtained during a classical delayed alternation working memory task 41 ., This article considers simple discrete-time piecewise-linear ( PL ) recurrent neural networks ( RNN ) of the form, zt=Azt−1+Wmax{0 , zt−1−θ}+Cst+εt , εt∼N ( 0 , Σ ) ,, ( 1 ), where zt = ( z1t…zMt ) T is the ( M×1 ) -dimensional ( latent ) neural state vector at time t = 1…T , A = diag ( a11…aMM ) is an M×M diagonal matrix of auto-regression weights , W = ( 0 w12…w1M , w21 0 w23…w2M , w31 w32 0 w34…w3M , … ) is an M×M off-diagonal matrix of connection weights , θ = ( θ1…θM ) T is a set of ( constant ) activation thresholds , st is a sequence of ( known ) external K-dimensional inputs , weighted by ( M×K ) matrix C , and εt denotes a Gaussian white noise process with diagonal covariance matrix Σ=diag ( σ112…σMM2 ) ., The max-operator is assumed to work element-wise ., In physiological terms , latent variables zmt are often interpreted as a membrane potential ( or current ) which gives rise to spiking activity as soon as firing threshold θm is exceeded ( e . g . 42 , 43 ) ., According to this interpretation , the diagonal elements in A may be seen as the neurons’ individual membrane time constants , while the off-diagonal elements in W represent the between-neuron synaptic connections which multiply with the presynaptic firing rates ., In statistical terms , ( 1 ) has the form of an auto-regressive model with a nonlinear basis expansion in variables zmt ( e . g . 44;45 ) , which retains linearity in parameters W for ease of estimation ., Restricting model parameters , e . g . Σ , to be of diagonal form , is common in such models to avoid over-specification and help identifiabiliy ( e . g . 26; 46; see also further below ) ., For instance , including a diagonal in W would be partly redundant to parameters A ( strictly so in a pure linear model ) ., For similar reasons , and for ease of presentation , in the following we will focus on a model for which K = M and C = I ( i . e . , no separate scaling of the inputs ) , although the full model as stated above , Eq 1 , was implemented as well ( and code for it is provided; of course , the case K>M could always be accommodated by pre-multiplying st by some predefined matrix C , obtained e . g . by PCA on the input space ) ., While different model formulations are around in the computational neuroscience and machine learning literature , often they may be related by a simple transformation of variables ( see 47 ) and , as long as the model is powerful enough to express the whole spectrum of basic dynamical phenomena , details of model specification may also not be overly crucial for the present purposes ., A particular advantage of the PLRNN model is that all its fixed points can be obtained easily analytically by solving ( in the absence of external input ) the 2M linear equations, z*= ( A+WΩ−I ) −1WΩθ ,, ( 2 ), where Ω is to denote the set of indices of units for which we assume zm ≤ θm , and WΩ the respective connectivity matrix in which all columns from W corresponding to units in Ω are set to 0 ., Obviously , to make z* a true fixed point of ( 1 ) , the solution to ( 2 ) has to be consistent with the defined set Ω , that is z*m ≤ θm has to hold for all m ∈ Ω and z*m > θm for all m ∉ Ω ., For networks of moderate size ( say M<30 ) it is thus computationally feasible to explicitly check for all fixed points and their stability ., For estimation from experimental data , latent state model ( 1 ) is then connected to some N-dimensional observed vector time series X = {xt} via a simple linear-Gaussian model ,, xt=Bϕ ( zt ) +ηt , ηt∼N ( 0 , Γ ) ,, ( 3 ), where ϕ ( zt ) ≔ max{0 , zt−θ} , {ηt} is the ( white Gaussian ) observation noise series with diagonal covariance matrix Γ=diag ( γ112…γNN2 ) , and B an N×M matrix of regression weights ., Thus , the idea is that only the PL-transformed activation ϕ ( zt ) reaches the ‘observation surface’ as , e . g . , with spiking activity when the underlying membrane dynamics itself is not visible ., We further assume for the initial state ,, z1∼N ( μ0+s1 , Σ ) ,, ( 4 ), which has , for simplicity , the same covariance matrix as the process noise in general ( reducing the number of to be estimated parameters ) ., In the case of multiple , temporally separated trials , we allow each one to have its own individual initial condition μk , k = 1…K ., The general goal here is to determine both the model’s unknown parameters Ξ = {μ0 , A , W , Σ , B , Γ} ( assuming fixed thresholds θ for now ) as well as the unobserved , latent state path Z ≔ {zt} ( and its second-order moments ) from the experimentally observed time series {xt} ., These could be , for instance , properly transformed multivariate spike time series or neuroimaging data ., This is accomplished here by the Expectation-Maximization ( EM ) algorithm which iterates state ( E ) and parameter ( M ) estimation steps and is developed in detail for model ( 1 ) and ( 3 ) in the Methods ., In the following I will first discuss state and parameter estimation separately for the PLRNN , before describing results from the full EM algorithm in subsequent sections ., This will be done along two toy problems , a higher-order nonlinear oscillation ( stable limit cycle ) , and a simple working memory paradigm in which one of two discrete stimuli had to be retained across a temporal interval ., Finally , the application of the validated PLRNN EM algorithm will be demonstrated on multiple single-unit recordings obtained from rats on a standard working memory task ( delayed alternation; data from 41 , kindly provided by Dr . James Hyman , University of Nevada , Las Vegas ) ., The latent state distribution , as explained in Methods , is a high-dimensional ( piecewise ) Gaussian mixture with the number of components growing as 2T×M with sequence length T and number of latent states M . Here a semi-analytical , approximate approach was developed that treats state estimation as a combinatorial problem by first searching for the mode of the full distribution ( cf ., 16; 48; in contrast , e . g . , to a recursive filtering-smoothing scheme that makes local ( linear-Gaussian ) approximations , e . g . 15; 26 ) ., This approach amounts to solving a high ( 2M×T ) -dimensional piecewise linear problem ( due to the piecewise quadratic , in the states Z , log-likelihood Eqs 6 and 7 ) ., Here this was accomplished by alternating between ( 1 ) solving the linear set of equations implied by a given set of linear constraints Ω ≔ { ( m , t ) |zmt ≤ θm} ( cf . Eq 7 in Methods ) and ( 2 ) flipping the sign of the constraints violated by the current solution z* ( Ω ) to the linear equations , thus following a path through the ( M×T ) -dimensional binary space of linear constraints using Newton-type iterations ( similar as in 49 , see Methods; note that here the ‘constraints’ are not fixed as in quadratic programming problems ) ., Given the mode and state covariance matrix ( evaluated at the mode from the negative inverse Hessian ) , all other expectations needed for the EM algorithm were then derived analytically , with one exception that was approximated ( see Methods for full details ) ., The toy problems introduced above were used to assess the quality of these approximations ., For the first toy problem , an order-15 limit cycle was produced with a PLRNN consisting of three recurrently coupled units , inputs to units #1 and #2 , and parameter settings as indicated in Fig 1 and provided Matlab file ‘PLRNNoscParam’ ., The limit cycle was repeated for 50 full cycles ( giving 750 data points ) and corrupted by process noise ( cf . Fig 1 ) ., These noisy states ( arranged in a ( 3 x 750 ) matrix Z ) were then transformed into a ( 3 x 750 ) output matrix X , to which observation noise was added , through a randomly drawn ( 3 x 3 ) regression weight matrix B . State estimation was started from a random initial condition ., True ( but noise-corrupted ) and estimated states for this particular problem are illustrated in Fig 1A , indicating a tight fit ( although some fraction of the linear constraints were still violated , ≈0 . 27% in the present example and <2 . 3% in the working memory example below; see Methods on this issue ) ., To examine more systematically the quality of the approximate-analytical estimates of the first and second order moments of the joint distribution across states z and their piecewise linear transformations ϕ ( z ) , samples from p ( Z|X ) were simulated using bootstrap particle filtering ( see Methods ) ., Although these simulated samples are based only on the filtering ( not the smoothing ) steps ( and ( re- ) sampling schemes may have issues of their own; e . g . 26 , analytical and sampling estimates were in tight agreement , correlating almost to 1 for this example , as shown in Fig 2 ., Fig 3A illustrates the setup of the ‘two-cue working memory task’ , chosen for later comparability with the experimental setup ., A 5-unit PLRNN was first trained by conventional gradient descent ( ‘real-time recurrent learning’ ( RTRL ) , see 50; 51 ) to produce a series of six 1’s on unit #3 and six 0’s on unit #4 five time steps after an input ( of 1 ) occurred on unit #1 , and the reverse pattern ( six 0’s on unit #3 and six 1’s on unit #4 ) five time steps after an input occurred on unit #2 ., A stable PLRNN with a reasonable solution to this problem was then chosen for further testing the present algorithm ( cf . Fig 3C ) ., ( While the RTRL approach was chosen to derive a working memory circuit with reasonably ‘realistic’ characteristics like a wider distribution of weights , it is noted that a multi-stable network is relatively straightforward to construct explicitly given the analytical accessibility of fixed points ( see Methods ) ; for instance , choosing θ = ( 0 . 5 0 . 5 0 . 5 0 . 5 2 ) , A = ( 0 . 9 0 . 9 0 . 9 0 . 9 0 . 5 ) , and W = ( 0 ω − ω − ω − ω , ω 0 − ω − ω – ω , − ω − ω 0 ω – ω , − ω − ω ω 0 − ω , 11110 ) with ω = 0 . 2 , yields a tri-stable system ., ) Like for the limit cycle problem before , the number of observations was taken to be equal to the number of latent states , and process and observation noise were added ( see Fig 4 and Matlab file ‘PLRNNwmParam’ for specification of parameters ) ., The system was simulated for 20 repetitions of each trial type ( i . e . , cue-1 or cue-2 presentations ) with different noise realizations and each ‘trial’ started from its own initial condition μk ( see Methods ) , resulting in a total series length of T = 20×2×20 = 800 ( although , importantly , in this case the time series consisted of distinct , temporally segregated trials , instead of one continuous series , and was treated as such an ensemble of series by the algorithm ) ., As before , state estimation started from random initial conditions and was provided with the correct parameters , as well as with the observation matrix X . While Fig 3B illustrates the correlation between true ( i . e . , simulated ) and estimated states across all trials and units , Fig 3C shows true and estimated states for a representative cue-1 ( left ) and cue-2 ( right ) trial , respectively ., Again , our procedure for obtaining ( or approximating ) the maximum a-posteriori ( MAP ) estimate of the state distribution appears to work quite well ( in general , only locally optimal or approximate solutions may be achieved , however , and the algorithm may have to be repeated with different state initializations; see Methods ) ., Given the true states , how well would the algorithm retrieve the parameters of the PLRNN ?, To assess this , the actual model states ( which generated the observations X ) from simulation runs of the oscillation and the working memory task described above were provided as initialization for the E-step ., Based on these , the algorithm first estimated the state covariances for z and ϕ ( z ) ( see above ) , and then the parameters in a second step ( i . e . , the M-step ) ., Note that the parameters can all be computed analytically given the state distribution ( see Methods ) , and , provided the state covariance matrices ( summed across time ) as required in Eq 17A , 17D and 17F are non-singular , have a unique solution ., Hence , in this case , any misalignment with the true model parameters can only come from one of two sources:, i ) estimation was based on one finite-length noisy realization of the PLRNN process ,, ii ) all second order moments of the state distribution were still estimated based on the true state vectors ., However , as can be appreciated from Fig 1B ( oscillation ) and Fig 4 ( working memory ) , for the two ( relatively low-noise ) example scenarios studied here , all parameter estimates still agreed tightly with those describing the true underlying model ., In the more general case where both the states and the parameters are unknown and only the observations are given , note that the model as stated in Eqs 1 & 3 is over-specified as , for instance , at the level of the observations , additional variance placed into Σ may be compensated for by adjusting Γ accordingly , and by rescaling W and , within limits , A ( cf . 52; 53 ) ., In the following we therefore always arbitrarily fixed Σ ( to some scalar; see Methods ) , as common in many latent variable models ( like factor analysis ) , including state space models ( e . g . 27; 46 ) ., It may be worth noting here that the relative size of Σ vs . Γ determines how much weight is put on temporal consistency among states ( “Σ<Γ” ) vs . fitting of the observations ( “Σ>Γ” ) within the likelihood , Eq 5 ., The observations above confirm that our algorithm finds satisfactory approximations to the underlying state path and state covariances when started with the right parameters , and—vice versa—identifies the correct parameters when provided with the true states ., Indeed , the M-step , since it is exact , can only increase the expected log-likelihood Eq 5 with the present state expectancies fixed ., However , due to the systems piecewise-defined discrete nature , modifying the parameters may lead to a new set of constraint violations , that is may throw the system into a completely different linear subspace which may imply a decrease in the likelihood in the E-step ., It is thus not guaranteed that a straightforward EM algorithm converges ( cf . 54; 55 ) , or that the likelihood would even monotonically increase with each EM iteration ., To examine this issue , full EM estimation of the WM model ( as specified in Fig 4 , using N = 20 outputs in this case ) was run 240 times , starting from different random , uniformly distributed initializations for the parameters ., Fig 5B ( Δt = 0 ) gives , for the five highest likelihood solutions across all 240 runs ( Fig 5A ) , the mean squared error ( MSE ) avg ( xit−x^it ) 2 between actual neural observations xit and model predictions x^it , which is close to 0 ( and , correspondingly , correlations between predicted and actual observations were close to 1 ) ., With respect to the inferred states , note that estimated and true model states may not be in the same order , as any permutation of the latent state indices together with the respective columns of observation matrix B will be equally consistent with the data X ( see also 27 ) ., For the WM model examined here , however , partial order information is implicitly provided to the EM algorithm through the definition of unit-specific inputs sit ., For the present example , true and estimated states for the highest likelihood solution were nicely linearly correlated for all 5 latent variables ( Fig 6 ) , but some of the regression slopes significantly differed from 1 , indicating a degree of freedom in the scaling of the states ., Note that if the system were strictly linear , the states would be identifiable only up to a linear transformation in general , since any multiplication of the latent states by some matrix V could essentially be reversed at the level of the outputs by back-multiplying B with V-1 ( cf . 27 ) ., Likewise , in the present piecewise linear system , one may expect that there is a class of piecewise-linear transformations of the states which is still compatible with the observed outputs , and hence that the model is only identifiable up to this class of transformations ( a general issue with state space models , of course , not particular to the present one; cf . 53 ) ., However , this might not be a too serious issue , if one is primarily interested in the latent dynamics ( rather than in the exact parameters ) ., Fig 7 illustrates the distribution of initial and final parameter estimates around their true values across all 240 runs ( before and after reordering the estimated latent states based on the rotation that would be required for achieving the optimal mapping onto the true states , as determined through Procrustes analysis ) ., Fig 7 reveals that, a ) the EM algorithm does clearly improve the estimates and, b ) these final estimates seemed to be relatively ‘unbiased’ ( i . e . , with deviations centered around 0 ) ., How do the computational costs of the algorithm grow as the number of latent variables in the model is increased ?, As pointed out in Paninski et al . 16 , exploiting the block-tridiagonal nature of the covariance matrices , the numerical operations within one iteration of the state inference algorithm ( i . e . , solving ∂QΩ* ( Z ) /∂Z=0 , Eq 7 ) can be done in linear , O ( M×T ) , time , just like with the Kalman filter ( due to the model’s Markov properties , full inversion of the Hessian is also not necessary to obtain the relevant moments of the posterior state distribution ) ., This leaves open the question of how many more mode search iterations , i . e . linear equation solving ( Eq 7 ) and constraint-flipping ( vector dΩ ) steps , are required as the number of latent variables ( through either M or T ) increases ., The answer is provided in Fig 8A which is based on the experimental data set discussed below ., Although a full computational complexity analysis is beyond the scope of this paper , at least for these example data ( and similar to what has sometimes been reported for the somewhat related Simplex algorithm; 56 ) , the increase with M appears to be at most linear ., Likewise , the total number of iterations within the full EM procedure , i . e . the number of mode-search steps summed across all EM iterations ( thus reflecting the overall scaling of the full algorithm ) , was about linear ( Fig 8B; in this case , single-constraint instead of complete flipping ( see Methods ) was used which , of course , increases the overall number of iterations but may perform more stably; note that in general the absolute number of iterations will also depend on detailed parameter settings of the algorithm , like the EM convergence criterion and error tolerance ) ., Thus , overall , the present state inference algorithm seems to behave quite favorably , with an at most linear increase in the number of iterations required as the number of latent variables is ramped up ., I next was interested in what kind of structure the present PLRNN approach would retrieve from experimental multiple ( N = 19 ) single-unit recordings obtained while rats were performing a simple and well-examined working memory task , namely spatial delayed alternation 41 ( see Methods ) ., ( Note that in the present context this analysis is mainly meant as an exemplification of the current model approach , not as a detailed examination of the working memory issue itself . ), The delay was always initiated by a nose poke of the animal into a port located on the side opposite from the response levers , and had a minimum length of 10 s ., Spike trains were first transformed into kernel density estimates by convolution with a Gaussian kernel ( see Methods ) , as done previously ( e . g . 12; 57; 58 ) , and binned with 500 ms resolution ., This also renders the observed data more suitable to the Gaussian noise assumptions of the present observation model , Eq 3 ., Models with different numbers of latent states were estimated , with M = 5 or M = 10 chosen for the examples below ., Periods of cue presentation were indicated to the model by setting external inputs sit = 1 to units i = 1 ( left lever ) or i = 2 ( right lever ) for three 500 ms time bins surrounding the event ( and sit = 0 otherwise ) , and the response period was indicated by setting s3t = 1 for 3 consecutive time bins irrespective of the correct response side ( i . e . , non-discriminatively ) ., The EM algorithm was started from a range of different initializations of the parameters ( including thresholds θ ) , and the 5 highest likelihood solutions were considered further for the examples below ., Fig 10A gives the log-likelihoods across EM iterations for these 5 highest-likelihood solutions ( starting from 36 different initializations ) for the M = 5 model ., Interestingly , there were single neurons whose responses were predicted quite well by the estimated model despite large trial-to-trial fluctuations ( Fig 9A , top row ) , while there were others with similar trial-to-trial fluctuations for which the model only captured the general trend ( Fig 9A , bottom row; to put this into context , Fig 9B gives the full distribution of correlations between actual and predicted observations across all 19 neurons ) ., This may potentially indicate that trial-to-trial fluctuations in single neurons could be for very different reasons: For instance , in those cases where strongly varying single unit responses are nevertheless tightly reproduced by the estimated model , a larger proportion of their trial-to-trial fluctuations may have been captured by the latent state dynamics , ultimately rooted in different ( trial-unique ) initializations of the states ( recall that the states are not completely free to vary in accounting for the observations , but are constrained by the model’s temporal consistency requirements ) ., In contrast , when only the average trend is captured , the neuron’s trial-to-trial fluctuations may be more likely to represent true intrinsic ( or measurement ) noise sources that the model’s deterministic part cannot account for ., In practice , such conclusions would have to be examined more carefully to rule out that no other factors in the estimation procedure , like different local maxima , initializations , or over-fitting issues ( see below ) , could account for these differences ., Although this was not further investigated here , this observation nevertheless highlights the potential of ( nonlinear ) state space models to possibly provide new insights also into other long-standing issues in neurophysiology ., Cross-validation is an established means to address over-fitting 45 , although due to the presence of both unknown parameters and unknown states , its application to state space models and its interpretation in this context may be a bit less straightforward ., Here the cross-validation error was first assessed by leaving out each of the 14 experimental trials in turn , estimating model parameters Ξ from the remaining 13 trials , inferring states zt given these parameters on the left-out trial , and computing the squared prediction errors ( xit−x^it ) 2 between actual neural observations xit and model predictions x^it on the left-out trial ., As shown in Fig 10B , this measure steadily ( albeit sub-linearly ) decreases as the number M of latent states in the model is increased ., At first sight , this seems to suggest that with M = 5 or even M = 10 the over-fitting regime is not yet reached ., On the other hand , the latent states are , of course , not completely fixed by the transition equations , but have some freedom to vary as well ( the true effective degrees of freedom for such systems are in fact very hard to determine , cf . 59 ) ., Hence , we also examined the Δt-step-ahead prediction errors , that is , when the transition model were iterated Δt steps forward in time , and x^i , t+Δt=bi⋅ϕ ( z^t+Δt ) estimated from the deterministically predicted states z^t+Δt=HΔt ( Ezt ) ( with HΔt the Δt-times iterated map H ( zt ) = Azt + Wϕ ( zt ) + Cst ) , not from the directly inferred states ( that is , predictions were made on data points which were neither used to estimate parameters nor to infer the current state ) ., These curves are shown for Δt = 1 and Δt = 3 in Fig 10C , and confirm that M = 5 might be a reasonable choice at which over-fitting has not yet ensued ., ( Alternatively , the predictive log-likelihood , logp ( Xtest|Ξ^train ) =log∫p ( Xtest|Z^ ) p ( Z^|Ξ^train ) dZ^ , may be used for model selection ( i . e . , choice of M ) , with p ( Z^|Ξ^train ) either approximated through the E-step algorithm ( with all X-dependent terms removed ) , or bootstrapped by generating Z^-trajectories from the model with parameters Ξ^train ( note that this is different from particle filtering since p ( Z^|Ξ^train ) does not depend on test observations Xtest ) ., This is of course , however , computationally more costly to evaluate than the Δt-step-ahead prediction error ., ) Fig 11 shows trial-averaged latent states for both left- and right-lever trials , illustrated in this case for one of the five highest likelihood solutions ( starting from 100 different initializations ) for the M = 10 model ., Recall that the first 3 PLRNN units received external inputs to indicate left cue ( i = 1 ) , right cue ( i = 2 ) , or response ( i = 3 ) periods , and so , not too surprisingly , reflect these features in their activation ., On the other hand , the cue response is not very prominent in unit i = 1 , indicating that activity in the driven units is not completely dominated by the external regressors either , while unit i = 10 ( not externally driven ) shows a clear left-cue response ., Most importantly , many of the remaining state variables clearly distinguish between the left and right lever options throughout the delay period of the task , in this sense carrying a memory of the cue ( previous response ) within the delay ., Some of the activation profiles appear to systematically climb or decay across the delay period , as reported previously ( e . g . 1; 60 ) , but are a bit harder to read ( at least in the absence of more detailed behavioral information ) , such that one may want to stick with the simpler M = 5 model discussed above ., Either way , for this particular data set , the extracted latent states appear to summarize quite well the most salient computational features of this simple working memory task ., Further insight about the dynamical mechanisms of working memory might be gained by examining the system’s fixed points and their eigenvalue spectrum ., For this purpose , the EM algorithm was started from 400 different initial conditions ( that is , initial parameter estimates and threshold settings θ ) with maximum absolute eigenvalues ( of the corresponding fixed points ) drawn from a relatively uniform distribution within the interval 0 3 ., Although the estimation process rarely returned truly multi-stable solutions ( just 2 . 5% of all cases ) , one frequently discussed candidate mechanism for working memory ( e . g . 29; 32 ) , there was a clear trend for the final maximum absolute eigenvalues to aggregate around 1 ( Fig 12 ) ., For the discrete-time dynamical system ( 1 ) this implies it is close to a bifurcation , with fixed points on the brink of becoming unstable , and will tend to produce ( very ) slow dynamics as the degree of convergence shrinks to zero along the maximum eigenvalue direction ( strictly , a single eigenvalue near 1 does not yet guarantee a slow approach , but makes it very likely , especially in a ( piecewise ) linear system ) ., Indeed , effectively slow dynamics is all that is needed to bridge the delays ( see also 1 ) , while true multi-stability may perhaps even be the physiologically less likely scenario ( e . g . 61; 62 ) ., ( Reducing the bin width from 500 ms to 100 ms appeared to produce solutions with eigenvalues even closer to 1 while retaining stimulus selectivity across the delay , but this observation was not followed up more systematically here ) ., Linear dynamical systems ( LDS ) have frequently and successfully been used to infer smooth neural trajectories from spike train recordings 15; 16; 20; 22 or other measurement modalities 63 ., However , as noted before , they 
The computational and cognitive properties of neural systems are often thought to be implemented in terms of their ( stochastic ) network dynamics ., Hence , recovering the system dynamics from experimentally observed neuronal time series , like multiple single-unit recordings or neuroimaging data , is an important step toward understanding its computations ., Ideally , one would not only seek a ( lower-dimensional ) state space representation of the dynamics , but would wish to have access to its statistical properties and their generative equations for in-depth analysis ., Recurrent neural networks ( RNNs ) are a computationally powerful and dynamically universal formal framework which has been extensively studied from both the computational and the dynamical systems perspective ., Here we develop a semi-analytical maximum-likelihood estimation scheme for piecewise-linear RNNs ( PLRNNs ) within the statistical framework of state space models , which accounts for noise in both the underlying latent dynamics and the observation process ., The Expectation-Maximization algorithm is used to infer the latent state distribution , through a global Laplace approximation , and the PLRNN parameters iteratively ., After validating the procedure on toy examples , and using inference through particle filters for comparison , the approach is applied to multiple single-unit recordings from the rodent anterior cingulate cortex ( ACC ) obtained during performance of a classical working memory task , delayed alternation ., Models estimated from kernel-smoothed spike time data were able to capture the essential computational dynamics underlying task performance , including stimulus-selective delay activity ., The estimated models were rarely multi-stable , however , but rather were tuned to exhibit slow dynamics in the vicinity of a bifurcation point ., In summary , the present work advances a semi-analytical ( thus reasonably fast ) maximum-likelihood estimation framework for PLRNNs that may enable to recover relevant aspects of the nonlinear dynamics underlying observed neuronal time series , and directly link these to computational properties .
Neuronal dynamics mediate between the physiological and anatomical properties of a neural system and the computations it performs , in fact may be seen as the ‘computational language’ of the brain ., It is therefore of great interest to recover from experimentally recorded time series , like multiple single-unit or neuroimaging data , the underlying stochastic network dynamics and , ideally , even equations governing their statistical evolution ., This is not at all a trivial enterprise , however , since neural systems are very high-dimensional , come with considerable levels of intrinsic ( process ) noise , are usually only partially observable , and these observations may be further corrupted by noise from measurement and preprocessing steps ., The present article embeds piecewise-linear recurrent neural networks ( PLRNNs ) within a state space approach , a statistical estimation framework that deals with both process and observation noise ., PLRNNs are computationally and dynamically powerful nonlinear systems ., Their statistically principled estimation from multivariate neuronal time series thus may provide access to some essential features of the neuronal dynamics , like attractor states , generative equations , and their computational implications ., The approach is exemplified on multiple single-unit recordings from the rat prefrontal cortex during working memory .
neural networks, applied mathematics, neuroscience, learning and memory, simulation and modeling, algorithms, cognitive neuroscience, systems science, mathematics, cognition, algebra, memory, research and analysis methods, computer and information sciences, animal cells, dynamical systems, nonlinear dynamics, approximation methods, working memory, cellular neuroscience, cell biology, linear algebra, neurons, biology and life sciences, cellular types, physical sciences, cognitive science, eigenvalues