Variational online learning of neural dynamics
Yuan Zhao Il Memming Park
stat.ML




Bandit Convex Optimization for
 Scalable and Dynamic IoT Management
    
Tianyi Chen and Georgios B. Giannakis
Work in this paper was supported by NSF 1509040, 1508993, and 1711471.

T. Chen and G. B. Giannakis are with the Department of Electrical and Computer Engineering and the Digital Technology Center, University of Minnesota, Minneapolis, MN 55455 USA. Emails: {chen3827, georgios}@umn.edu




    December 30, 2023
=============================================================================================================================================================================================================================================================================================================================================




	New technologies for recording the activity of large neural populations during complex behavior provide exciting opportunities for investigating the neural computations that underlie perception, cognition, and decision-making.
	Nonlinear state space models provide an interpretable signal processing framework by combining an intuitive dynamical system with a probabilistic observation model, which can provide insights into neural dynamics, neural computation, and development of neural prosthetics and treatment through feedback control.
	It brings the challenge of learning both latent neural state and the underlying dynamical system because neither is known for neural systems a priori.
	We developed a flexible online learning framework for latent nonlinear state dynamics and filtered latent states. Using the stochastic gradient variational Bayes approach, our method jointly optimizes the parameters of the nonlinear dynamical system, the observation model, and the black-box recognition model.
	Unlike previous approaches, our framework can incorporate non-trivial distributions of observation noise and has constant time and space complexity. These features make our approach amenable to real-time applications and the potential to automate analysis and experimental design in ways that testably track and modify behavior using stimuli designed to influence learning.




§ INTRODUCTION


Discovering interpretable structure from a streaming high-dimensional time series has many applications in science and engineering.
Since the invention of the celebrated Kalman filter, state space models have been successful in providing a succinct, hence more interpretable, description of the underlying dynamics that explains the observed time series as trajectories in a low-dimensional state space.
Taking a step further, state space models equipped with nonlinear dynamics provide an opportunity to describe the latent “laws” of the system that is generating the seemingly entangled time series <cit.>.
Specifically, we are concerned with the problem of identifying a continuous nonlinear dynamics in the state space (t) ∈ℝ^d that captures the spatiotemporal structure of a noisy observation (t):

	
    = F_θ((t), (t))
    		
    		    (state dynamics)
    (t)       ∼ P((t) | G_θ((t), (t)))     (observation model)


where F and G are continuous functions that may depend on parameter θ, (t) is control input, and P denotes a probability distribution that captures the noise in the observation; e.g., Gaussian distribution for field potentials or Poisson distribution for spike counts.

In practice, the continuous-time state dynamics is more conveniently formulated in discrete time as,

    _t+1    = f_θ(_t, _t) + ϵ_t     (discrete time state dynamics)

where ϵ_t is intended to capture the unobserved (latent) perturbations of the state _t.
Such (spatially) continuous state space models are natural in many applications where the changes are slow and the underlying system follows physical laws and constraints (e.g., object tracking), or where learning the laws are of great interest (e.g. in neuroscience and robotics) <cit.>.
Specifically, in the context of neuroscience, the state vector _t represents the instantaneous state of the neural population, while f captures the time evolution of the population state.
Further interpretation of f can provide understanding as to how neural computation is implemented <cit.>.

If the nonlinear state space model is fully specified, Bayesian inference methods can be employed to estimate the current state <cit.>.
Conventionally, the estimation of latent states using only the past observation is referred to as filtering – inference of the filtering distribution, p(_t |_≤ t).
If both past and future observations are used, then the quantity of interest is usually the smoothing distribution, p(_≤ t|_≤ t)).
We are also interested in predicting the distribution over future states, p(_t:t+s|_≤ t) and observations, p(_t+1:t+s|_≤ t) for s > 0.
However, in many applications, the challenge is in learning the parameters θ of the state space model (a.k.a. the system identification problem).
We aim to provide a method for simultaneously learning both the latent trajectory _t and the latent (nonlinear) dynamical and observational system θ, known as the joint estimation problem <cit.>.

Expectation maximization (EM) based methods have been widely used in practice <cit.>, and more recently variational autoencoder methods <cit.> have been proposed, all of which are designed for offline analysis, and not appropriate for real-time applications.
Recursive stochastic variational inference has been successful in streaming data assuming independent samples <cit.>, however, in the presence of temporal dependence, proposed variational algorithms (e.g. <cit.>) remain theoretical and lack testing.

