The canonical measure for quantifying information transmission via time-varying signals is the information transmission rate [1], [2]. It is defined as the rate at which the path mutual information between the input and output trajectories of the system increases with the trajectory duration. It quantifies the speed at which distinct messages are transmitted through the system, taking into account the correlations that are present in the input and output signals. The information rate has been used to to quantify biochemical signaling performance [3]–[5], to perform model reduction [6], to detect the causality of interactions [7], [8], to test for nonlinearities in time series [9], to assess dependencies between stock prices or market indices over time [10]–[12], or to quantify information exchange between different regions of the brain [13], [14]. In the absence of feedback the rate also equals the multi-step transfer entropy [15], [16].

Yet, computing the path mutual information between input and output trajectories is a notoriously difficult problem because trajectories are high-dimensional objects. Conventional methods to compute the mutual information rely on non-parametric estimates of the joint probability distribution of input and output, e.g., via histograms or kernel density methods [17], [18]. These methods are, however, infeasible for high-dimensional data [19], [20]. More advanced non-parametric estimators such as the k-nearest-neighbor (KNN) estimator [21] are better suited for high-dimensional data, but suffer from uncontrolled biases as the dimensionality of data increases [20], [22], [23]. Consequently, the rate is often computed using approximate schemes, such as the Gaussian framework [2], [24], moment-closure approximations [5], [25], or techniques limited to specific systems [26]–[28]. Schemes based on machine learning have also been introduced, such as decoding-based approaches [20], Information Noise-Contrastive Estimation (InfoNCE) [29] and the Difference-of-Entropies (DoE) estimator [30].

We recently presented Path Weight Sampling (PWS) [31], a computational technique that makes it possible to compute the information transmission rate exactly for any stochastic model, be it described by a master equation or stochastic differential equations. The principal idea is to use the stochastic model to evaluate the path likelihood, i.e. the conditional probability of an individual output trajectory for a given input trajectory, and then average this quantity via Monte Carlo sampling in trajectory space to obtain the path mutual information. While the scheme is exact and does not rely on uncontrolled approximations, it does require a stochastic model of the system of interest to evaluate the path likelihood. The method therefore cannot be directly applied to experimental time series data.

In this manuscript, we show how PWS can be combined with Machine Learning (ML) to obtain the information rate directly from experimental data. The new idea of our technique, called ML-PWS, is to first use machine learning to develop a generative model which describes the time series data, allowing us to compute the path likelihood. As in conventional PWS, this path likelihood is then averaged to obtain the information transmission rate. Here, we develop the generative model based on a neural autoregressive sequence prediction, as is used e.g. in speech synthesis [32] or text generation [33], but the principal idea is more generic. We demonstrate the power of ML-PWS by applying it to synthetic time series data generated from a non-linear model, and by comparing the rate thus obtained against the ground-truth result obtained by applying PWS directly to this model. We further illustrate the utility of ML-PWS by applying it to neuronal time-series data [34].

