Designing optimal perturbation inputs for system identification in neuroscience
- Graduate School of Arts and Sciences, The University of Tokyo, Japan
- RIKEN Center for Brain Science, RIKEN, Japan
- Graduate School of Data Science, Yokohama City University, Japan
Figures
Figure 1
Passive and perturbation states and the resulting neural dynamics models.
(a) Passive state of the neural network and the corresponding ground-truth model. (b) Recorded neural activities without perturbation. (c) Estimated model obtained without perturbation. (d) Perturbation state of the neural network and corresponding ground-truth model. (e) Recorded neural activity with perturbation. (f) Estimated model obtained with perturbation.
Figure 2
Visualization of continuous-time system dynamics by perturbations.
(a) System matrix representing the dynamics and its eigenvalues. (b) Temporal changes in the state variables of the passive state without perturbation. (c) Time-varying trajectories of each PCA component. The color gradient indicates temporal progression. (d) Estimated connectivity matrix in the passive condition. (e) Error matrix for the passive-state estimation. (f–i) Temporal changes and estimation results in response to an impulse input. (j–m) Temporal changes and estimation results in response to a sinusoidal input.
Figure 3
Effect of perturbation on the eigenvalues of the state covariance matrix.
Ellipsoidal plots show the eigenvalues and reciprocal eigenvalues of the covariance matrix of the state trajectories , comparing passive dynamics and the case with perturbation. All ellipsoidal plots (c–f) are drawn along the first three eigenvector directions of ; panels (c, e) scale the axes by the eigenvalues , while panels (d, f) scale them by the reciprocals . (a) Examples of a passive signal, a perturbation input, and a perturbed signal. (b) Analytical relationship linking the eigenvalues of the state covariance matrix to system identification error. (c) Passive: The covariance matrix is dominated by a single eigenvalue direction. (d) Passive (reciprocal): A single dominant eigenvalue produces a widely spread reciprocal-space ellipsoid (e) Perturbation: Additional dynamical modes are excited, increasing the smaller eigenvalues. (f) Perturbation (reciprocal): Perturbations yield a more uniform reciprocal-space ellipsoid, reducing estimation error.
Figure 4
Error in the system matrix and eigenvalue-scaled ellipsoids.
(a) State transition matrix and its eigenvalues . (b) System frequency responses, showing the sum of eigenvalues of the state covariance matrix (left), the sum of inverse eigenvalues of the state covariance matrix (middle), and the system matrix error (right). (c) Eigenvalue ellipsoids of under passive, 5 , 10 , and 15 Hz conditions. (d) Reciprocal-eigenvalue ellipsoids of under passive, 5 , 10 , and 15 Hz conditions.
Figure 5
Simulation of a multi-mode system and eigenvalue-scaled ellipsoids.
(a) System matrix and its eigenvalues; the system exhibits two damped oscillatory modes at 10 and 20 Hz. (b) Theoretical estimation error, defined as for two-frequency input combinations; cooler colors indicate smaller estimation errors. (c) Comparison of the reciprocal eigenvalue sum for flat-spectrum versus composite-frequency inputs. (d) Eigenvalue-scaled ellipsoids (μ-scaled) under four input conditions (columns): Passive, 10 , 20 , and 10 Hz +20 Hz; ellipsoid axes align with the eigenvectors. (e) Relative-scale view of (d). (f) Eigenvalue-scaled ellipsoids (inverse-scaled, ) for the same four conditions. (g) Relative-scale view of (f).
Figure 6
Optimal stimulation location (a hub node or the heavily-damped subnetwork) minimizes estimation error.
A 7-node network with three oscillatory modes (15, 25, and 35 Hz) is used to identify effective perturbation locations. The 15 Hz mode (Nodes 1–2) is heavily damped; Node 7 is a hub with only outgoing edges. (a) The 7×7 connectivity matrix . (b) Graph representation. (c) Damping rates of the eigenvalues of . (d–f) 3D ellipsoids whose axes are proportional to , , (reciprocals of the three smallest eigenvalues of ; larger axes indicate greater estimation error): (d) passive, (e) impulse to Node 1, (f) impulse to Node 7. (g) Heatmap of across input nodes (P=passive); uniformly small values appear only for Node 7. (h) Full-matrix estimation error per input node. (i) Estimation error for the submatrix excluding Node 7; Nodes 1, 2, and 7 achieve comparable errors, confirming that direct stimulation of the damped subnetwork matches hub stimulation.
Figure 7
Neural state classification enhanced by perturbation-aided model estimation.
(a) Ground-truth dynamics models defining five distinct task conditions, each with unique eigenvalue distributions and network structures of matrix . The optimally excited node is highlighted in red. (b) Simulated neural states (activity) induced by the five task conditions, shown separately for the passive (top) and perturbed (bottom) regimes. (c) Estimated model matrices visualized using linear discriminant analysis (LDA). (d) Confusion matrices showing classification accuracy for the passive condition (left) and the optimally designed perturbation condition (right). (e) ROC curves for the passive condition (left) and a comparison of random and optimal perturbation conditions (right). (f) LDA plots for passive condition with increasing time window length, showing clustering and classification accuracy as the time window length increases from 20 to 200.
Figure 8
Effects of system identification on network control theory.
(a) State transitions using models identified from the passive and perturbed conditions. (b) Controllability Gramians with two models. (c) Control costs computed from the estimated controllability Gramians.
Figure 9
Accurate system identification via perturbation enables precise optimal control.
(a) True matrix and its network structure. Nodes 1, 2, and 5 are controllable nodes. Nodes 3 and 4 are targets controlled by the determined control input. When the connectivity matrix is misestimated, the structure has edges between nodes 1 and 2, and nodes 3 and 4. (b) Matrix estimation failure under passive conditions and controlled state transitions. The control objectives are plotted as white circles. (c) Successful matrix estimations and controlled state transitions with sufficiently strong impulse inputs (). (d) Errors in the estimated matrices and and the control cost (average controllability).
Figure 10
Iterative perturbation design converges to the theoretical optimum in a 32-channel neural network simulation.
(a) True matrix and its network structure, which was designed to include both oscillatory and attenuating components. (b) Comparison of estimation errors. The purple line indicates the iterative design process initialized from the passive-state estimate (Iteration 1). The blue line shows the iterative design starting from a flat-spectrum input (Iteration 1). Both processes converge towards the theoretically optimal error (dashed line) as the design is refined (design-1, design-2). The third line (labeled ‘Random’) represents a baseline in which a single randomly chosen node is stimulated with a flat-spectrum input throughout all iterations. (c) Schematic of the iterative procedure. An initial estimated system matrix () is used to design the first perturbation input (), which yields a better estimate (), and the cycle repeats.
Appendix 2—figure 1
Robustness of the perturbation framework to nonlinear state dynamics.
(a) The nonlinearity with the min–max envelope of within the ordinary least squares (OLS) window shaded for . (b) Relative 1-step prediction RMSE vs. β. (c) Relative free-run (rollout) RMSE vs. β. (d) Training trajectories within the OLS estimation window for , under passive observation and impulse perturbation. (e) The same trajectories shown in full; the gold shade marks the OLS window and the dotted line marks the impulse time. (f) Held-out test trajectories (black) overlaid with the per-segment 1-step prediction from . (g) The same held-out test trajectories overlaid with the free-run (rollout) prediction. (h) The OLS-window data in the top three principal components (independent principal component analysis (PCA) per cell; viridis = sample index).
Appendix 2—figure 2
Perturbation framework can still improve the estimation of A even when B is unknown, provided that A and B are recovered jointly by ordinary least squares (OLS).
(a) Top row: which channel(s) each condition stimulates (filled nodes are active; gray nodes are silent). Middle and bottom rows: recovered and at from a single trial. Seven conditions are compared: ground truth, passive observation, a naive baseline that uses the all-nodes input but assumes , and joint estimation under four perturbation patterns (node 1 only, node 3 only, nodes 1+3, all nodes). (b) Corresponding sample state trajectories; the gold shade marks the OLS window and dotted lines mark the impulse times. (c) as a function of ε (log scale). (d) as a function of ε.
Appendix 2—figure 3
Variation of eigenvalues and system matrix error with impulse and step input strengths.
(a) Assumed system matrix and its eigenvalues. (b) Changes in the sum of eigenvalues, the sum of their reciprocals, and the estimation error of with respect to impulse input strength. (c) Changes in the same quantities with respect to step input strength.
Appendix 2—figure 4
Comparison between white noise input and optimized frequency input.
(a) Matrix A with 16 modes, assuming 32 channels. (b) Channel-averaged power spectrum in the passive state. (c) Comparison of the sum of reciprocal eigenvalues () of the covariance matrix of the state vector under white noise and optimized frequency input. (d) Comparison of eigenvalue-scaled ellipsoids obtained from white noise and optimized frequency input, visualizing the components corresponding to the 1st–3rd and 30th–32nd eigenvalues.
Tables
Appendix 2—table 1
Conditions for generating neural dynamics.
| Section | Sampling frequency | Time length | Data length T | Number of trials | Noise | Input strength |
|---|---|---|---|---|---|---|
| Figure 2 (Failure of passive estimation) | 1000 | 10 | 10,000 | 100 | 10−32 | 1 |
| Figures 4 and 5 (Sinusoidal inputs) | 100 | 50 | 5000 | 10,000 | 10−2 | 1 |
| Figure 6 (Location tuning) | 1000 | 1 | 1000 | 200 | 10−2 | 10 |
| Figure 7 (Neural state classification) | 10 | 2 | 20 | 100 | 10−2 | 102 |
| Figure 9 (Optimal control) | 50 | 1 | 50 | 100 | 10−2 | 103 |
| Figure 10 (Practical designing procedure) | 100 | 25 | 2500 | 500 | 10−2 | 0.1 |
| Appendix 2—figure 1 (Nonlinear dynamics) | 1000 | 1 | 1000 | 30 | 10−32 | 0.15–0.3 |
| Appendix 2—figure 2(B-Matrix uncertainty) | 1000 | 1 | 1000 | 30 | 10−32 | 0.15–0.3 |
| Appendix 2—figure 3 (Impulse/step) | 100 | 50 | 5000 | 1 | 10−2 | 0.5–5 |
| Appendix 2—figure 4 (White noise/designed) | 100 | 100 | 10,000 | 1 | 10−2 | 0.1 |
Additional files
Download links
A two-part list of links to download the article, or parts of the article, in various formats.
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Designing optimal perturbation inputs for system identification in neuroscience
eLife 15:RP110030.
https://doi.org/10.7554/eLife.110030.3