In this study, we are interested in real-time signal processing and state space control setting <cit.> where online algorithms are needed that can recursively solve the joint estimation problem on streaming observations.
A popular solution to this problem exploits the fact that online state estimators for nonlinear state space models such as extended Kalman filter (EKF) or unscented Kalman filter (UKF) can be used for nonlinear regression formulated as a state space model.







By augmenting the state space with the parameters, one can build an online dual estimator using nonlinear Kalman filtering <cit.>.
However, they involve coarse approximation of Bayesian filtering, involve many hyperparameters, do not take advantage of modern stochastic gradient optimization, and are not easily applicable to arbitrary observation likelihoods.
There are also closely related online version of EM-type algorithms <cit.> that share similar concerns.

In hopes of lifting these concerns, we derive an online black-box variational inference framework, referred to as variational joint filtering (VJF), applicable to a wide range of nonlinear state dynamics (dynamic models) and observation models, that is, the computational demand of the algorithm is constant per time step.
Our approach aims at

	
  * Online adaptive learning: Our target application scenarios are streaming data. This allows the inference during an experiment or as part of a neural prosthetics. If the system changes, the inference will catch up with the altered system parameters.
	
  * Joint estimation: The proposed method is supposed to simultaneously learn the latent states p(_t |_≤ t), state dynamics f(_t, _t) and the observation model G(, ). No offline training is necessary to learn the system parameters.
	
  * Interpretability: Under the framework of state space modeling, rather than interpret the system via model parameters, we employ the language of dynamical systems and capture the characteristics of the system qualitatively via fixed point, limit cycle, strange attractor, bifurcation and so on which are key components of theories of neural dynamics and computation.

We focus on low-dimensional latent dynamics that often underlie many neuroscientific experiments and allow for producing interpretable visualizations of complex collective network dynamics in this study.



§ VARIATIONAL PRINCIPLE FOR ONLINE JOINT ESTIMATION


The crux of recursive Bayesian filtering is updating the posterior over the latent state one step at a time:


    p( _t |_≤ t )     =
    	p( _t |_t )_likelihood p( _t |_<t )_prior at time t 
    	/p(_t |_<t)_marginal likelihood


where the input _t and parameters θ are omitted for brevity.
Unfortunately, the exact calculations of Eq. (<ref>) are not tractable in general, especially for nonlinear dynamic models and/or non-conjugate distributions.
We thus turn to approximate inference and develop a recursive variational Bayesian filter by deriving an evidence lower bound for the marginal likelihood as the objective function.
Let q(_t) denote an arbitrary probability measure which will eventually approximate the filtering density p(_t |_≤ t).
From Eq. (<ref>), we can rearrange the marginal log-likelihood as

    log p(_t|_<t)                                                                                                                                                           
        = logp(_t|_t) p(_t|_<t)/p(_t|_≤ t)    for any _t
        = _q(_t)[logp(_t|_t) p(_t|_<t)/p(_t|_≤ t)]     the marginal is constant to q(_t)
        = _q(_t)[
    			logp(_t|_t) p(_t|_<t) q(_t)/p(_t|_≤ t) q(_t)]
    		
        =
    		_q(_t)[ log p(_t|_t) ]_reconstruction log-likelihood
    		- (q(_t)    p(_t|_<t))
    		+ (q(_t)    p(_t|_≤ t))_variational gap #1
       ≥_q(_t)[
    				log p(_t|_t)
    				]
    		+ (q(_t))_entropy
    		+ _q(_t)[
    			log p(_t|_<t) ]
    		
        = _q(_t)[ log p(_t|_t) ]
    		+ (q(_t))
    		+ _q(_t)[
    			log_p( _t-1|_<t)[
    				p(_t |_t-1)
    				]
    			]
    		
        = _q(_t)[ log p(_t|_t) ]
    		+ (q(_t))
    		+ _q(_t)[
    			_p(_t-1 | _<t)[
    				log p(_t|_t-1)
    				]]
          +
    		_q(_t)[
    				(p(_t-1|_<t)  p(_t-1|_t,_<t))
    				]
    		_variational gap #2
       ≥_q(_t)[ log p(_t|_t) ]
    		+ (q(_t))
    		+ _q(_t)_ p( _t-1|_<t) [
    			log p(_t |_t-1)
    			]

where  denotes Shannon's entropy and  denotes the Kullback-Leibler (KL) divergence <cit.>.
Maximizing this lower bound would result in a variational posterior
q(_t) ≈ p(_t |_≤ t) w.r.t. q(_t).
Naturally we plug in the previous step's solution to the next time step, obtaining a loss function suitable for recursive estimation:

    ℒ   _q(_t)[ log p(_t|_t) ]
    		+ (q(_t))
    		+ _q(_t)_ q( _t-1 ) [
    				log p(_t |_t-1)
    				]
    		_dynamics log-likelihood