Path Weight Sampling—The information rate \(R(\mathcal{S}, \mathcal{X})\) is defined as the speed at which the mutual information \(I(S_{1:n}, X_{1:n})\) between an input trajectory \(S_{1:n}\) and output trajectory \(X_{1:n}\) increases with the trajectory duration \(T=n \delta t\) \[R(\mathcal{S}, \mathcal{X}) = \lim_{n\to\infty} \frac{1}{n \delta t} I(S_{1:n}, X_{1:n}),\] where \[I(S_{1:n},X_{1:n}) = \left \langle \ln\frac{\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n})}{\mathop{\mathrm{\mathrm{P}}}(x_{1:n})} \right\rangle_{\mathop{\mathrm{\mathrm{P}}}(s_{1:n},x_{1:n})} \label{eq:trajectory95mi} \,.\tag{1}\] \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n})\) is the marginal distribution of trajectories \(X_{1:n}\), and \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n})\) is the conditional distribution of the output given the input trajectory, discretized with a timestep \(\delta t\). PWS is based on the idea that for systems described by a master equation or a stochastic differential equation, the path likelihood \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n})\) can be computed on the fly, which makes it possible to estimate the mutual information as a Monte Carlo average over trajectory space: \[\begin{align} \hat{I}_{\rm MC} \approx \frac{1}{N} \sum^N_{i=1} \ln \frac{\mathop{\mathrm{\mathrm{P}}}(x^i_{1:n} | s^i_{1:n})}{\mathop{\mathrm{\mathrm{P}}}(x^i_{1:n})},\label{eqn:IMC} \end{align}\tag{2}\] with \[\begin{align} \mathop{\mathrm{\mathrm{P}}}(x^i_{1:n}) \approx \frac{1}{M} \sum^M_{j=1} \mathop{\mathrm{\mathrm{P}}}(x^i_{1:n} | s^j_{1:n}) .\label{eqn:PxMC} \end{align}\tag{3}\] In the brute-force version of PWS, called Direct PWS (DPWS), the mutual information is indeed obtained via two nested Monte Carlo averages, in which \(N\) pairs of \((s^i_{1:n}, x^i_{1:n})\) are generated from \(\mathop{\mathrm{\mathrm{P}}}(s_{1:n}, x_{1:n})\), and for each output trajectory \(x^i_{1:n}\) \(M\) input trajectories \(s^j_{1:n}\) are generated from the input distribution \(\mathop{\mathrm{\mathrm{P}}}(s_{1:n})\) [31]. The main challenge is the computation of \(\mathop{\mathrm{\mathrm{P}}}(x^i_{1:n})\), Eq. 3 , for which efficient techniques are presented in Ref. [31] .

While PWS is an exact Monte Carlo scheme, providing an unbiased statistical estimate of the mutual information, it does require a stochastic model from which the path likelihood \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n})\) can be obtained, see Eqs. 2  and 3 . While the list of systems for which an accurate mechanistic or a phenomenological model is available is rapidly growing, in many cases no such model is available. To overcome this problem, ML-PWS first learns a stochastic model from experimental time-series data, enabling the computation of the path likelihood \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n})\), from which the mutual information is then computed using PWS.

Machine Learning PWS—Our generative stochastic model is based on autoregressive neural networks which have been used for sequence prediction [35], image modeling [32], natural language processing [36] and other tasks. All of these models factorize the joint probability of a sequence \(x_{1:n}\) as \[\begin{align} \mathop{\mathrm{\mathrm{P}}}(x_{1:n}) = \prod_{i=1}^n \mathop{\mathrm{\mathrm{P}}}(x_i | x_{1:i-1}). \end{align}\] To obtain the path likelihood \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n})\), we also need to include the conditioning on the input signal, taking into account the causal relationship between \(s_{1:n}\) and \(x_{1:n}\). Since \(s_{1:n}\) and \(x_{1:n}\) represent time-series, \(x_i\) does not depend on any of the future inputs \(s_{i+1},\ldots,s_n\). Therefore, instead of conditioning on the full input \(s_{1:n}\), we only condition on \(s_{1:i}\): \[\begin{align} \mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n}) &= \prod^n_{i=1} \mathop{\mathrm{\mathrm{P}}}(x_i|x_{1:i-1}, s_{1:i}). \label{eqn:autoregressive_conditional_density} \end{align}\tag{4}\] Up to this point, no approximation has been made. To make further progress, we model each conditional distribution \(\mathop{\mathrm{\mathrm{P}}}(x_i|x_{1:i-1}, s_{1:i})\) in Eq. 4 using a parametric distribution. An autoregressive neural network then predicts the values of the parameters of this distribution at time step \(i\) based on the observed sequence of inputs \(s_{1:i}\) and past outputs \(x_{1:i-1}\).

To train the model, we exploit that the training data consists of \(N\) pairs of trajectories \((s^k_{1:n}, x^k_{1:n})\) for \(k=1,\ldots,N\) that represent independent draws from the true data distribution. The loss function is then given by \[\mathcal{L}(\boldsymbol{\theta}) = -\sum^N_{k=1} \ln \mathop{\mathrm{\mathrm{P}}}(x^k_{1:n}|s^k_{1:n}, \boldsymbol{\theta}) \,. \label{eqn:loss_function}\tag{5}\] and is minimized with respect to the neural weights \(\boldsymbol{\theta}\).