This also results in consistent q(_t) for all time steps as they are in the same family of distribution.

Meanwhile, as it is aimed to jointly estimate the observation model p(_t |_t) and state dynamics p(_t |_t-1), we achieve online inference by maximizing this objective ℒ w.r.t. their parameters (omitted for brevity) and the variational posterior distribution q(_t) simultaneously provided that q(_t-1) takes some parameterized form and has been estimated from the previous time step.
Maximizing the objective ℒ is equivalent to minimizing the two variational gaps:
(1) the variational filtering posterior has to be close to the true filtering posterior, and
(2) the filtering posterior from the previous step needs to be close to p(_t-1|_t, _<t).
Note that the second gap is invariant to q(_t) if p(_t-1|_t, _<t) = p(_t-1|_<t), that is, the one-step backward smoothing distribution is identical to the filtering distribution.

On the flip side, intuitively, there are three components in ℒ that are jointly optimized: (1) reconstruction log-likelihood which is maximized if q(_t) concentrates around the maximum likelihood estimate given only _t, (2) the dynamics log-likelihood which is maximized if q(_t) concentrates at around the maximum of _q(_t-1)[ log p(_t|_t-1) ], and (3) the entropy that expands q(_t) and prevents it from collapsing to a point mass.

In order for this recursive estimation to be real-time, we choose q(_t) to be a multivariate normal with diagonal covariance 𝒩(μ_t, _t) where μ_t is the mean vector and _t is the diagonal of the covariance matrix in this study.
Moreover, to amortize the computational cost of optimization to obtain the best q(_t) on each time step, we employ the variational autoencoder architecture <cit.> to parameterize q(_t) with a recognition model.
Intuitively, the recognition model embodies the optimization process of finding q(_t), that is, it performs an approximate Bayesian filtering computation (in constant time) of Eq. (<ref>) according to the objective function ℒ.
We use a recursive recognition model that maps q(_t-1) and _t to q(_t).
In particular, the recognition model takes a deterministic recursive form:

    [μ_t, _t] =
    	h(_t, u_t-1, μ_t-1, _t-1)

Specifically h takes a simple the form of the multi-layer perceptron (MLP) <cit.> in this study, and we refer to its parameters as the recognition model parameters.
Note that the recursive architecture of the recognition model reflects the Markovian structure of the assumed dynamics (c.f., smoothing networks often use bidirectional recurrent neural network (RNN) <cit.> or graphical models <cit.>).

The expectations appearing in the reconstruction log-likelihood and dynamics log-likelihood are not always tractable in general. For those intractable cases, one can use the reparameterization trick and stochastic variational Bayes <cit.>: rewriting the expectations over q as expectation over a standard normal random variable, i.e., μ_t + _t^1/2𝒩(0, 1), and using a single sample for each time step.
Hence, in practice, we optimize the following objective function (the other variables and parameters are omitted for brevity),

    ℒ̂ =
    	log p(_t|_t, θ) + _q(_t)log p(_t |_t-1, θ)
    	+ H(q(_t))

where _t and _t-1 represent random samples from q(_t) and q(_t-1) respectively.
Note that the remaining expectation over q(_t) has closed form solution under our Gaussian state noise, ϵ_t, assumption.
Thus, our method can handle arbitrary observation and dynamic model unlike dual form nonlinear Kalman filtering methods that usually suffer difficulties in sampling, e.g. transforming Gaussian random numbers into point process observations.

The schematics of the proposed inference algorithm is summarized by two passes in Figure <ref>.
In the forward pass, the previous latent state generates the new state through the dynamic model, and the new state transforms into the observation through the observation model. In the backward pass, the recognition model recovers the current latent state from the observation, and the observation model, recognition model and dynamic model, are updated by backpropagation.


Algorithm <ref> is an overview of the recursive estimation algorithm.
Denoting the set of all parameters by Θ of the observation model, recognition model and dynamic models, the objective function in Eq. (<ref>) is differentiable w.r.t. Θ, and thus we employ empirical Bayes and optimize it through stochastic gradient ascent (using Adam <cit.>).
We outline the algorithm for a single vector time series, but we can filter multiple time series with a common state space model simultaneously, in which case the gradients are averaged across the instantiations.
Note that this algorithm has constant time and space complexity per time step.

In practice, the measurements _t and input u_t are sampled at a regular interval.
Algorithm <ref> is called after every such observation event, which will return the state estimate along with the parameters and the dynamical system.
One can visualize these for real-time for monitoring, and/or have it streamed to another system for further automated processing (e.g. detect anomalies and raise an alarm or deliver feedback controls).





§ APPLICATION TO LATENT NEURAL DYNAMICS


Our primary applied aim is real-time neural interfaces where a population of neurons are recorded while a low-dimensional stimulation is delivered <cit.>.
State-space modeling of such neural time series have been successful in describing population dynamics <cit.>.
Moreover, models of neural computation are often described as dynamical systems <cit.>.
For example, attractor dynamics where the convergence to one of the attractors represents the result of computation <cit.>.
Here we propose a parameterization and tools for visualization of the model suitable for studying neural dynamics and building neural interfaces <cit.>.
In this section, we provide methodological details for the results presented in the next section.




 §.§ Parameterization of the state space model


Having in mind high-temporal resolution neural spike trains where each time bin has at most one action potential per channel, we describe the case for point process observation.
However, note that our method generalizes to arbitrary observation likelihoods that are appropriate for other modalities, including calcium imaging or local field potentials.
The observed point process time series _t is a stream of sparse binary vectors.
All experiments of point process observation were binned finely so that the time bins contain one event each at most in this study.

Our observation model, Eq. (<ref>), assumes that the observation vector _t is sampled from a probability distribution P determined by the latent state _t though a linear-nonlinear map possibly together with extra parameters at each time t,

    _t ∼ P(g(_t + ))

where gℝ→ℝ is a point-wise map.
We use the canonical link g(·) = exp(·) for Poisson likelihood and identity for Gaussian likelihood in this study.
Note that this observation model is not identifiable since _t = () (^-1_t) where  is an arbitrary invertible matrix.

We normalize the loading matrix  in each iteration.
It is straightforward to include more additive exogenous variables, history-filter for refractory period, coupling between processes, and stimulation artifacts <cit.>.

For state dynamic model, we propose to use a specific additive parameterization with  state transition function and input interaction as a special case of Eq. (<ref>),
 
	
    _t+1    = _t + f(_t) + _t _t + ϵ_t+1
    
    		f(_t)                     = ϕ(_t)                                     
    _0, ϵ_t    ∼𝒩(0, σ^2 )


where ϕ(·) is a vector of r continuous basis functions, i.e.
ϕ(·) = (ϕ_1(·), …, ϕ_r(·))^⊤,  is the weight matrix of the radial basis functions, and _t is the interaction with the input _t. The interaction _t can be globally linear, parameterized as a matrix independent from _t, or locally linear, parameterized as a matrix-valued function of _t using also RBF networks. i.e. ((_t)) = _B ϕ(_t) where _B is the respective weight matrix.
In this study, we used squared exponential radial basis functions <cit.>,

    ϕ_i() = exp(-γ_i ‖ - _i ‖_2^2)

with centers _i and corresponding inverse squared kernel width γ_i.
Though the dynamics can be modeled by other universal approximators such as percepton and RNN, we chose the radial basis function network for the reasons of non-wild extrapolation (zero velocity when the state is far away from data) and fast computation.



The time complexity of our algorithm is 𝒪(mpr + n(m + p + q)) where n, m, p, q, r denote the dimensions of observation, latent space, input, the numbers of hidden units and radial basis functions for this specific parameterization.
Practically to achieve realistic computation time for real-time applications in neuroscience, the number of radial basis functions and hidden units are constrained by the requirement.
Note that the time complexity does not grow with time that enable efficient online inference.

If we compare this to an efficient offline algorithm such as PLDS <cit.> run repeatedly for every new observation (“online mode”), its time complexity is 𝒪(t · (m^3 + mn)) per time step at time t which grows as time passes, making it impractical for real-time application.




 §.§ Phase portrait analysis

Phase portrait displays key qualitative features of dynamics, and with a little bit of training, it provides a visual means to interpreting dynamical systems.
The law that governs neural population dynamics captured in the inferred function f() directly represents the velocity field of an underlying smooth dynamics (<ref>) in the absence of input <cit.>.
In the next section, we visualize the estimated dynamics as phase portrait which consists of the vector field, example trajectories, and estimated dynamical features (namely fixed points) <cit.>.
We can numerically identify candidate fixed points ^∗ that satisfy f(^∗) ≈ 0.
For the synthetic experiments, we performed an affine transformation to orient the phase portrait to match the canonical equations in the main text when the simulation is done through the proposed observation model if the observation model is unknown and estimated.






 §.§ Prediction