The above ML approach makes it possible to develop a generative model \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n})\), which can be combined with PWS to obtain the mutual information (Eqs. 2  and 3 ). Yet, ML can also be used to optimize PWS itself, by improving its most challenging step, which is the computation of \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n})\) (see Eq. 3 ). The improvement is based on the observations that (a) \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n})\) is also given by \[\begin{align} \mathop{\mathrm{\mathrm{P}}}(x_{1:n}) = \left\langle\frac{\mathop{\mathrm{\mathrm{P}}}(s_{1:n})\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n})}{q(s_{1:n}|x_{1:n})}\right\rangle_{q(s_{1:n}|x_{1:n})}, \label{eqn:ml_importance_sample_mc} \end{align}\tag{6}\] where \(q(s_{1:n}|x_{1:n})\) is an importance sampling distribution over \(s_{1:n}\), and (b) the variance in the estimate of \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n})\) is minimized when \(q(s_{1:n}|x_{1:n})\) equals the true posterior distribution \(\mathop{\mathrm{\mathrm{P}}}(s_{1:n}|x_{1:n})\). In the spirit of variational autoencoders, the key idea is then to train a second neural network that parametrizes \(q(s_{1:n}|x_{1:n})\), the inference model, by minimizing the Kullback-Leibler difference between \(q(s_{1:n}|x_{1:n})\) and \(\mathop{\mathrm{\mathrm{P}}}(s_{1:n}|x_{1:n})\), via the Evidence Lower Bound Objective (ELBO) [37]. Although the estimate of \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n})\) is always unbiased, independent of the choice of \(q(s_{1:n}|x_{1:n})\), optimizing \(q(s_{1:n}|x_{1:n})\) using this approach improves the efficiency of PWS. For this inference model, we use an autoregressive normalizing flow, which is optimized jointly with the forward network [see Supplementary Material (SM)].

Benchmarks—We test ML-PWS by applying it to synthetic data generated from a minimal non-linear model and comparing the result against the exact result obtained by directly applying PWS to the same model. The model to generate the synthetic data \((s_{1:50}, x_{1:50})\) consists of an Auto-Regressive input of order 3, AR(3), and a nonlinear output with a nonlinearity that is characterized by a gain parameter \(\gamma\), see [fig:ARtest](a) (see SM).

To obtain the information rate using ML-PWS, we use this data to train a Gaussian autoregressive neural network, where each conditional probability distribution \(\mathop{\mathrm{\mathrm{P}}}(x_i|x_{1:i-1}, s_{1:i})\) in Eq. 4 is Gaussian with mean \(\mu_i(x_{1:i-1}, s_{1:i})\) and standard deviation \(\sigma_i(x_{1:i-1}, s_{1:i})\). While each conditional distribution is Gaussian, the whole sequence is not, due to the nonlinear nature of the neural network.

In [fig:ARtest](b) we compare the mutual information estimates of ML-PWS for this dataset against various benchmarks. The green dots display the ML-PWS estimate of the mutual information \(I(S_{1:50}, X_{1:50})\) as a function of the gain \(\gamma\). As expected, for small \(\gamma\), the mutual information grows with \(\gamma\) as the gain enhances the signal-to-noise ratio. For larger values of \(\gamma\), we observe a saturation and even a decline in the information rate due to the saturation effect of the logistic function. This behavior is indicative of the nonlinearity of the system.

First, we compare the ML-PWS result against the “ground truth” mutual information obtained by applying PWS directly to the model. [fig:ARtest](b) shows that the ML-PWS result matches the ground truth very well across all values of \(\gamma\). This demonstrates that the autoregressive neural network can accurately learn the stochastic dynamics of the nonlinear model and reliably estimate the path likelihood, which is required for the Monte Carlo estimate of mutual information. These results confirm that combining PWS with machine learning is a feasible and promising approach for computing the mutual information rate in complex nonlinear systems.