For state space models, we can predict both future latent trajectory and future observations.
The s-step ahead prediction can be sampled from the predictive distributions:

	
    p(_t+1:t+s|_≤ t) =    𝔼_q(_t) [p(_t+1:t+s|_t)]
    		
    
    		p(_t+1:t+s|_≤ t) =    𝔼_p(_t+1:t+s|_≤ t) [p(_t+1:t+s|_t+1:t+s)]


given estimated parameters by current time t without seeing the data _t+1:t+s during these steps.
In the figures of experiments, we plot the mean of the predictive distribution as trajectories.



§ EXPERIMENTS ON THEORETICAL MODELS OF NEURAL COMPUTATION


We demonstrate our method on a range of nonlinear dynamical systems relevant to neuroscience.
Many theoretical models have been proposed in neuroscience to represent different schemes of computation.
For the purpose of interpretable visualization, we choose to simulate from two or three dimensional dynamical systems. We apply the proposed method to four such low-dimensional models: a ring attractor model as a model of internal head direction representation, an nonlinear oscillator as a model of rhythmic population-wide activity, a biophysically realistic cortical network model for a visual discrimination experiment and a chaotic attractor.

In the synthetic experiments, we first simulated state trajectories by respective differential equations, and then generated either Gaussian or point process observations (to mimic spikes) via Eq. (<ref>) with corresponding distributions. The parameters  and  were randomly drawn, and they were constrained to keep firing rate <60 Hz on average for realistic spiking behavior. All observations are spatially 200-dimensional unless otherwise mentioned.
We refer to their conventional formulations under different coordinate systems, but our simulations and inferences are all done in Cartesian coordinates.
Note that we focus on online learning in this study and always train our model with streaming data, even while comparing with offline methods.

The approximate posterior distribution is defined recursively in Eq. (<ref>) as diagonal Gaussian with mean and variance determined by corresponding observation, input and previous step via a recurrent neural network.
We used a one-hidden-layer MLP in this study. Typically the state noise variance σ^2 is unknown and has to be estimated from data. To be consistent with Eq. (<ref>), we set the starting value of σ^2 to be 1, and hence μ_0 = 0, _0 =.
We initialize the loading matrix  by factor analysis, and column-wisely normalize it by ℓ_2 norm every iteration to keep the system identifiable.






 §.§ Ring attractor


Continuous attractors are often used as models for neural representation of continuous variables <cit.>. For example, a bump attractor network with ring topology is proposed as the dynamics underlying the persistently active set of neurons that are tuned for the angle of the animal's head direction <cit.>. Here we use the following 2-variable reduction of the ring attractor system:
First, we study the following two-variable ring attractor system:

    τ_r ṙ    = r_0 - r 
    τ_φφ̇    = I

where φ represents the direction driven by input I, and r is the radial component representing an internal circular variable, such as head direction.
We simulated 100 trajectories (1000 steps) with step size Δ t = 0.1, r_0=1, τ_r = 1, τ_φ=1 with Gaussian state noise (std=0.005) added each step.
Though the ring attractor is defined in polar coordinate system, we transformed it into Cartesian system for simulation and training. In simulation we used strong input (tangent drift) to keep the trajectories flowing around the ring clockwise or counter-clockwise. The point process observations were generated by passing the states through a linear-nonlinear map (Eq. (<ref>)) and sampling from a Poisson distribution.
We streamed the observations into the proposed algorithm that consists of point process likelihood, dynamic model with 20 radial basis functions and locally linear input interaction in Eq. (<ref>) and a recognition MLP with 100 hidden units.



Figure <ref>A illustrates one latent trajectory (black) and its variational posterior mean (blue).
These two trajectories start at green circle and diamond respectively and end at the red markers. The inference starts near the center (origin) that is relatively far from the true location because the initial posterior mean is set at zero. The final states are very close which implies that the recognition model works well. Figure <ref>B shows the reconstructed velocity field by the model.
We visualized the velocity as colored directional streamlines. We can see the velocity toward the ring attractor and the speed is smaller closer to the ring.
The model also identifies a number of fixed points arranged around the ring attractor via numerical roots finding.
Figure <ref>C shows the distribution of posterior means of all data points in the state space. We have more confidence of the inferred dynamical system in the denser area.