Second, we compute the mutual information using the Gaussian approximation, which is widely used for directly estimating mutual information rates from time-series data [24]. To make a fair comparison of our ML-PWS technique against the Gaussian approximation, we use the same dataset for the Gaussian approximation as for training the ML model. We refer to this benchmark as “Gaussian I” in [fig:ARtest](b). The Gaussian approximation suffers from two sources of bias: a finite sample size bias and a bias arising from the assumption of linearity which does not hold at large \(\gamma\). To distinguish between these two sources of bias, we created another benchmark, called “Gaussian II”. It is similar to Gaussian I but is obtained not from 1000 but 100000 trajectory pairs \((s_{1:50}, x_{1:50})\), thus allowing for very precise estimates of the covariance matrices required for the Gaussian approximation, effectively eliminating the sample size bias.

[fig:ARtest](b) shows that both Gaussian descriptions closely match the ground truth for \(\gamma \lesssim 1\), but deviate for higher gain. In the latter regime, the Gaussian model fails to correctly capture the nonlinear dynamics of the system, in contrast to ML-PWS. The panel also shows that the estimate of Gaussian I is consistently higher than that of Gaussian II because of the finite sample size bias. Yet, ML-PWS, applied to the same original data set as that used for Gaussian I, accurately estimates the mutual information. Clearly, ML-PWS is more sample efficient than the Gaussian approximation.

[fig:ARtest](c) compares ML-PWS against two other ML-based schemes, InfoNCE [29] and DoE [30]. ML-PWS shares similarities with DoE: both train a generative model to estimate \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n})\). Yet, while DoE trains a second neural network to estimate \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n})\), ML-PWS obtains \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n})\) via an exact marginalisation of \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n})\mathop{\mathrm{\mathrm{P}}}(s_{1:n})\) using PWS (see Eq. 3 ). NCE estimates the mutual information by training a neural network to maximize the contrast between positive pairs of input and output trajectories drawn from the joint distribution \(\mathop{\mathrm{\mathrm{P}}}(s_{1:n},x_{1:n})\) versus negative pairs drawn from the marginal distributions \(\mathop{\mathrm{\mathrm{P}}}(s_{1:n})\) and \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n})\); in contrast to ML-PWS and DoE, it provides a lower bound on the mutual information. [fig:ARtest](c) shows that the estimates of ML-PWS and DoE are both very close to the exact PWS result for all trajectory lengths. In contrast, NCE only accurately estimates the mutual information for relatively short trajectories. Since the NCE estimate is upper bounded by \(I_{\mathrm{NCE}}(S, X) \leq \ln N\) [29], [38], this method severely underestimates the mutual information for long trajectories.

Application to neural data—To demonstrate the utility of ML-PWS, we apply it to neuronal time series data [34]. The data was recorded from retinal ganglion cells of a salamander, responding to a dark horizontal bar, displayed on a screen, moving stochastically in a vertical direction [34]. The dataset consists of spike trains from 230 neurons, sampled at a temporal resolution of \(0.1 {\rm ms}\). These cells were stimulated by 136 repeated trials of \(T=30.08\,{\rm s}\) during which the same stimulus was presented. The dynamics of the stimulus, the moving bar, is that of an underdamped particle in a harmonic well, which means that the distribution of inputs, \(\mathop{\mathrm{\mathrm{P}}}(s_{1:n})\), is known. For analysis, the time series of the output was discretized at a resolution of \(\Delta t = 10\,{\rm ms}\), yielding \(n = T /\Delta t = 3008\) time bins per trial. The number of spikes per time bin is assumed to be Poisson distributed. These discretized trajectories were used to train a neural network, yielding a generative model that was then combined with PWS to obtain the information rate. We trained a neural network for each cell separately, yielding a generative model for the path likelihood \(\mathop{\mathrm{\mathrm{P}}}(x^{(i)}_{1:n}|s_{1:n})\) for each cell \(i\), from which we then computed the path mutual information \(I(S_{1:n},X^{(i)}_{1:n})\) for each individual cell. In addition, we trained a neural network for the collective response of the population of cells, yielding the path likelihood \(\mathop{\mathrm{\mathrm{P}}}({\boldsymbol{x}}_{1:n}|s_{1:n})\) for the population cells, with \({\boldsymbol{x}}_{1:n}= (x^{(1)}_{1:n},\dots,x^{(50)}_{1:n})\); this yields the mutual information \(I(S_{1:n},{\boldsymbol{X}}_{1:n})\) between the input \(S_{1:n}\) and the collective response \({\boldsymbol{X}}_{1:n}\) of the population of cells.