Figure <ref>D shows the three components of Eq. (<ref>) and the objective lower bound clearly, demonstrating the convergence of the algorithm.
We can see each component reaches a plateau within 400 sec.
As the reconstruction and dynamics log-likelihoods increase, the recognition model and dynamical model are getting more accurate while the decreasing entropy indicates the increasing confidence (inverse posterior variance) on the inferred latent states. The average computation time of a joint estimation step is 1.1 ms
(hardware specification: Intel Xeon E5-2680 2.50G Hz, 128GB RAM, no GPU).







 §.§ Nonlinear oscillator


Dynamical systems have been a successful application in the biophysical models of single neuron in neuroscience. We used a relaxation oscillator, the FitzHugh-Nagumo (FHN) model <cit.>, which is a 2-dimensional reduction of the Hodgkin-Huxley model:
We used a  with the following nonlinear state dynamics:

    v̇    = v(a - v)(v - 1) - w + I, 
    ẇ    = bv - cw,

where v is the membrane potential, w is a recovery variable and I is the magnitude of stimulus current in modeling single neuron biophysics.
This model was also used to model global brain state that fluctuates between two levels of excitability in anesthetized cortex <cit.>.
We use the following parameter values a = -0.1, b = 0.01, c = 0.02 and I = 0.1 to simulate 100 trajectories of 1000 steps with step size 0.5 and Gaussian noise (std=0.002).
At this regime, unlike the ring attractor, the spontaneous dynamics is a periodic oscillation, and the trajectory follows a limit cycle.
The point process observations were also sampled via the observation model of the same parametric form as that of the ring attractor example.
We used 20 radial basis functions for dynamic model and 100 hidden units for recognition model. While training the model, the input was clamped to zero, and expect the model to learn the spontaneous oscillator.

We compare the state estimation with the standard particles filtering (PF) which are powerful online methods theoretically capable of producing arbitrarily accurate filtering distribution.
We run two variants of the particle filter with different proposal distributions.
One used diffusion as the proposal, i.e. _t = _t-1 + ϵ_t where  is the vector of state variables v and w, and the other, a.k.a. bootstrap particle filter <cit.>, used the true dynamics in Eq. (<ref>).
We provided the true parameters for the observation model and noise term to PF which gives them an advantage.
Both particle filters and VJF were run on 50 realizations of 5000-step long observation series.
Figure <ref> shows the root mean squared deviations (RMSE) (mean and standard error over 50 realizations). It is expected that the bootstrap particle filter outperformed the diffusion particle filter since the former utilized the true dynamics. One can see the state estimation by VJF improved as learning carrying on and eventually outperformed both particle filters. Note that VJF had to learn the parameters of likelihood, dynamic model and recognition model during the run.
We varied the number of RBFs (20 and 30) but the results are not substantially different.



We also reconstructed the phase portrait (Fig. <ref>B) comparing to the truth (Fig. <ref>C). The two dashed lines are the theoretical nullclines of the true model on which the velocity of corresponding dimension is zero. The reconstructed field shows a low speed valley overlapping with the nullcline especially on the right half of the figure. At the intersection of the two nullclines there is an unstable fixed point. We can see the identified fixed point is close to the intersection. As most of the trajectories lie on the oscillation path (limit cycle) with merely few data points elsewhere, the inferred system shows the oscillation dynamics similar to the true system around the data region. The difference mostly happens in the region far from the trajectories because of the lack of data.

We run a long-term prediction using VJF without seeing the future data _t+1:T during these steps (T=1000 steps = 1 sec) beginning at the final state of training data. We show the truth and prediction in figure <ref>D. The upper row is the true latent trajectory and corresponding observations. The lower row is the filtered trajectory and prediction by the proposed method. The light-colored parts are the 500 steps of inference before prediction and the solid-colored parts are 1000 steps truth and prediction. We also show the sample observations from the trained observation model during the prediction period.

One of the popular latent process modeling tools for point process observation that can make prediction is the Poisson Linear Dynamical System (PLDS) <cit.> which assumes latent linear dynamics.
We compared PLDS fit with EM on its long-term prediction on both the states and spike trains (Fig. <ref>). This demonstrates the nonlinear dynamical model outperforming the linear model even in the unfair online setting.