[fig:NR](a) shows, for 10 randomly selected cells, that the generative model can accurately describe the experimental time series data. In particular, the validation data in the right panel shows that the model can predict the measured spiking rate. [fig:NR](b) shows the path mutual information as a function of trajectory duration, for 50 individual cells (blue lines), their sum (green line), as well as between the input and the collective response of the population of 50 cells (red line). After a transient caused by the delayed response of the system, the mutual information rises linearly with the trajectory duration, with a slope that defines the information transmission rate. It is seen that the mutual information varies substantially from cell to cell, as was also observed in Ref. [34]. Moreover, the mutual information \(I(S_{i:n};{\boldsymbol{X}}_{1:n})\) between the input \(S_{1:n}\) and the collective response \({\boldsymbol{X}}_{1:n}\) (red line) is significantly lower than the sum of the mutual information of the individual cells: correlations between the response of the individual cells lead to redundant coding, lowering the information that is encoded in the collective response [34]. Interestingly, the information rate estimates of ML-PWS, given by the slopes of the curves in Fig. [fig:NR], are about 50% higher than those between the true input and the decoded bar input [34], which indeed provides a lower bound on the true rate. This difference underscores that obtaining precise information estimates requires an accurate scheme like ML-PWS.

Discussion—We demonstrated how autoregressive sequence prediction models can be trained on time-series data to learn a generative model, which can be combined with PWS to compute the information rate. By applying ML-PWS to a nonlinear model, we showed that it provides more accurate mutual information estimates than the Gaussian approximation. While this example serves as a proof of concept, it shows the potential of advanced machine learning techniques to automatically derive stochastic models from experimental data, and to enable the computation of information-theoretic measures for complex, high-dimensional data.

While ML-PWS shares similarities with the DoE estimator [30], ML-PWS has the benefit that \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n})\) is obtained via an exact marginalisation of \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n}) \mathop{\mathrm{\mathrm{P}}}(s_{1:n})\) using PWS, and not by training a second neural network as in DoE; this ensures that the marginal statistics \(x_{1:n}\) are consistent with the input distribution \(\mathop{\mathrm{\mathrm{P}}}(s_{1:n})\), as determined by the setup of the experiment, and the conditional input distribution \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n})\) of the generative model. Moreover, ML-PWS also enables the efficient computation of the mutual information \(I(\boldsymbol{S}^\prime, \boldsymbol{X}^\prime)\) between a input signal \(\boldsymbol{S}^\prime\) with different statistics and the corresponding output \(\boldsymbol{X}^\prime\), without needing to re-train the marginal model. This is particularly useful if one is interested in the channel capacity, as determined by the input distribution that maximizes the mutual information or information rate for the system of interest. In systems without feedback, \(\mathop{\mathrm{\mathrm{P}}}(\boldsymbol{x}|\boldsymbol{s})\) is a property of the system and does not change upon changing the input. The generative model then remains the same, such that ML-PWS can directly recompute the mutual information for different input statistics.

Acknowledgments—This work is part of the Dutch Research Council (NWO) and was performed at the research institute AMOLF. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 885065), and was financially supported by NWO through the “Building a Synthetic Cell (BaSyC)” Gravitation grant (024.003.019). GT acknowledges the support of WWTF grant LS23-026 “Understanding pancreas biology with AI/ML".

Code Availability—The code used to generate the results of this study is openly available [39].

1 End Matter↩︎

Appendix A: Gaussian Autoregressive RNN model—Let the input sequence be \(s_{1:n} = (s_1, \ldots, s_n)\). A recursive neural network (RNN) takes this sequence as well as an initial state \(h_0\in\mathbb{R}^d\) and generates a sequence \(h_{1:n}=(h_1,\ldots,h_n)\) from a recursive relation \[h_i = f_\theta(s_i, h_{i-1}) \label{eq:rnn}\tag{7}\] where \(h_i\in\mathbb{R}^d\) for \(i\in\{1,\ldots,n\}\) and an activation function \(f_\theta: \mathbb{R}\times\mathbb{R}^d\mapsto\mathbb{R}^d\). The activation function \(f_\theta\) could, for instance, be a fully connected neural network layer. We instead use LSTM cells [40] for \(f_\theta\) which allow the model to better learn long-term dependencies. In any case, given an initial state \(h_0\) and the parameters \(\theta\), an RNN represents a deterministic map from one sequence \(s_{1:n}\) to another sequence \(h_{1:n}\).

From the sequence \(h_{1:n}\) we can obtain a stochastic representation of the output sequence \(x_{1:n}\). We decompose \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n})\) as a product of conditional probabilities \[\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n}) = \prod^n_{i=1} \mathop{\mathrm{\mathrm{P}}}(x_i|x_{1:i-1}, s_{1:i}) \,.\] Specifically, we extend 7 by adding a sampling step to obtain \(x_i\) from \(s_{i}\) and \(x_{i-1}\) \[\begin{align} h_i &= f_\theta\left( \begin{bmatrix} s_i \\ x_{i-1} \end{bmatrix} , h_{i-1} \right) \label{eq:hi} \\ x_i \mid h_i &\sim \mathcal{N}(\mu(h_i), \sigma(h_i)) \end{align}\tag{8}\] such that each \(x_i\) is a normally distributed random variable with mean \(\mu(h_i)\) and standard deviation \(\sigma(h_i)\). \(\mu\) and \(\sigma\) are typically modeled as neural networks.

The autoregressive model represents a generative model for the distribution \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n})\). Generating \(x_{1:n} \mid s_{1:n}\) is done sequentially. For \(i=1,\ldots,n\), the sampling procedure alternates between computing \(h_i\) using 8 to get the parameters of the conditional distribution, \(\hat{\mu}_i = \mu(h_i)\), \(\hat{\sigma}_i = \sigma(h_i)\), and sampling the next \(x_i \sim \mathcal{N}(\hat{\mu}_i, \hat{\sigma}_i)\). Thus, the conditional probability of the resulting output sequence is given by \[\begin{align} \mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n}) &= \prod^n_{i=1} \mathop{\mathrm{\mathrm{P}}}(x_i|x_{1:n-1}, s_{1:n})\\ &= \prod^n_{i=1} \frac{1}{\sqrt{2\pi}\hat{\sigma}_i} \exp\left(\frac{(x_i-\hat{\mu}_i)^2}{2\hat{\sigma}^2_i}\right)\,. \end{align} \label{eq:autoregressive95conditional95density}\tag{9}\] In practice, we use the following form for \(\mu\) and \(\sigma\) \[\begin{align} \mu(h) &= W_\mu h + b_\mu \\ \sigma(h) &= \exp (W_\sigma h + b_\sigma) \end{align}\] where \(\boldsymbol{\theta}=(\theta, W_\mu, b_\mu, W_\sigma, b_\sigma)\) are model weights that are trained.

There are a few practical considerations for efficiently training the model. Training is performed in iterations and, as typically done for training neural networks, the loss function in Eq. 5 is only computed for a subset of the training data, in mini-batches of size \(N_\text{batch}=25\), instead of the whole training set of size \(N=\num{1000}\). At the beginning of each iteration, the subset that is used is then randomly resampled.