To compare to the methods with nonlinear dynamical models, we also run latent factor analysis via dynamical systems (LFADS) <cit.> offline using the same data. LFADS implements its dynamical model with the gated recurrent unit (GRU) <cit.> that requires high dimensions. For this 2D system, we tried different GRU dimensionalities. We made minimal changes to its recommended setting including only the generator dimensionality, batch and no controller. The result shows that LFADS requires much higher dimension than the true system to capture the oscillation (Fig. S1). (The figure of its inferred trajectories is shown in the supplement.) We report the fitted log-likelihood per time bin -0.1274, -0.1272 and -0.1193 for 2D, 20D and 50D GRU respectively. In comparison, the log-likelihood of the proposed approach is -0.1142 with a 2D dynamical model (higher the better).




 §.§ Fixed point attractor for decision-making


Perceptual decision-making paradigm is a well-established cognitive task where typically a low-dimensional decision variable needs to be integrated over time, and subjects are close to optimal in their performance.
To understand how the brain implements such neural computation, many competing theories have been proposed <cit.>.
We test our method on a simulated biophysically realistic cortical network model for a visual discrimination experiment <cit.>.
In the model, there are two excitatory subpopulations that are wired with slow recurrent excitation and feedback inhibition to produce attractor dynamics with two stable fixed points (Fig. <ref>A).
Each fixed point represents the final perceptual decision, and the network dynamics amplify the difference between conflicting inputs and eventually generates a binary choice.