Appendix B: Details of the First Study—To generate training data for the neural network, we combine a linear auto-regressive input, with a stochastic nonlinear output model. Specifically, we considered an input that evolves according to \(\mathrm{AR}(p)\) statistics: \[S_t = \sum^p_{j=1} \phi_j S_{t-j} + \xi_t \label{eq:AR1}\tag{10}\] where \(\xi_t\) are iid random variables from a unit Gaussian distribution, and the \(\phi_j\in[0, 1)\) are model parameters. The output \(X_t\) is governed by the equation \[X_t = \sigma(\gamma S_t) + \rho X_{t-1} + \vartheta \eta_t \label{eq:nonlinear95model}\tag{11}\] where \(\eta_t\) are iid Gaussian random numbers, \(\gamma\), \(\rho\) and \(\vartheta\) are positive real parameters, and \[\sigma(x) = \frac{1}{1+e^{-x}}\] is the logistic function. The gain \(\gamma\) effectively controls the strength of the nonlinearity; see Fig. [fig:ARtest]a. This process models a response that saturates as the input grows. In fact, \(\sigma(x)\) is equivalent to the Hill function commonly used in biochemistry to describe saturating enzyme kinetics [41].

For this case study, we trained our machine learning model with synthetic data generated according to 10 11 for various values of the gain, denoted by \(\gamma\). The other parameters are specified in 1. For each value of \(\gamma\), we created a distinct training set of \(N=1000\) pairs of time series \((s_{1:50}, x_{1:50})\) and trained one autoregressive model per training set. Once the models were trained, we estimated the mutual information for each of them using PWS, employing variational inference to perform the marginalization; see also Appendix D: PWS simulation details below.

For the generative model, we use a recurrent neural network (RNN). Using LSTM cells with a hidden size of 64 along with a dense layer that outputs two values, representing the mean \(\mu\) and log-variance \(\ln\sigma^2\) of a Gaussian distribution. For each time step \(i\), the model receives input signals \(s_i\) and \(x_{i-1}\) and predicts the next output value \(x_i\) by sampling from this Gaussian. The model is trained by iteratively optimizing for 100 epochs over the training set in mini-batches of size 64 which are shuffled after each epoch. We optimize the model using the Adam optimizer [42] (\(b_1 = \num{0.9}, b_2 = \num{0.999}\)) with weight decay regularization [43] of \(\lambda=\num{1e-4}\) and using a cosine decay learning rate schedule [44] that smoothly decreases the learning rate from 1 × 10⁻² to 1 × 10⁻³ throughout the training process.

The inference model used for marginalization is also modeled as a LSTM and represents a variational approximation of the posterior \(\mathop{\mathrm{\mathrm{P}}}(s_{1:n}|x_{1:n})\). It can be seen as an autoregressive normalizing flow that transforms a sequence of iid Gaussian random numbers \(\epsilon_{1:n}\) into the samples \(s_{1:n}\) [45]. At each time step \(i\) it receives the current output observation \(x_i\) and a noise sample \(\epsilon_i \sim \mathcal{N}(0, 1)\). First, \(x_i\) is transformed by a fully connected neural network layer into a feature vector \(y_i\in\mathbb{R}^k\) where we choose \(k=16\). This vector \(y_i\) is concatenated with the previous prediction \(s_{i-1}\) and passed to an LSTM cell which outputs two parameters \(u_i\) and \(v_i\). These parameters parametrize the affine transformation \(s_i=u_i \epsilon_{i} + v_i\). The network is trained for 200 epochs with mini-batches of size 64. For each \(x_{1:n}\) in the training set, 16 Monte Carlo draws from the inference network \(\tilde{s}_{1:n}\sim q(\tilde{s}_{1:n}|x_{1:n})\) are used to estimate the ELBO loss. The loss function gradient is estimated using SGVB [37]. The model is optimized using the ADAM optimizer with weight decay regularization (same parameters as above). We use an learning rate of 1 × 10⁻².

Appendix C: Details of the Second Study—We analyze neuronal responses to a moving bar stimulus, recorded from retinal ganglion cells of a salamander. The dataset consists of spike trains from 230 neurons, sampled at a temporal resolution of 0.1 ms. Recordings span a total duration of 8148 s (approximately 2 hours) and are structured into 136 repeated trials, each lasting \(T = \SI{30.08}{\second}\). In each trial, the same stimulus is presented.