Note that, different from the former examples that use a linear-nonlinear map of latent states, the point process observations (spikes) of this experiment were directly sampled from the spiking neural network [The detail of the spiking neural network can be found in <cit.> and the code can be found at <https://github.com/xjwanglab/book/tree/master/wang2002>.](1 ms binwidth) that was governed by its own high-dimensional intrinsic dynamics. It is filling the gap between fully specified state space models and real neuron populations.

We subsampled 480 selective neurons out of 1600 excitatory neurons from the simulation to be observed by our algorithm.
The simulated data is organized into decision-making trials where each trial lasts for 2 sec and with different strength of visual evidence, controlled by “coherence”.
Our method with 20 radial basis functions learned the dynamics from 140 training trials (20 per coherence level c, c=-1, -0.2, -0.1, 0, 0.1, 0.2, 1).

Figure <ref>C shows the velocity field at zero coherence stimulus as colored streamlines.
Note that our approach did not have prior knowledge of the network dynamics as the mean-field reduction <cit.> in Figure <ref>B.
Although the absolute arrangement is dissimilar, the topology and relation of the five identified fixed points show correspondence with the mean-field reduction.
The inference was completely data-driven (partial observation of spike trains) while the mean-field method required knowing the true dynamical model of the network and careful approximation by <cit.>.
We showed that our method can provide a qualitatively similar result to the theoretical work which reduces the dimensionality and complexity of the original network.






 §.§ Chaotic dynamics


Chaotic dynamics (or edge-of-chaos) have been postulated to support asynchronous states in the cortex, and neural computation over time by generating rich temporal patterns <cit.>.
We consider the 3-dimensional standard Lorenz attractor as an example chaotic system to demonstrate the flexibility of our method.
We simulated 216 latent trajectories from:

    ẋ = 10(y - x),   ẏ = x(28 - z) - y,   ż = x y - 8/3 z.

The each coordinate of the initial states are on the uniform grid of 6 values in [-50, 50] inclusively, of which the combination results in 216 unique states.
We discarded the first 500 transient steps of each trajectory and then use the following 1000 steps.
We generated 200-dimensional Gaussian observations driven by the trajectories.
Figure <ref>A shows estimated latent trajectory and the ground truth. One can see that the estimation lies in a similar manifold.
In addition, we predicted 500 steps of future latent states without knowing the respective observations. Figure <ref>B shows 4 predicted trajectories starting from different initial states. One can see that the inferred system could generate qualitatively similar trajectory at most initial states but not perfectly for the true system is chaotic.







 §.§ Nonstationary system

Another feature of our method is that its state dynamics estimate never stops.
As a result, the algorithm is adaptive, and can potentially track slowly varying (nonstationary) latent dynamics.
To test this feature, we compared a dual EKF and the proposed approach on nonstationary linear dynamical system. A spiral-in linear system was suddenly changed from clockwise to counter-clockwise at the 2000th step and the latent state was perturbed (Fig. <ref>).
To adapt EKF, we used Gaussian observations that were generated through linear map from 2-dimensional state to 200-dimensional observation with additive noise (𝒩(0, 0.5)). To focus on the dynamics, we fixed all the parameters except the transition matrix for both methods, while our approach still has to learn the recognition model in addition.
Figure <ref> shows that our approach achieved better online performance as dual EKF in this experiment.






§ REAL NEUROPHYSIOLOGICAL APPLICATION


We applied the proposed method to a large scale recording to validate that it picks up meaningful dynamics. The dataset <cit.> consists of 148 simultaneously recorded single units from the primary cortex (V1) while directional drifting gratings were presented to an anesthetized monkey for around 1.3s per trial (Fig. <ref>A). We used the spike trains from 63 well-tuned units. The spike times were binned with a 1ms window (max 1 spike per bin).
There is one continuous circular variables in the stimuli space: temporal phase of oscillation induced by the drifting gratings.

A partial warm-up helps with the training. We chose a good initialization for the observation model, specifically the loading matrix and bias. There are 72 motion directions in total, each repeated 50 trials.
We used the trials corresponding to 0 deg direction to initialize the observation model with dimensionality reduction methods such as variational latent Gaussian processes <cit.>, and then trained VJF with a 2D dynamic model fully online on the trials corresponding to 180 deg direction that it had not seen before. Since we do not have long enough continuously-recorded trials, we concatenated the trials (equivalent to 500s) as if they were continuously recorded to mimic an online setting.
As expected, Figure <ref>B and C shows the inferred dynamical system is able to implement the oscillation.
The two goodness of fit measures (log-likelihood and ELBO) in Figure <ref>D shows that our method benefits from but does not necessarily require such a warm-up. The model with warm-up initialization has better starting goodness of fit than the random initialized model but the random initialized model eventually achieved similar goodness of fit with adequate amount of data.







§ DISCUSSION


Neurotechnologies for recording the activity of large neural populations during meaningful behavior provide exciting opportunities for investigating the neural computations that underlie perception, cognition, and decision-making.
However, the datasets provided by these technologies currently require sophisticated offline analyses that slow down the scientific cycle of experiment, data analysis, hypothesis generation, and further experiment.
Moreover, in closed-loop neurophysiological setting, real-time adaptive algorithms are extremely valuable <cit.>.

To fulfill this demand, we proposed an online algorithm for recursive variational Bayesian inference that simultaneously performs system identification and state filtering under the framework of state space modeling, in hope that it can greatly impact neuroscience research and biomedical engineering.
There is no other method capable of all features, hence we compared several methods in different measures, often giving them the advantage.
We showed that our proposed method consistently outperforms the state-of-the-art methods.

Using the language of dynamical systems, we interpret the target system not via model parameters but via dynamical features: fixed points, limit cycles, strange attractors, bifurcations and so on.
In our current approach, this interpretation heavily relies on visual inspection of the qualitative nonlinear dynamical system features.
In contrast, most popular state space models assume linear dynamics <cit.> which is appropriate for smoothing latent states, but not expressive enough to recover the underlying vector field.
Recently the Koopman theory that allows representation of general nonlinear dynamics as linear operators in infinite dimensional spaces <cit.> has gained renewed interest in modeling nonlinear dynamics.
Although elegant in theory, in practice, however, the Koopman operators need to be truncated to a finite dimensional space with linear dynamics <cit.>.
We note that the resulting linear models do not allow for topological features such as multiple isolated fixed points, nonlinear continuous attractors, stable limit cycles—features critical for non-trivial neural computation.


Our algorithm is highly flexible and general—it allows a wide range of observation models (likelihoods) and dynamic models, is computationally tractable, and produces interpretable visualizations of complex collective network dynamics.
Our key assumption is that the dynamics consists of a continuous and slow flow, which enable us to parameterize the velocity field directly.
This assumption reduces the complexity of the nonlinear function approximation, thus it is easy to identify the fixed/slow points.
Specifically we chose the radial basis function network to model the dynamics for our experiments, which regularizes and encourages the dynamics to occupy a finite phase volume around the origin.





Our method has a number of hyperparameters.
In the experiments, the differentiable hyperparameters were learnt via gradient descent while the selection of the other hyperparameters were made simple.
In general, our method was robust; Perturbing the number of RBFs did not produce qualitatively different results (Fig. <ref>).
<cit.> discussed growing radial basis function network adaptively which could be incorporated in our method to enable online tuning of the number of RBFs.
The depth and width of neural networks were chosen empirically to improve the interpretability of resulting dynamical systems, but tuning did not result in large changes in the results.


This work opens many avenues for future work.
One direction is to apply this model to large-scale neural recording from a behaving animal.
We hope that further development would enable on-the-fly analysis of high-dimensional neural spike train during electrophysiological experiments.
Clinically, a nonlinear state space model provides a basis for nonlinear feedback control as a potential treatment for neurological diseases that arise from diseased dynamical states.



hunsrt