The stimulus is a vertically moving horizontal bar displayed on a screen. The \(x\)-position of the bar is varied stochastically. Its dynamics are given by a stochastically driven harmonic oscillator model. We rescale the input data, such that the variance \(\langle x^2 \rangle\) is 1. This model is defined by the following equations \[\begin{align} \tau \dot{x} &= v \\ \tau \dot{v} &= - \frac{1}{4 \gamma^2} x - v + \frac{1}{\sqrt{2} \gamma} \eta(t) \end{align}\] where the time-scale is given by \(\tau = \SI{50}{\milli\second}\) and the damping coefficient is given by \(\gamma = (2 \omega_0 \tau)^{-1} = \num{1.06}\) for a value of \(\omega_0 = \SI{9.42}{\per\second}\). Since \(\gamma > 1\), the oscillator is in the over-damped regime (albeit close to critically damped).

For analysis, we discretize the data at a resolution of \(\delta t = \SI{20}{\milli\second}\), yielding \(n = T/\delta t = \num{1504}\) time bins per trial. Sampling times are denoted as \(t_i = i\delta t\) for \(i = 0, \dots, n\). The spike count of neuron \(k\) in the interval \([t_i, t_{i+1})\) is represented as \(x^k_i\).

Our ML network models neuron spike probabilities as follows. In each time-interval \([t_n, t_{n+1})\) we assume the spike counts of each neuron are Poisson distributed. The firing intensity of neuron \(k\) at time \(t_i = i \Delta t\) is denoted by \(\lambda^k_i\). The random variable \(X^k_i\) denotes the number of spikes of neuron \(k\) in the time-interval \(t_{i:i+1} = [t_i, t_{i+1})\) and is assumed to be Poisson-distributed. The negative log likelihood is given by \[\mathcal{L} = \sum_i - n_i \ln \lambda_i \Delta t + \lambda_i \Delta t + \ln(n_i !)\]

We use a one-dimensional convolutional neural network for predicting the spiking activity of retinal cells in response to a visual stimulus. The network processes the time-varying stimulus \(s(t)\), which represents the visual input over time, using a series of causal convolutional layers. Each convolutional layer applies a 1D convolution with a specified kernel size, followed by batch normalization and a ReLU activation, allowing the network to extract temporal features from the input. To enforce causality—ensuring that predictions at a given time step only depend on past inputs—a zero-padding layer is added to the data before the convolutions, offsetting the receptive field appropriately. The final layer is a 1x1 convolution that maps the extracted features to spike intensity predictions for each neuron by modeling spiking activity as a Poisson process, where the predicted intensities determine the firing probability of each neuron.

Optimization is performed using the Adam optimizer, and the learning rate is dynamically adjusted using a OneCycleLR scheduler [46], which initially increases the learning rate before gradually decaying it to improve convergence.

Appendix D: PWS simulation details—The central idea of ML-PWS is to develop a generative model that can learn the conditional probability distribution \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n}|s_{1:n})\) and then compute the marginal output \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n})\) via Monte Carlo averaging in trajectory space. We find that for relatively short signal trajectories, as in the application of ML-PWS to neuronal data (Fig. [fig:NR]), it is sufficient to generate input trajectories \(s_{1:n}\) from an importance sampling distribution \(q(s_{1:n}|x_{1:n})\) using the ELBO variational procedure, and then obtain \(\mathop{\mathrm{\mathrm{P}}}(x_{1:n})\) by directly averaging over the full trajectories using Eq. 6 . In contrast, for longer signal trajectories, it is necessary to improve the sampling distribution \(q(s_{1:n}|x_{1:n})\) via a sequential Monte Carlo approach, generating the trajectories and performing the concomitant marginalization segment by segment, as in RR-PWS [31]. Indeed, for the application of ML-PWS to the autoregressive model (Fig. [fig:ARtest]), we combined RR-PWS with ELBO.

References↩︎