Designing optimal perturbation inputs for system identification in neuroscience

eLife Assessment

This important study introduces a method, based on active control, for determining the functional connectivity between neurons or groups of neurons. While the theory is derived for a linear model with a perfectly observed control signal, the authors present compelling evidence that their method will work even when the system is nonlinear and the control signal is not fully observed. Furthermore, they provide a practical recipe for doing so.

https://doi.org/10.7554/eLife.110030.3.sa0

Significance of the findings:

Important: Findings that have theoretical or practical implications beyond a single subfield

Landmark
Fundamental
Important
Valuable
Useful

Strength of evidence:

Compelling: Evidence that features methods, data and analyses more rigorous than the current state-of-the-art

Exceptional
Compelling
Convincing
Solid
Incomplete
Inadequate

During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments

Abstract
Introduction
Results
Discussion
Methods
Appendix 1
Appendix 2
Data availability
References
Article and author information
Metrics

Abstract

Investigating the dynamics of neural networks, which are governed by connectivity between neurons, is a fundamental challenge in neuroscience. Because passive (spontaneous) activity provides only limited information for estimating connectivity, perturbation-based approaches are widely applied in neuroscience as they can evoke underlying hidden dynamics. However, the characteristics of such perturbations have typically been designed based on empirical or biological intuition. To enable more accurate estimation of connectivity, we propose a data-driven and theoretically grounded framework for optimally designing perturbation inputs, based on formulating the neural model as a control system. The core theoretical insight underlying our approach is that neural signals observed in the passive state lack sufficient latent information, which leads to failures in the system identification. Perturbations reveal these hidden dynamics and lead to improved estimation. Guided by these insights, we derive a theoretical basis for optimizing perturbation inputs that minimize estimation errors in neural system identification. Building upon this, we further explore the relationship of this theory with stimulation patterns commonly used in neuroscience, such as frequency, impulse, and step inputs. We demonstrate the effectiveness of this framework for neuroscience through simulations grounded in experimental paradigms, such as neural state classification and optimal control of neural states. Our theoretical analysis, together with multiple simulations, consistently shows that perturbations designed according to our framework achieve substantially more accurate system identification compared to the conventional, intuition-based inputs. This study provides a theoretical foundation for designing perturbation inputs to achieve accurate estimation of neural dynamics. This, in turn, enables reliable discrimination of neural states such as levels of consciousness and pathological conditions, and facilitates precise control of their transitions toward recovery from abnormal states.

Introduction

Much recent interest has focused on how interactions between individual neurons and neuronal populations support cognitive functions and behavior, and on how the disruption of these interactions contributes to various neurological and psychiatric disorders (Waites et al., 2006; Friston, 2011; Boly et al., 2012; Sporns, 2014; Petersen and Sporns, 2015). This interaction—called neuronal connectivity (Sporns, 2016; Seguin et al., 2023)—is often represented as a network structure or a connectivity matrix. In particular, a commonly studied aspect is functional connectivity, which captures the temporal dependencies between neural activities (Stevens, 2009; Honey et al., 2009; Rubinov and Sporns, 2010; Friston et al., 2013). The functional connectivity is estimated through recording techniques, such as neural spike recording techniques (Siegle et al., 2021; Yan and Murphy, 2024), electrocorticography (ECoG) (Park et al., 2016; Kucyi et al., 2018), functional magnetic resonance imaging (fMRI) (Greicius et al., 2007; Rocchi et al., 2022) and electroencephalography (EEG) (Nentwich et al., 2020; Bogéa Ribeiro and da Silva Filho, 2023). While various methods have been proposed and used to estimate these connections (see for example Friston, 2011; Reid et al., 2019 for a comprehensive review), one typical method is based on dynamic modeling, among which is a simple but widely used linear auto-regressive model (Figure 1a; Goebel et al., 2003; Ding et al., 2006; Ting et al., 2015; Seth et al., 2015; Cekic et al., 2018). Connectivity is statistically estimated from the time series data of neural activity (Figure 1b) by fitting the model parameters (Figure 1c). A neural connectivity matrix helps visualize the connections between different brain regions, and provides a detailed map of functional interactions within the brain. By comparing connectivity matrices across individuals or groups, researchers can deepen their understanding of neuronal architectures and communication across various disciplines (Stevens, 2009; Ruiz et al., 2014; Bassett and Sporns, 2017; Lee et al., 2020; Bergmann et al., 2021; Siddiqi et al., 2022).

Figure 1

Download asset Open asset

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.

A fundamental problem with passive observation is that it fails to reveal hidden dynamics, which leads to an invalid estimate of the corresponding parts of the model. This limitation stems from the attenuation of latent dynamical modes, such as transient or damped components. When neural activity is modeled as a linear dynamical system, as shown in Figure 1a, the true connectivity matrix governs the generation of time-series signals (Figure 1b). Passive observation records these spontaneous neural activities within an arbitrarily predefined temporal window and attempts to estimate the connectivity matrix using parameter estimation techniques. However, the resulting matrix $A_{passive}$ deviates significantly from the true model (Figure 1c), because some modes quickly decay and vanish from the observable data (see also Figure 2 for an intuitive explanation). Consequently, relying solely on passive recordings leads to misinterpretations of the neural system—such as overlooking fast transient dynamics or underestimating connectivity strength—which, in turn, may result in inaccurate conclusions in both basic neuroscience and clinical contexts.

Figure 2

Download asset Open asset

Visualization of continuous-time system dynamics by perturbations.

(a) System matrix $A$ 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.

To overcome these limitations, the application of perturbations has emerged as a powerful approach for improving the estimation of connectivity in neuroscience (Bergmann et al., 2021; Casali et al., 2013; Stepaniants et al., 2020; Lepperød et al., 2023). Its effectiveness lies in its ability to actively manipulate the neural system, rather than merely observing it, and to thereby uncover dynamic and causal interactions that are otherwise hidden (Stepaniants et al., 2020). Recent studies using stimulation techniques, such as optogenetics (Lepperød et al., 2023; Steinberg and Janak, 2013) and transcranial magnetic stimulation (TMS) Bergmann et al., 2021; Hallett et al., 2017; Tik et al., 2023 have demonstrated that these interventions can significantly improve the detection of causal relationships within the neural system. Such approaches have been applied to the study of cognitive processes and states of consciousness (Bergmann et al., 2021; Casali et al., 2013). These studies underscore the value of perturbation inputs in driving state transitions in the neural system, thereby enabling more accurate and reliable connectivity estimation.

However, the parameters of perturbation protocols, such as the location, intensity, and shape of stimulation, have often been determined based on anatomical and physiological insights and on empirically established techniques (Bergmann et al., 2021; Massimini et al., 2007; George, 2019; Bernal-Casas et al., 2017; Grimm et al., 2024). While these approaches have provided practical utility, they may not optimally exploit the underlying neural dynamics or maximize the informativeness of the perturbation. To overcome these limitations and enhance the reliability of connectivity inference, there is a pressing need to design data-driven and theoretically grounded perturbations.

In this paper, we propose a framework for designing the optimal perturbation input through control theory in neuroscience. We interpret neural dynamics as a control system (Gu et al., 2015; Deng et al., 2022; Kawakita et al., 2022; Kamiya et al., 2023; Moradi Amani et al., 2024; Shikauchi et al., 2025), and treat external perturbations as control inputs to design properties of neural stimulation (Figure 1d). If the optimal perturbation input can be systematically designed, it becomes possible to steer the neural system toward states that are maximally informative (Figure 1e), thereby enhancing the accuracy of the inferred connectivity (Figure 1f). While recent studies have begun to develop algorithmic approaches to active stimulation design for specific experimental platforms (Minai et al., 2024; Wagenmaker et al., 2024), a general theoretical framework for why and which perturbation inputs are effective has yet to be established. We first describe how to formulate neural dynamics as a control system and how to estimate the model parameters from observed data. Building upon this formulation, we derive a theoretical basis that enables us to design the optimal perturbation inputs for the neural system identification. We demonstrate the validity and utility of this theoretical basis by exploring its implications for optimizing parameters of common neurostimulation techniques and by applying it to practical examples, including neural state classification (Casali et al., 2013; Hallett et al., 2017; Shikauchi et al., 2025; Angulo-Sherman et al., 2017) and control of neural states (Gu et al., 2015; Kawakita et al., 2022; Kamiya et al., 2023; Karrer et al., 2020). In these demonstrations, we define concrete problems and apply the theory to validate its practical utility.

Our research offers a comprehensive framework for system identification with perturbation in neuroscience, and paves the way for more precise and effective analysis of neural dynamics. By providing clear guidelines on the design and application of perturbations, this framework serves as a practical reference for experimental settings, helping researchers determine the most effective stimulation parameters for their studies.

Background

To evaluate how external perturbations enhance the accuracy of system identification, this section reviews the estimation formulations for linear dynamical systems, comparing the cases with and without perturbation inputs.

System identification with external perturbation

This section reviews a well-established framework for system identification of linear dynamical systems with external perturbations. Researchers have modeled brain dynamics as a discrete-time linear dynamical system with external inputs and stochastic system noise (Gu et al., 2015; Deng et al., 2022; Betzel et al., 2016; Kim et al., 2018; Braun et al., 2021). We assume that the brain state evolves according to the following linear dynamics:

x (t + 1) = A x (t) + B u (t) + ξ (t),

where $x (t) \in R^{n}$ denotes the neural state vector at time t, such as the activity of neurons, populations, or brain regions, and n represents its dimension. The connectivity matrix $A \in R^{n \times n}$ characterizes intrinsic interactions reflecting functional connectivity. The input matrix $B \in R^{n \times m}$ specifies how external inputs $u (t) \in R^{m}$ —including sensory stimuli or interventions, such as TMS—affect the system, where $m$ represents the number of perturbation channels. The input $u (t)$ is freely designed and known to the experimenter. Finally, $ξ (t) \in R^{n}$ represents stochastic fluctuations (e.g. synaptic noise or unobserved inputs), modeled as i.i.d. zero-mean Gaussian noise across time (i.e. delta-correlated in time) with covariance matrix $Σ_{ξ} = E [ξ (t) ξ (t)^{⊤}] ≻ 0$ . We remark that neural activities measured at the macroscopic level, such as functional magnetic resonance imaging (fMRI) or intracranial electroencephalography (iEEG) have been experimentally and theoretically validated to follow approximately linear dynamical systems in Ahmed and Nozari, 2022; Nozari et al., 2024.

Here, we assume the input matrix B is known, and we are only estimating the connectivity matrix A from the time series data of x and u. We construct the data matrices as

X = [x (0), x (1), \dots, x (T - 1)], Y = [x (1), \dots, x (T)], U = [u (0), \dots, u (T - 1)],

where T represents the length of the time series data. Using these matrices, the parameter matrix A can be estimated by minimizing the reconstruction error in the sense of ordinary least squares (OLS):

\hat{A} = \underset{A}{\arg min} ‖ Y - (A X + B U) ‖_{F}^{2}

= (Y - B U) X^{†},

where $‖ \cdot ‖_{F}$ represents the Frobenius norm and the symbol ^† denotes the pseudoinverse of a matrix. This manuscript derives a theoretical framework based on this formulation, which assumes linear dynamics with the input matrix B known. The relaxation of these assumptions—nonlinear state dynamics and the case of unknown B—is addressed in Appendix 2.2 and Appendix 2.3.

System identification without external perturbation

To compare the quality of system identification with and without perturbations, we also consider the passive case without external input. This corresponds to setting $u (t) = 0$ in the formulation in the perturbed case, yielding the following dynamics:

x_{passive} (t + 1) = A x_{passive} (t) + ξ (t) .

Here, A and $ξ (t)$ are identical to those defined in Equation 1. Because this system evolves differently from the perturbed case, we denote its state trajectory as $x_{passive} (t)$ to explicitly distinguish it from the dynamics under perturbation.

Similarly, the parameter matrix A in the passive condition can be estimated by applying the ordinary least squares (OLS) method:

{\hat{A}}_{passive} = Y_{passive} X_{passive}^{†},

where

X_{passive} = [\begin{matrix} x_{passive} (0), x_{passive} (1), \dots, x_{passive} (T - 1) \end{matrix}], Y_{passive} = [\begin{matrix} x_{passive} (1), \dots, x_{passive} (T) \end{matrix}] .

The subscript ‘passive’ is used again to explicitly distinguish variables associated with the unperturbed dynamics from those obtained under external perturbations.

Neural stimulation as control inputs

This section describes how commonly used neural stimulation techniques can be related to input signals in control theory. Their adjustable parameters vary depending on how the stimulation inputs are modulated.

Three non-invasive electrical stimulation methods illustrate how stimulation paradigms map onto basic control inputs. Transcranial magnetic stimulation (TMS) induces brief and transient perturbations via electromagnetic pulses (Momi et al., 2023), which are naturally represented as a sequence of impulse-like inputs, where the timing and intensity of each pulse are the primary controllable parameters. Transcranial direct current stimulation (tDCS) primarily modulates neural activity through approximately constant inputs (Tapsell et al., 2025), which can be viewed as a step-like signal whose main controllable parameter is the amplitude of the applied current. Transcranial alternating current stimulation (tACS) delivers oscillatory inputs (Grover et al., 2023), corresponding to sinusoidal signals characterized by amplitude, frequency, and phase. In control theory, impulse, step, and sinusoidal inputs are the basic components used to characterize system responses and dynamics (Ogata, 2010; Nise, 2020).

The control input framework extends beyond non-invasive techniques to invasive and optogenetic stimulation. Invasive electrical stimulation, including intracranial microstimulation and deep brain stimulation (DBS), enables direct delivery of electrical inputs to neural tissue (Lozano et al., 2019), providing flexible control over amplitude and timing through pulse trains or temporally structured waveforms. Optogenetic stimulation allows genetically targeted activation or inhibition of specific neurons using light (Deisseroth, 2011), providing fine-grained control over multiple input dimensions, including amplitude (light intensity), temporal pattern, and cell-type specificity. In particular, recent developments enable stimulation at the level of individual neurons with high temporal precision (Tong et al., 2023; LaFosse et al., 2023), allowing flexible construction of spatiotemporal input patterns.

These stimulation examples demonstrate that the theoretical framework developed in this paper connects to practical experimental settings. While a substantial gap remains between idealized control inputs in theory and experimentally realizable stimulation, the core principles established in the following sections provide a foundation that naturally extends to these practical stimulation paradigms.

Results

We present a theoretical framework for determining the optimal perturbation inputs to enhance neural system identification. We begin by presenting an intuitive example in which the failure of passive estimation is demonstrated, and the benefit of perturbation inputs in re-exciting weakly observable dynamics becomes evident. Following this, we provide a guiding theoretical principle that estimation error of A is reduced by excitation of the state $x$ by perturbation input u. We then validate this theoretical foundation by applying it to canonical perturbation signals used in neuroscience—sinusoidal, impulse, and step inputs, as typically delivered by tACS, TMS, and tDCS, respectively—and derive analytical relationships between input parameters and estimation errors. Furthermore, we demonstrate the applicability of the proposed framework to practical problems in neuroscience, such as neural state classification and optimal control for neural state transitions. Finally, we outline a practical framework for designing optimal perturbation inputs to estimate neural dynamics.

Typical failure of passive-state model estimation

We demonstrate why passive state observation, a common experimental condition in system neuroscience, fails to estimate the system model of neural dynamics. Although passive observation is widely employed due to its convenience and non-invasiveness, this practice has fundamental limitations: it inevitably overlooks causal and dynamical information that vanishes under spontaneous conditions but can be recovered through external perturbations, leading to degraded connectivity estimates. Such limitations lead to misinterpretations of the underlying neural functions when neural dynamics are modeled as a control system.

We estimate the connectivity matrix from simulated time series data generated by the true model. The true model, characterized by a connectivity matrix A, is illustrated in Figure 2a. Neural data are generated based on the connectivity matrix, yielding the time series shown in Figure 2b. These data are not influenced by external input (passive state), and thus the oscillations of $x_{3}$ and $x_{4}$ gradually attenuate over time. This attenuation is further confirmed by principal component analysis (PCA), which reveals that the first and second principal components capture the 10 Hz dynamics, while the third principal component, associated with the 40 Hz mode, contributes negligibly to the total variance (Figure 2c). Consequently, the matrix estimated by OLS deviates substantially from the true connectivity matrix (Figure 2d), particularly showing errors in the lower-right components, as highlighted in Figure 2e. This failure arises because the true matrix A contains two distinct dynamical modes: one is a continuous oscillatory mode at 10 Hz, and the other is a damped mode at 40 Hz. The damped mode becomes unobservable in the time series once its contribution has attenuated and vanished. Passive recordings capture only superficial aspects of the connectivity and fail to estimate the true underlying neural dynamics.

The limitation of the passive state can be resolved by applying perturbation inputs. Figure 2f–i shows the application of an impulse input u to the system. This perturbation input primarily affects the nodes corresponding to $x_{3}$ and $x_{4}$ , reintroducing variations in the third principal component direction. Estimating matrix A from this perturbed time-series data yields a more accurate estimate than in the passive case. Similarly, an improvement is observed when different types of input—a sinusoidal input, for example—are used instead of the impulse input, as shown in Figure 2j–m. These results illustrate that appropriate perturbation inputs can effectively re-excite weak dynamics, thereby enhancing system identification accuracy. This can be theoretically supported by the concept of persistent excitation (Green and Moore, 1986; Shimkin and Feuer, 1987; Jenkins et al., 2018; Lu and Cannon, 2023). According to this theory, a perturbation input is said to be persistently exciting if it causes the system’s state to sufficiently explore the state space over time. These insights underscore the importance of incorporating controlled stimulation in experimental design, particularly when accurate system identification is desired. A similar trend under nonlinear dynamics is shown in Appendix 2.2.

Intuitive perturbation design informed by covariance eigenvalues

In this section, we show how perturbations reduce the estimation error of the system matrix A by exciting the state dynamics, thereby providing a theoretical basis for designing effective perturbation inputs. When the data length T is sufficiently large ( $T \to \infty$ ), the estimation error of A asymptotically converges to the following expression (Hamilton, 1994):

\begin{array}{ll} E [‖ \hat{A} - A ‖_{F}^{2}] \approx \frac{1}{T} t r (Σ_{ξ}) t r (Σ_{X}^{- 1}) . \end{array}

where $Σ_{X} = \frac{1}{T} (X - \bar{X}) (X - \bar{X})^{⊤}$ with $\bar{X} = \frac{1}{T} \sum_{t = 0}^{T - 1} x (t)$ represents the covariance matrix of the state vector x, and tr denotes the trace of a matrix. The detailed derivation of this equation is provided in Appendix 1.1. This equation indicates that the estimation error of A can be reduced either by increasing the measurement duration T or by minimizing the trace of the inverse covariance matrix $Σ_{X}^{- 1}$ , which depends on the perturbation input $u (t)$ . Intuitively, this relationship resembles a signal-to-noise ratio, where increasing the covariance of the state dynamics relative to the noise covariance leads to improved estimation accuracy.

The asymptotic expression in Equation 8 reveals that the estimation error of A depends inversely on the eigenvalues of the covariance matrix $Σ_{X}$ . Expressing this relationship explicitly in terms of the eigenvalues yields

\begin{array}{ll} E [‖ \hat{A} - A ‖_{F}^{2}] \approx \frac{1}{T} t r (Σ_{ξ}) \sum_{i = 1}^{n} \frac{1}{μ_{i}} . \end{array}

where $μ_{i}$ denotes the i-th eigenvalue of $Σ_{X}$ . This relationship indicates that small eigenvalues dominate the estimation error through their inverse contributions. Therefore, the perturbation input should be designed not only to increase the overall variance of the state x, but to enlarge the eigenvalues of $Σ_{X}$ more uniformly across all directions of the state space. In other words, exciting all dynamical modes of the system—rather than amplifying a limited subset—is essential for improving estimation accuracy.

This theoretical insight can be illustrated intuitively in Figure 3. The figure visualizes ellipsoids constructed from the eigenvalues $μ_{i}$ and their reciprocals $1 / μ_{i}$ of the covariance matrix $Σ_{X}$ . The covariance matrices $Σ_{X}$ are calculated from the time series of the state vector x obtained under passive and perturbation conditions (Figure 3a). The length of each axis of the ellipsoid corresponds to the variance of the state along that direction, which is associated with the eigenvalue $μ_{i}$ . Its reciprocal counterpart $1 / μ_{i}$ represents the contribution to the estimation error (Figure 3b). As shown in Figure 3c, when the neural dynamics is dominated by a single eigenvalue $μ_{1}$ , the other eigenvalues $μ_{2}$ and $μ_{3}$ remain small, resulting in an elongated reciprocal ellipsoid (Figure 3d) and a large estimation error. When perturbation inputs that excite the directions associated with $μ_{2}$ and $μ_{3}$ are applied, the ellipsoid expands and becomes closer to a sphere (Figure 3e). As a result, the estimation error decreases in all directions, as illustrated in Figure 3f. This uniform enlargement of the eigenvalues provides a clear guideline: perturbations should be designed so that the variance becomes large in all directions, rather than being confined to specific modes.

Figure 3

Download asset Open asset

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 $X$ , comparing passive dynamics and the case with perturbation. All ellipsoidal plots (c–f) are drawn along the first three eigenvector directions of $Σ_{X}$ ; panels (**c, e**) scale the axes by the eigenvalues $μ_{i}$ , while panels (**d, f**) scale them by the reciprocals $1 / μ_{i}$ . (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.

Theoretical formulation of sinusoidal inputs for an overall increase in eigenvalues

Based on the guideline presented in the previous section, we derive analytical expressions to investigate which frequency of perturbation input will efficiently excite the dynamical modes of a system. In neuroscience, such a perturbation is analogous to tACS. Since neural dynamics possess modal frequencies, there should exist perturbation input frequencies that resonate with these modes. Such resonance leads to an amplification of $x (t)$ , which generally results in larger eigenvalues $μ_{i}$ of the covariance matrix $Σ_{X}$ , reflecting increased variability along the corresponding modes. To clarify this relationship, we derive an analytical expression for $x (t)$ to investigate how the input frequency affects its amplitude. The solution $x (t)$ to the model in Equation 1 can be obtained by iterating the state transition matrix. Specifically,

x (t) = A^{t} x (0) + \sum_{k = 0}^{t - 1} A^{t - 1 - k} B u (k) + \sum_{k = 0}^{t - 1} A^{t - 1 - k} ξ (k)

We assume a cosine input $u (k) = \sum_{l = 1}^{L} \cos (ω_{l} k) u_{0}$ , where $ω_{l}$ are the angular frequencies. By substituting this into the second term, we can obtain the following equation by performing diagonalization of $A$ and expressing the eigenvalues in polar form.

x (t) = x_{passive} (t) + x_{diff} (t),

x_{passive} (t) = A^{t} x (0) + \sum_{k = 0}^{t - 1} A^{t - 1 - k} ξ (k)

x_{diff} (t) = \sum_{l = 1}^{L} \sum_{d = 1}^{n} Re [\frac{r_{d}^{t + 1} e^{i θ_{d} (t - 1)} - r_{d}^{t} e^{i (θ_{d} t - ω_{l})} - r_{d} e^{i (ω_{l} t - θ_{d})} + e^{i ω_{l} (t - 1)}}{r_{d}^{2} - 2 r_{d} \cos (θ_{d} - ω_{l}) + 1} v_{d} w_{d}^{⊤} B u_{0}]

where each eigenvalue of $A$ is written in polar form as $λ_{d} = r_{d} e^{i θ_{d}}$ , with $r_{d} = | λ_{d} |$ denoting the magnitude and $θ_{d} = \arg (λ_{d})$ denoting the argument. $x_{passive}$ represents the state x that would be realized if no perturbation were applied, i.e., the state that would evolve according to Equation 5. The vector $v_{d}$ denotes the right eigenvector of A associated with $λ_{d}$ , and $w_{d}^{⊤}$ denotes the corresponding left eigenvector. The detailed derivation is provided in Appendix 1.2.1.

We focus on $x_{diff} (t)$ simply because it is the component directly driven by the control input. When the input frequencies $ω_{l}$ coincide with the mode’s angular component $θ_{d}$ , a resonance-like amplification occurs in $x_{diff} (t)$ —that is, the amplitude of $x (t)$ increases markedly. Increasing the magnitude of $x (t)$ through resonance naturally leads to an increase in the overall variance in $Σ_{X}$ . This enhancement directly contributes to increasing all eigenvalues $μ_{i}$ of the covariance matrix rather than leaving some unexcited. Therefore, the analysis of Equation 13 provides the necessary guidelines to target and amplify specific dynamical modes, suggesting that oscillatory inputs, such as tACS, should be designed with frequencies that correspond to the true dynamical mode frequencies.

Theoretical formulation of impulse and step inputs

The eigenvalues of the covariance matrix can also be increased by manipulating the intensity of the perturbation input. In neuroscience, external perturbations, such as TMS and tDCS, can be modeled as impulse-like or step-like inputs to neural systems. The state vector x(t) for impulse inputs can be written as:

x_{passive} (t) = A^{t} x (0) + \sum_{k = 0}^{t - 1} A^{t - 1 - k} ξ (k), x_{diff} (t) = α A^{t - 1} B u_{0} .

For a step input, the state vector is described by,

x_{passive} (t) = A^{t} x (0) + \sum_{k = 0}^{t - 1} A^{t - 1 - k} ξ (k), x_{diff} (t) = β \sum_{k = 0}^{t - 1} A^{t - 1 - k} B u_{0} .

where α and β are the intensities of the impulse and step inputs. The detailed derivation is provided in Appendix 1.2.2. From these equations, it is clear that the contribution of the perturbation to the state covariance increases proportionally to the squares of the input strengths, $α^{2}$ and $β^{2}$ . However, estimation accuracy is not determined by input intensity alone. Equations 14; 15 show that the resulting dynamics $x_{diff} (t)$ are critically dependent on the term $B u_{0}$ . This term represents how the input vector $u_{0}$ (which defines the spatial pattern of the stimulation, i.e. which nodes are targeted) interacts with the system’s input matrix B. Therefore, while increasing intensity (α, β) within experimental constraints is beneficial, the spatial pattern $u_{0}$ is the key design parameter that determines which dynamical modes are excited. An improperly chosen $u_{0}$ may excite only the modes already dominant in the passive state, failing to enlarge the smaller eigenvalues that are critical for system identification. The problem of how to design the optimal spatial pattern $u_{0}$ to most effectively excite the weakly observable dynamics will be addressed in a later section.

Demonstration of sinusoidal inputs

In this section, we demonstrate how to design the frequency of oscillatory input by analyzing the eigenvalues of the covariance matrix. Theoretical formulations indicate that tuning the frequency of oscillatory perturbation inputs is more complex compared to that of impulse and step inputs. We generate neural dynamics governed by the connectivity matrix A, which is characterized by designated dynamical modes and compare its eigenvalues $λ_{i}$ with eigenvalues $μ_{i}$ of the covariance matrix of the perturbed state vectors. The demonstration for impulse and step inputs can be found in Appendix 2.4.

The simulation results show that an oscillatory input with the same frequency as a mode of the system matrix minimized the estimation error when the system had a single mode. Figure 4a shows the system matrix and its eigenvalues, which correspond to a single 10 Hz mode. We generated time series data using this system matrix and estimated the matrix. Figure 4b illustrates the changes in the sum of eigenvalues of the state covariance matrix $Σ_{X}$ , denoted as $\sum_{i} μ_{i}$ , the sum of reciprocal eigenvalues $\sum_{i} (1 / μ_{i})$ , and the estimation error defined by the Frobenius norm. The sum of eigenvalues is maximized by a 10 Hz sinusoidal input, while the sum of reciprocal eigenvalues is minimized. Consistent with these results, the estimation error is minimized by the 10 Hz sinusoidal input. As shown in Equation 9, the estimation error is proportional to the sum of the inverse eigenvalues. This tendency can be explained more clearly by plotting the eigenvalues as ellipsoids. We visualized the eigenvalues using eigenvalue-scaled ellipsoids, as shown in Figure 4c and d. Both the major and minor axes of the ellipsoids in Figure 4c are extended by the sinusoidal inputs, especially by the 10 Hz input. In contrast, the major and minor axes of the ellipsoids in Figure 4d are substantially shortened by the 10 Hz sinusoidal input.

Figure 4

Download asset Open asset

Error in the system matrix and eigenvalue-scaled ellipsoids.

(a) State transition matrix $A$ and its eigenvalues $λ_{i}$ . (b) System frequency responses, showing the sum of eigenvalues of the state covariance matrix $Σ_{X}$ (left), the sum of inverse eigenvalues of the state covariance matrix $Σ_{X}$ (middle), and the system matrix error (right). (c) Eigenvalue ellipsoids of $Σ_{X}$ under passive, 5 , 10 , and 15 Hz conditions. (d) Reciprocal-eigenvalue ellipsoids of $Σ_{X}$ under passive, 5 , 10 , and 15 Hz conditions.

When the system exhibits multiple modes, the input frequencies should be designed to reflect all these modes. We demonstrate this using the system matrix shown in Figure 5a, along with simulated time series and the corresponding matrix estimation. Because the system matrix contains two distinct pairs of eigenvalues (10 and 20 Hz), we constructed the evaluation input u as a combination of two frequencies. The sum of eigenvalue reciprocals is minimized only when the input contained both 10 and 20 Hz components; inputs with a single frequency yield larger values (Figure 5b). To further evaluate the effect of input design, we compared the two-frequency input with a flat-spectrum input (Figure 5c). When the total input energy was normalized, the two-frequency input achieves better identification performance than the flat-spectrum input, which is commonly employed in control and system identification studies. This result highlights the importance of tailoring the input spectrum to the system’s intrinsic dynamics rather than relying on uniform excitation. The superior performance of the two-frequency combination can be understood through the eigenvalue-scaled ellipsoids in Figure 5d–g. Figure 5d illustrates ellipsoids determined by the three dominant eigenvalues. The ellipsoid volume is expanded by three types of sinusoidal inputs (10 Hz, 20 Hz, and the combined input). However, the third axis (blue line) is not extended under single-frequency inputs, as shown in Figure 5e, leading to larger estimation errors. Small eigenvalues disproportionately increase the sum of reciprocal eigenvalues, as seen in Figure 5f and g. In these cases, the reciprocal eigenvalue-scaled ellipsoids are enlarged along the third axis (blue line) when only a single-frequency input is applied. By contrast, the combined input uniformly increases all eigenvalues, including the third, thereby successfully minimizing the reciprocal eigenvalue-scaled ellipsoid volume. These findings indicate that sinusoidal inputs should be designed to match multiple system modes rather than single modes. It should be noted that the true system modes cannot generally be known a priori. They must be identified through iterative experiments and refinement of the estimated system matrix. This practical procedure for optimal perturbation design is described in a later section (see Practical designing procedure for optimal perturbation inputs).

Figure 5

Download asset Open asset

Simulation of a multi-mode system and eigenvalue-scaled ellipsoids.

(a) System matrix $A$ and its eigenvalues; the system exhibits two damped oscillatory modes at 10 and 20 Hz. (b) Theoretical estimation error, defined as $\sum_{i} (1 / μ_{i})$ 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, $1 / μ_{i}$ ) for the same four conditions. (g) Relative-scale view of (f).

Demonstration of location tuning

This section illustrates how our theoretical framework enables location-specific tuning of perturbation inputs within neural networks. As established earlier, perturbations should be designed so that the variance becomes large in all directions, rather than being confined to specific modes. Two candidate strategies arise naturally: directly stimulating the nodes associated with heavily damped modes to target those specific modes, or stimulating a hub node whose outgoing connections propagate the perturbation across the entire network.

To examine this relationship in practice, we constructed a network as shown in Figure 6a and b. The network consists of seven nodes, structured into three oscillatory modes (Nodes 1–6) and one hub node (Node 7). It is designed to have three oscillatory modes at 15, 25, and 35 Hz. Figure 6c shows the damping rates $| λ_{A} |$ of the eigenvalues of A: the mode formed by Nodes 1 and 2 (15 Hz) is heavily damped, while those formed by Nodes 3–6 are moderately damped. Node 7 serves as a hub with only outgoing edges. Each node was individually perturbed by an impulse input. Here, an impulse input is defined as a Kronecker delta at $t = 0$ with fixed amplitude $α = 10$ , with no external input applied at any subsequent time step.

Figure 6

Download asset Open asset

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 $A$ . (b) Graph representation. (c) Damping rates $| λ_{A} |$ of the eigenvalues of $A$ . (d–f) 3D ellipsoids whose axes are proportional to $1 / μ_{5}$ , $1 / μ_{6}$ , $1 / μ_{7}$ (reciprocals of the three smallest eigenvalues of $Σ_{X}$ ; larger axes indicate greater estimation error): (d) passive, (e) impulse to Node 1, (f) impulse to Node 7. (g) Heatmap of $1 / μ_{i}$ 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.

Stimulating Node 7 (hub node) minimizes the total estimation error of A (Figure 6h), because its outgoing connections propagate the perturbation to all nodes and re-excite every dynamical mode. Figure 6d–f illustrate how the choice of input node shapes the smallest eigenvalues of $Σ_{X}$ . Because the 15 Hz mode (Nodes 1–2) is heavily damped, it decays rapidly and contributes almost no variance to $Σ_{X}$ under passive conditions; as a result, the fifth through seventh eigenvalues $μ_{5}$ – $μ_{7}$ are small and their reciprocals $1 / μ_{5}$ – $1 / μ_{7}$ are large, producing the elongated ellipsoid seen in Figure 6d. When Node 1 receives an impulse (Figure 6e), the eigenvalue associated with the 15 Hz mode grows and the $1 / μ_{5}$ and $1 / μ_{6}$ axes shrink; however, the $1 / μ_{7}$ axis—which corresponds to directions that Node 1 alone does not excite—remains large. Stimulating Node 7 (Figure 6f) re-excites all six oscillatory nodes through its outgoing connections, simultaneously reducing all three axes. Figure 6g confirms this pattern across all candidate nodes: the contribution $1 / μ_{i}$ remains uniformly small only when Node 7 receives the impulse. Figure 6h and i show that Node 7 achieves the smallest full-matrix estimation error, followed by Nodes 1 and 2, which directly drive the heavily damped 15 Hz mode.

As this simulation represents only one example of location design, its limitations and the corresponding countermeasures should be stated. In actual experiments, the most effective perturbation location depends on factors, such as hub-node connectivity and modal damping rates, and B itself may not always be known a priori. In such cases, approaches, such as iterative optimization of $t r (Σ_{X}^{- 1})$ (described in a later section) and joint estimation of A and B (Appendix 2.3) provide systematic alternatives. Nevertheless, the results presented here provide an intuitive guideline: stimulation directed at hub nodes or at nodes driving heavily damped modes effectively excites the full set of dynamical modes and minimizes estimation error.

Neural state classification

We demonstrate the effectiveness of our framework by applying it to a neural state classification problem. Neural state classification is crucial for diagnosing illnesses and assessing brain conditions in many experimental neuroscience settings (Casali et al., 2013; Hallett et al., 2017; Shikauchi et al., 2025; Angulo-Sherman et al., 2017). Some of these studies have already applied arbitrary perturbation inputs for classification, but not optimal perturbation inputs. By simulating such real-world applications, we demonstrate how well optimal perturbation design can contribute to neuroscience research.

We designed a neural network with clearly distinct task conditions and considered a simulation setting in which these conditions are classified using signals of a fixed duration. These distinct task conditions consist of five types, each defined by a unique linear dynamical system characterized by differing eigenvalue spectra and connectivity topologies of matrix A (Figure 7a). These task conditions are intended to mimic different cognitive or behavioral contexts. For example, in a typical motor task experiment, such conditions could correspond to motor execution or imagery involving the left or right hand, or resting state (Shikauchi et al., 2025; Angulo-Sherman et al., 2017). The neural signals were simulated under five different task conditions and two stimulation conditions: passive observation and external perturbation. Perturbation was applied as impulse-type inputs, such as TMS. The stimulus location was determined for each task condition by applying an impulse to each node and selecting the one that minimized $t r (Σ_{X}^{- 1})$ . The resulting time-series data are shown in Figure 7b. Using this data, we estimated the underlying dynamical system via a control-based identification approach presented in Equation 4, which corresponds to an estimation of functional connectivity.

Figure 7

Download asset Open asset

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 $A$ . 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 $A$ 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 $T$ increases from 20 to 200.

The simulation demonstrates that classification performance under perturbation is superior to that under the passive condition, both in visual inspection and in quantitative evaluation. As revealed by linear discriminant analysis (LDA) projection, task clusters in the passive condition exhibit substantial overlap, whereas perturbation leads to clear separation among the tasks (Figure 7c). Classification with LDA and 10-fold cross-validation shows the average accuracy is substantially higher in the perturbation condition (78.20%) compared to the passive condition (24.20%), as shown in Figure 7d. ROC curve analysis further confirms that perturbation substantially improves discriminability across all tasks and outperforms the random perturbation condition (Figure 7e). These results can be explained by Equation 8: the designed perturbation inputs increase the eigenvalues of the covariance matrix, including even the smallest ones, which in turn leads to a decrease in the estimation error of A.

To obtain estimates under the passive condition that are comparable to those derived under perturbation, it is necessary to experimentally observe extensive time-series data. Figure 7f illustrates the LDA projection and classification accuracy for different time-series lengths. The leftmost LDA plot $(T = 20)$ corresponds to the passive condition shown in Figure 7c, indicating that the estimation performance in the passive condition becomes comparable to that in the perturbation condition only when the time window reaches approximately $T = 200$ , a 10-fold increase compared to $T = 20$ . While the absolute duration depends on the interpretation of the time unit, such time windows may not be prohibitive in some experimental settings. Nevertheless, our results consistently show that passive observation requires substantially longer recordings to achieve comparable performance, highlighting the efficiency of the perturbation-based approach when the available data length is limited.

Neural state transitions via optimal control

This section demonstrates the practical utility of our framework by applying it to a neural state control problem (Gu et al., 2015; Karrer et al., 2020). Here, we estimate the system parameters from both passive and perturbation states. Then, we determine optimal control inputs which control the neural states to desired targets, and associated control costs. By using perturbation inputs, A is accurately estimated, which in turn allows precise determination of the optimal inputs and the corresponding controlled state transitions (Figure 8a). Moreover, use of the perturbation approach to derive the model allows the controllability Gramian (Gu et al., 2015; Shikauchi et al., 2025), as well as the estimation of control costs for each network area, to be more reliably assessed (Figure 8b and c). This simulation underscores the importance of accurate system identification for achieving optimal control and reliable estimation of neural functions. Throughout this section, we use the term perturbation input to refer to the input used for system identification and control input to refer to the input used for state transitions.

Figure 8

Download asset Open asset

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.

We designed a network as shown in Figure 9a to clearly show the differences in system identification between passive and perturbation states. The network consists of two disconnected groups. Nodes 1 and 2, as well as nodes 3, 4, and 5, are each connected separately, with no connections between these groups, as shown in matrix A. We assumed a known input matrix B in which nodes 3 and 4 cannot be directly controlled. In this scenario, the optimal control strategy to manipulate nodes 3 and 4 is to apply inputs to node 5. However, if the system identification is inaccurate, there is the possibility that control will be attempted through nodes 1 and 2, which are disconnected from nodes 3 and 4. These matrices were estimated from the time series signals generated by simulation. We generated passive and perturbation states to estimate the parameter matrices. The perturbation condition uses an impulse input with sufficient intensity ( $α = 10^{3}$ ). By using the estimated system model, we determined optimal control inputs, then evaluated the state transitions and associated costs. The controlled transition test was run with $T = 50$ . The control objective was to set nodes 3 and 4–25 while keeping all other nodes at 0 without any movement.

Figure 9

Download asset Open asset

Accurate system identification via perturbation enables precise optimal control.

(a) True matrix $A$ 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 ( $α = 10^{3}$ ). (d) Errors in the estimated matrices $A$ and $W$ and the control cost (average controllability).

The simulation results show that perturbation-based estimation enables the precise control of neural states. If the estimation of A is inaccurate, the state transitions under optimal control cannot be executed correctly (Figure 9b). Nodes 3 and 4 fail to reach the target, and unnecessary activations occur in nodes 1 and 2. Furthermore, unnecessary control inputs are required in nodes 1 and 2. On the other hand, when the model is correctly estimated using a high-strength impulse, the resulting control inputs enable accurate state transitions (Figure 9c). In this case, the strength of the impulse is correctly concentrated on node 5. As the error in matrix A increases, the estimation errors of the controllability Gramian W lead to inaccuracies in estimating the control cost (Figure 9d). Here, we adopt average controllability as an example measure of control cost (Gu et al., 2015; Moradi Amani et al., 2024). The key point here is that, since the controllability Gramian involves multiple products of A (Equation 21), even small errors in A can result in significant discrepancies. Therefore, when estimating optimal control input for neural transitions or control costs, it is more appropriate to identify those model parameters with at least arbitrary perturbations, and ideally with designed perturbations.

Practical designing procedure for optimal perturbation inputs

This section presents a practical framework for designing optimal perturbation inputs to estimate neural dynamics. The practical framework progressively refines the input signal through repeated cycles of model identification and perturbation input design. By leveraging this alternating scheme, each iteration utilizes the current model estimate to design a more informative input for the next round of identification. We also focus on perturbations characterized by a broad and uniform frequency distribution, a form commonly adopted in control engineering and conceptually analogous to transcranial random noise stimulation (tRNS) in neuroscience (Paulus, 2011; Paulus et al., 2016). As shown in Figure 5 and Appendix 2—figure 4, the designed composite-wave input allows for more accurate model estimation than the uniform-frequency input when their total input energies are equal.

Iterative refinement of both the perturbation design and the estimation process progressively improves the accuracy of A. The time-series data is collected from 32 points, where the matrix A (Figure 10a) is designed to have 16 oscillatory modes (i.e. 16 complex-conjugate eigenvalue pairs, yielding 32 eigenvalues in total). In this simulation, one stimulation session is applied to a single node at each design step. At each iteration, the stimulation is designed as a composite-frequency sinusoidal input encompassing all modes of the estimated $\hat{A}$ , and the target node of the perturbation is determined through numerical optimization that minimizes $t r (Σ_{X}^{- 1})$ . Appendix 2—table 1 summarizes the simulation conditions, including measurement duration and number of trials used in each simulation. The stimulation duration was set equal to the measurement duration.

Figure 10

Download asset Open asset

Iterative perturbation design converges to the theoretical optimum in a 32-channel neural network simulation.

(a) True matrix $A$ 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 ( $A_{passive}$ ) is used to design the first perturbation input ( $u_{design1}$ ), which yields a better estimate ( $A_{design1}$ ), and the cycle repeats.

The following procedure outlines a practical approach (Figure 10c) for this simulation, aiming to precisely estimate the matrix by designing an optimal perturbation input (assuming frequency stimulation, such as tACS).

Record spontaneous activity as the passive state and estimate A_passive.
Based on the estimated matrix, the optimal frequency and target node of the perturbation input are determined through numerical optimization following Equation 9.
The designed perturbation u_design1 is applied to the system, and a more precise matrix A_design1 is identified.
The obtained matrix A_design1 is again used to design a new perturbation input u_design2 and estimate a more accurate matrix A_design2. This iterative process can be continued until convergence criteria, such as changes in A, are met.

The proposed framework improved estimation accuracy iteration by iteration. As shown in Figure 10b, the estimation error converges toward the theoretical optimum (dashed line) with each design iteration. Here, the theoretical optimum is obtained from a simulation in which the input is constructed as a composite sinusoid covering all dynamical modes of the true matrix A, and the stimulation location is numerically selected to minimize $t r (Σ_{X}^{- 1})$ . The error from the second iteration (second design) is smaller than that from the first (first design), which in turn is substantially lower than the initial estimate. This convergence holds both when starting from the passive condition and when starting from a flat-spectrum input (a practical approach analogous to tRNS). While the flat-spectrum input (blue line, Iteration 1) provided a better starting point than passive data, the designed input from the initial flat-spectrum estimate further reduced the error by the second iteration. Furthermore, the converged estimation error was lower than that of a random stimulation baseline—in which a single randomly chosen node was stimulated using a flat-spectrum (all-frequency) input throughout all iterations (green line)—whose estimation error remained approximately constant across iterations. This contrast directly demonstrates that the benefit of the proposed framework arises specifically from model-based optimization of both the spatial (node selection) and spectral (frequency composition) stimulation parameters. These results should be interpreted as a case-specific illustration rather than a universal gain, as the magnitude of improvement depends on factors, such as network structure, recording duration, noise level, and stimulation design. This simulation demonstrates that our perturbation design framework enables the step-by-step refinement of system identification even without prior knowledge of the system.

Discussion

This study addressed the challenge of designing optimal perturbations for effectively identifying neural system dynamics. We introduced a framework for estimating the optimal perturbation input for identification of neural systems. Our findings demonstrated that incorporating perturbation inputs, including TMS, tDCS, and tACS, significantly improves identification accuracy. Specifically, alignment of their parameters with the intrinsic frequencies of matrix A, high input intensity, and targeted inputs directed at nodes that efficiently excite heavily damped dynamical modes (either nodes directly associated with those modes or hub nodes whose outgoing connections propagate the perturbation throughout the network) enhances system identification. Furthermore, we outlined an approach for designing the optimal perturbation input by iterating system identifications with perturbations, and demonstrated that the approach progressively decreases parameter estimation error. Beyond the parameter identification of neural systems, our findings provide insights into how these estimated parameters influence optimal control theory in neuroscience. These results emphasize the potential of perturbation input design to advance understanding of neural dynamics.

The assumption of linearity and first-order autoregressiveness in neural dynamics warrants discussion. Biological signals often exhibit complex and nonlinear interactions that cannot be fully captured by linear models. Nevertheless, prior studies have demonstrated that nonlinear behaviors can frequently be approximated using linear models (Honey et al., 2009; Gu et al., 2015; Fernández Galán, 2008). Linear models are particularly effective for providing accurate approximations of nonlinear systems within a specific operating range (Khalil, 2002). Therefore, establishing theoretical foundations based on linear assumptions can contribute to the development of theories for nonlinear systems or be effectively utilized in their advancement. To extend these theories into the nonlinear domain, it would be necessary to incorporate advanced frameworks, such as the Koopman operator (Koopman, 1931; Chow et al., 2024) and system identification methods leveraging machine learning techniques (Suk et al., 2016; Glaser et al., 2019). On the other hand, discussions regarding higher-order VAR models are relatively straightforward due to the simplicity of the estimation method. Previous studies involving actual functional data have employed VAR models of an order greater than one (Tseng et al., 1995; Chang et al., 2012; Shakeel et al., 2021). The framework proposed in this paper can be naturally and meaningfully extended to higher-order VAR models, broadening its applicability to more complex temporal dependencies.

Safety and experimental constraints play a pivotal role in the practical implementation of this framework, influencing both the scope and design of potential applications. In biological systems, stimulation parameters, such as amplitude, frequency, and location, must adhere to stringent safety thresholds to avoid adverse effects, such as tissue damage or unintended physiological responses (Rocchi et al., 2022; Bikson et al., 2016). For stimulus intensity, our proposed framework suggested that a stronger stimulus improves system identification; however, an excessively strong stimulus poses safety risks (Antal et al., 2017). With respect to sinusoidal input, stimuli are typically applied within ranges corresponding to neural activity. The optimal input frequencies derived from the proposed theory are expected to align with these ranges, suggesting that the proposed framework can be implemented without major complications (Antal et al., 2017; Mager et al., 2018). Furthermore, the selection of stimulation sites is often dictated by experimental accessibility or ethical considerations (Rossi et al., 2009; Heinrichs, 2012). Based on these constraints, it is necessary to determine the most appropriate stimulation parameters to achieve optimal system identification.

Two directions warrant further investigation: extending the framework to partial observability, and validating it through stimulation experiments. In experimental neuroscience, recordings are often limited to a subset of neural populations, resulting in partial observability. A growing body of work has leveraged delay-embedding techniques, represented by Takens' embedding theorem (Takens, 1980), to reconstruct hidden dynamics from partial observations (Arbabi and Mezić, 2017; Brunton et al., 2017; Ostrow et al., 2024; Huang et al., 2025). Applying such techniques enables the estimation of the full connectivity matrix, thereby extending our framework to settings with partial observability. The second direction concerns experimental validation. Validating a theoretical framework through experimental design is essential for bridging the gap between theory and practice. Recent research involving some of the present authors has demonstrated that TMS can facilitate the discrimination of neural states (Shikauchi et al., 2025). Future experiments can build on this finding by incorporating passive, arbitrary, and designed optimal stimuli to estimate the connectivity matrix, and then comparing the ability of these matrices to distinguish different neural states. If these proposed evaluations demonstrate the practical effectiveness of our theory, it could then be applied to actual experiments. As shown in Figure 10, a preliminary connectivity matrix is obtained through preliminary experiments to design optimal perturbations. These perturbations would then be applied in the main experiments, improving the estimation of connectivity matrices and neural dynamics. Experiments using our approach will enable the more precise and comprehensive analysis of brain and neural function, and in turn facilitate valuable new insights into human cognition and behavior.

Methods

State vector simulation and model estimation

To investigate appropriate perturbation inputs for the identification of neural systems, we constructed neural dynamics using matrices A and B to generate the temporal evolution of state variables. From the generated dynamics, we estimated the model matrix A, while assuming that the model matrix B was known.

Given an initial condition $x (0)$ , the neural dynamics were generated according to the true matrices A and B and the governing equation for $x (t)$ (Equation 1), with a predefined noise time series $ξ (t)$ added at each time step. For the discrete-time simulations, the sampling rate was set to an appropriate value for each simulation. The model matrix A was then estimated from the generated time-series states $x (t)$ and applied perturbation inputs $u (t)$ using the OLS.

Both the generation of neural dynamics and the model estimation were repeated for a predefined number of trials, with variations in initial conditions and noise sequences. The estimation errors of the matrices were quantified using the Frobenius norm. Details of the simulation parameters are provided in Appendix 2.1.

Eigenvector alignment for eigenvalue-scaled ellipsoids

To examine the relationship between the eigenvalues $μ_{i}$ of the state covariance matrix $Σ_{X}$ and the input frequencies, we plotted an eigenvalue-scaled ellipsoid, thereby providing an intuitive representation of optimal input design. However, the eigenvectors of $Σ_{X}$ are not fixed; they may vary depending on the applied input, which complicates direct comparison across conditions. To enable consistent comparison of eigenvalues across input conditions, we reordered eigenvalues according to the similarity of their associated eigenvectors to those obtained in the passive condition, which served as the reference basis. For each reference eigenvector $t_{i}$ , we identified the most similar eigenvector $u_{j}$ from the input condition by maximizing the absolute inner product

j^{*} = \arg max_{j} | t_{i}^{⊤} u_{j} | .

The eigenvalue $μ_{j^{*}}$ corresponding to $u_{j^{*}}$ was then assigned to the (i)-th position of the reordered list. Each $u_{j}$ was used only once, ensuring a one-to-one correspondence between reference and input eigenvectors. This alignment resolves the permutation and sign ambiguities inherent in eigendecomposition, and allows eigenvalues to be compared across conditions along a common axis.

Optimal control and network controllability

Accurately estimating dynamical models under perturbation not only facilitates precise inference of the neural dynamics but also contributes to the accurate estimation of control inputs for neural state transitions. Recently, the control theory which utilizes the controllability Gramian and control costs is applied for neuroscience, and provides insights into the efficiency and feasibility of inducing specific neural state transitions (Gu et al., 2015; Karrer et al., 2020) with some limitations (Tu et al., 2018; Suweis et al., 2019). For clarity, we refer to the input for system identification as the perturbation input and the input for state control in the optimal control theory as the control input.

The optimal control input is derived under the condition of assuming the linear system is controllable. We consider the sequence of control inputs that minimizes the following cost function (input energy minimization) while driving the state from the initial state $x_{0}$ to the terminal state $x_{T}$ :

M ({u (t)}_{t = 0}^{T - 1}) = \sum_{t = 0}^{T - 1} u (t)^{⊤} u (t),

under

x (0) = x_{0}, x (T) = x_{T} .

Here, the initial and terminal state conditions represent constraints on the optimization problem. This problem is expressed as follows:

u (t)^{*} = \underset{u}{arg min} M ({u (t)}_{t = 0}^{T - 1}) .

The constrained optimization problem can be solved using Lagrange multipliers. The optimal control input is given by

u (t)^{*} = B^{⊤} (A^{⊤})^{T - t - 1} W^{- 1} (x (T) - A^{T} x (0)),

where W is the controllability Gramian, defined as:

W = \sum_{t = 0}^{T - 1} A^{t} B B^{⊤} (A^{⊤})^{t} .

The controllability Gramian plays a pivotal role in determining the optimal control and control input cost. It has been used to identify the functional roles of individual brain regions (Gu et al., 2015; Shikauchi et al., 2025).

Average controllability for control cost

While the optimal control problem (Equation 19) calculates the specific input energy for a given state transition, a more general, state-independent metric is often used to characterize the system’s overall controllability. This metric, average controllability, is a state-independent metric used to quantify the system’s overall ease of control (Gu et al., 2015; Moradi Amani et al., 2024). It is defined as the trace of the controllability Gramian:

Average Controllability = tr (W) = \sum_{i = 1}^{N} λ_{i} (W)

A larger trace indicates that the system is, on average, more controllable (i.e. can be moved to various states with less input energy). This metric is, therefore, used as a measure of control efficiency.

Appendix 1

Theoretical details

1.1 Asymptotic property of the estimator for A

This section derives the asymptotic estimation error of the OLS estimator $\hat{A}$ for the system matrix A. The primary goal is to obtain a first-order asymptotic approximation for the Mean Squared Error (MSE) of the estimation, quantified by the expected squared Frobenius norm, $E [‖ \hat{A} - A ‖_{F}^{2}]$ .

We consider the linear dynamical system

x (t + 1) = A x (t) + B u (t) + ξ (t),

where the noise $ξ (t)$ is i.i.d. over time with zero mean and covariance matrix $Σ_{ξ} = E [ξ (t) ξ (t)^{⊤}] \in R^{n \times n} ≻ 0$ .

We estimate A by OLS from the data matrices $X \in R^{n \times T}$ , $\tilde{Y} := Y - B U \in R^{n \times T}$ , giving $\hat{A} = \tilde{Y} X^{†}$ .

The asymptotic distribution of the vectorized estimation error of the OLS estimator $\hat{A}$ is given by Hamilton, 1994, Proposition 11.1 and Appendix 11.A:

\sqrt{T} (v e c (\hat{A}) - v e c (A)) \overset{L}{\to} N (0, Σ_{A s y m}),

where T is the number of samples, $\overset{L}{\to}$ denotes convergence in distribution, and $Σ_{A s y m}$ is the asymptotic covariance matrix. Following the notation in the original text (and its assumptions), this covariance is:

Σ_{A s y m} = Σ_{ξ} \otimes Σ_{X}^{- 1},

where $\otimes$ denotes the Kronecker product [Hamilton, 1994, Appendix A.4, p.732], $Σ_{ξ} \in R^{n \times n}$ is the noise covariance, and $Σ_{X} = E [x (t) x (t)^{⊤}] \in R^{n \times n}$ is the state covariance matrix, so that $Σ_{A s y m} \in R^{n^{2} \times n^{2}}$ .

We are interested in the estimation error measured by the squared Frobenius norm, $‖ \hat{A} - A ‖_{F}^{2}$ . Let us define the normalized error vector $z_{T}$ based on (Equation A1.2):

z_{T} := \sqrt{T} (v e c (\hat{A}) - v e c (A)) .

We know that $z_{T} \overset{L}{\to} z$ , where $z \sim N (0, Σ_{A s y m})$ is an $n^{2}$ -dimensional Gaussian vector.

The squared Frobenius norm can be expressed as the squared Euclidean norm of the vectorized error:

‖ \hat{A} - A ‖_{F}^{2} = ‖ v e c (\hat{A} - A) ‖_{2}^{2} .

Multiplying by T, we can relate this quantity to $z_{T}$ :

T \cdot ‖ \hat{A} - A ‖_{F}^{2} = {(\sqrt{T} \cdot ‖ v e c (\hat{A} - A) ‖_{2})}^{2} = z_{T}^{⊤} z_{T} .

By the continuous mapping theorem, since $g (z) = z^{⊤} z$ is a continuous function, the convergence in distribution of $z_{T}$ implies the convergence in distribution of $z_{T}^{⊤} z_{T}$ :

T \cdot ‖ \hat{A} - A ‖_{F}^{2} \overset{L}{\to} z^{⊤} z, where z \sim N (0, Σ_{A s y m}) .

Its limit of the expected value equals the expectation of the limiting distribution $z^{⊤} z$ :

lim_{T \to \infty} E [T \cdot ‖ \hat{A} - A ‖_{F}^{2}] = E [z^{⊤} z]

= E [t r (z z^{⊤})]

= t r (E [z z^{⊤}]) .

By definition, $E [z z^{⊤}]$ is the covariance matrix of the limiting distribution z, which is $Σ_{A s y m}$ .

lim_{T \to \infty} E [T \cdot ‖ \hat{A} - A ‖_{F}^{2}] = t r (Σ_{A s y m})

= t r (Σ_{ξ} \otimes Σ_{X}^{- 1}) .

Using the trace property $t r (A \otimes B) = t r (A) t r (B)$ , we obtain:

lim_{T \to \infty} E [T \cdot ‖ \hat{A} - A ‖_{F}^{2}] = t r (Σ_{ξ}) t r (Σ_{X}^{- 1}) .

This result implies that for a sufficiently large T, the expected estimation error can be approximated as:

E [‖ \hat{A} - A ‖_{F}^{2}] \approx \frac{1}{T} t r (Σ_{ξ}) t r (Σ_{X}^{- 1}) .

This expression provides a simple interpretation of the estimation error. It confirms that the error is inversely proportional to the sample size T, meaning the error decreases as more data is collected. The error is also directly proportional to the total noise power, $t r (Σ_{ξ})$ , indicating that more noise leads to more error. Furthermore, the error is proportional to $t r (Σ_{X}^{- 1})$ , which acts as an inverse measure of signal strength; a weaker signal $x (t)$ (i.e. less variation) leads to a larger error. This relationship resembles the inverse of the Signal-to-Noise ratio, where the error is governed by the balance between noise, signal strength, and the amount of data.

1.2 State vectors considering typical inputs used in neuroscience

1.2.1 Sinusoidal Input

To determine which frequencies most effectively enhance system identification accuracy, we now analyze the time evolution of the state vector $x (t)$ in response to a sinusoidal input. This analysis is the necessary first step towards the goal established in the main text (i.e. sufficiently exciting the system’s dynamics to increase the eigenvalues of the covariance matrix). By understanding how $x (t)$ itself behaves under different frequencies, we can identify which inputs are most effective. We, therefore, examine frequency external inputs, represented in neuroscience by tACS, to derive the optimal input frequency based on the system’s dynamic response.

In the case where the perturbation consists of a composite sinusoidal input comprising multiple frequencies, the input can be written as

u (k) = \sum_{l = 1}^{L} \cos (ω_{l} k) u_{0},

where $ω_{l}$ are the angular frequencies.

The general solution which we denote as Equation 10 is:

x (t) = A^{t} x (0) + \sum_{k = 0}^{t - 1} A^{t - 1 - k} B u (k) + \sum_{k = 0}^{t - 1} A^{t - 1 - k} ξ (k)

Inserting Equation A1.15 into Equation 10 yields:

x (t) = A^{t} x (0) + \sum_{k = 0}^{t - 1} A^{t - 1 - k} B \sum_{l = 1}^{L} \cos (ω_{l} k) u_{0} + \sum_{k = 0}^{t - 1} A^{t - 1 - k} ξ (k)

Under the assumption that the system is stable (all eigenvalues $λ_{d}$ of A satisfy $| λ_{d} | < 1$ ), the second term of Equation A1.17 is the critical component to analyze.

By performing diagonalization of A as $A = V Λ V^{- 1}$ , where V is the matrix of eigenvectors, $Λ$ is the diagonal matrix of eigenvalues $λ_{d}$ ,... The spectral decomposition of $A^{m}$ is $A^{m} = \sum_{d = 1}^{n} λ_{d}^{m} v_{d} w_{d}^{⊤}$ , where $v_{d}$ denotes the d-th column vector of V, and $w_{d}^{⊤}$ denotes the d-th row vector of $V^{- 1}$ . The second term of Equation A1.17 can be expressed as follows:

\sum_{k = 0}^{t - 1} A^{t - 1 - k} B \sum_{l = 1}^{L} \cos (ω_{l} k) u_{0}

= \sum_{k = 0}^{t - 1} V Λ^{t - 1 - k} V^{- 1} B \sum_{l = 1}^{L} \cos (ω_{l} k) u_{0}

= \sum_{k = 0}^{t - 1} \sum_{d = 1}^{n} λ_{d}^{t - 1 - k} v_{d} w_{d}^{⊤} B \sum_{l = 1}^{L} \cos (ω_{l} k) u_{0}

Using Euler’s formula ( $\cos x = (e^{i x} + e^{- i x}) / 2$ ) and swapping the order of summation (which is permissible by linearity), we get:

\sum_{k = 0}^{t - 1} \sum_{d = 1}^{n} λ_{d}^{t - 1 - k} v_{d} w_{d}^{⊤} B \sum_{l = 1}^{L} \cos (ω_{l} k) u_{0}

= \sum_{k = 0}^{t - 1} \sum_{d = 1}^{n} \sum_{l = 1}^{L} \frac{e^{i ω_{l} k} + e^{- i ω_{l} k}}{2} λ_{d}^{t - 1 - k} v_{d} w_{d}^{⊤} B u_{0}

= \frac{1}{2} \sum_{l = 1}^{L} \sum_{d = 1}^{n} (\sum_{k = 0}^{t - 1} λ_{d}^{t - 1 - k} e^{i ω_{l} k} + \sum_{k = 0}^{t - 1} λ_{d}^{t - 1 - k} e^{- i ω_{l} k}) v_{d} w_{d}^{⊤} B u_{0}

We now introduce the polar form of the eigenvalues, $λ_{d} = r_{d} e^{i θ_{d}}$ . Given the stability assumption, we have $r_{d} = | λ_{d} | < 1$ . This yields:

= \frac{1}{2} \sum_{l = 1}^{L} \sum_{d = 1}^{n} (\frac{λ_{d}^{t} - e^{i ω_{l} t}}{λ_{d} - e^{i ω_{l}}} + \frac{λ_{d}^{t} - e^{- i ω_{l} t}}{λ_{d} - e^{- i ω_{l}}}) v_{d} w_{d}^{⊤} B u_{0}

= \frac{1}{2} \sum_{l = 1}^{L} \sum_{d = 1}^{n} (\frac{r_{d}^{t} e^{i θ_{d} t} - e^{i ω_{l} t}}{r_{d} e^{i θ_{d}} - e^{i ω_{l}}} + \frac{r_{d}^{t} e^{i θ_{d} t} - e^{- i ω_{l} t}}{r_{d} e^{i θ_{d}} - e^{- i ω_{l}}}) v_{d} w_{d}^{⊤} B u_{0}

= \frac{1}{2} \sum_{l = 1}^{L} \sum_{d = 1}^{n} (\frac{(r_{d}^{t} e^{i θ_{d} t} - e^{i ω_{l} t}) (r_{d} e^{- i θ_{d}} - e^{- i ω_{l}})}{(r_{d} e^{i θ_{d}} - e^{i ω_{l}}) (r_{d} e^{- i θ_{d}} - e^{- i ω_{l}})} + \frac{(r_{d}^{t} e^{i θ_{d} t} - e^{- i ω_{l} t}) (r_{d} e^{- i θ_{d}} - e^{i ω_{l}})}{(r_{d} e^{i θ_{d}} - e^{- i ω_{l}}) (r_{d} e^{- i θ_{d}} - e^{i ω_{l}})}) v_{d} w_{d}^{⊤} B u_{0}

\begin{array}{ll} = \frac{1}{2} \sum_{l = 1}^{L} \sum_{d = 1}^{n} (\frac{r_{d}^{t + 1} e^{i θ_{d} (t - 1)} - r_{d}^{t} e^{i (θ_{d} t - ω_{l})} - r_{d} e^{i (ω_{l} t - θ_{d})} + e^{i ω_{l} (t - 1)}}{r_{d}^{2} - 2 r_{d} \cos (θ_{d} - ω_{l}) + 1} \\ + \frac{r_{d}^{t + 1} e^{i θ_{d} (t - 1)} - r_{d}^{t} e^{i (θ_{d} t + ω_{l})} - r_{d} e^{i (- ω_{l} t - θ_{d})} + e^{- i ω_{l} (t - 1)}}{r_{d}^{2} - 2 r_{d} \cos (θ_{d} + ω_{l}) + 1}) v_{d} w_{d}^{⊤} B u_{0}, \end{array}

For brevity, we define the following terms:

Q (ω_{l}, θ_{d}) = \frac{r_{d}^{t + 1} e^{i θ_{d} (t - 1)} - r_{d}^{t} e^{i (θ_{d} t - ω_{l})} - r_{d} e^{i (ω_{l} t - θ_{d})} + e^{i ω_{l} (t - 1)}}{r_{d}^{2} - 2 r_{d} \cos (θ_{d} - ω_{l}) + 1}

F_{d} = v_{d} w_{d}^{⊤} B u_{0}

The second term is:

\sum_{k = 0}^{t - 1} A^{t - 1 - k} B \sum_{l = 1}^{L} \cos (ω_{l} k) u_{0} = \frac{1}{2} \sum_{l = 1}^{L} \sum_{d = 1}^{n} (Q (ω_{l}, θ_{d}) + Q (- ω_{l}, θ_{d})) F_{d}

Since the matrix A is real, its complex eigenvalues and eigenvectors must appear in complex conjugate pairs. If $λ_{d} = r_{d} e^{i θ_{d}}$ (with $θ_{d} \neq 0$ ) is an eigenvalue, then $λ_{d^{'}} = \bar{λ_{d}} = r_{d} e^{- i θ_{d}}$ is also an eigenvalue for some index $d^{'}$ . Thus, $θ_{d^{'}} = - θ_{d}$ . Furthermore, since B and $u_{0}$ are also real, the vector $F_{d} = v_{d} w_{d}^{⊤} B u_{0}$ associated with $λ_{d}$ and the vector $F_{d^{'}} = v_{d^{'}} w_{d^{'}}^{⊤} B u_{0}$ associated with $λ_{d^{'}}$ are also complex conjugates: $F_{d^{'}} = \bar{F_{d}}$ . We can use this property to simplify the sum.

Q (- ω_{l}, θ_{d^{'}}) F_{d^{'}} = Q (- ω_{l}, - θ_{d}) \bar{F_{d}}

= \frac{r_{d}^{t + 1} e^{- i θ_{d} (t - 1)} - r_{d}^{t} e^{i (- θ_{d} t + ω_{l})} - r_{d} e^{i (- ω_{l} t + θ_{d})} + e^{- i ω_{l} (t - 1)}}{r_{d}^{2} - 2 r_{d} \cos (- θ_{d} + ω_{l}) + 1} \bar{F_{d}}

= \bar{\frac{r_{d}^{t + 1} e^{i θ_{d} (t - 1)} - r_{d}^{t} e^{i (θ_{d} t - ω_{l})} - r_{d} e^{i (ω_{l} t - θ_{d})} + e^{i ω_{l} (t - 1)}}{r_{d}^{2} - 2 r_{d} \cos (θ_{d} - ω_{l}) + 1} F_{d}}

= \bar{Q (ω_{l}, θ_{d}) F_{d}}

Similarly, $Q (- ω_{l}, θ_{d}) F_{d} = \bar{Q (ω_{l}, θ_{d^{'}}) F_{d^{'}}}$

For real eigenvalues ( $θ_{d} = 0$ ), $F_{d}$ is real, and we find:

Q (- ω_{l}, 0) F_{d} = \frac{r_{d}^{t + 1} - r_{d}^{t} e^{i ω_{l}} - r_{d} e^{- i ω_{l} t} + e^{- i ω_{l} (t - 1)}}{r_{d}^{2} - 2 r_{d} \cos (ω_{l}) + 1} F_{d}

= \bar{\frac{r_{d}^{t + 1} - r_{d}^{t} e^{- i ω_{l}} - r_{d} e^{i ω_{l} t} + e^{i ω_{l} (t - 1)}}{r_{d}^{2} - 2 r_{d} \cos (ω_{l}) + 1} F_{d}}

= \bar{Q (ω_{l}, 0) F_{d}}

We divide the set of indices of eigenvalues into three subsets as follows. Let $D_{R}$ denote the set of indices corresponding to real eigenvalues ( $λ_{d} \in R$ ). Among complex eigenvalues, let $D_{C}^{+}$ denote the set of indices whose imaginary part is positive, and let $D_{C}^{-}$ denote the set of indices whose imaginary part is negative.

\frac{1}{2} \sum_{d = 1}^{n} (Q (ω_{l}, θ_{d}) + Q (- ω_{l}, θ_{d})) F_{d} = \frac{1}{2} \sum_{d \in D_{R}} (Q (ω_{l}, 0) F_{d} + \bar{Q (ω_{l}, 0) F_{d}})

+ \frac{1}{2} \sum_{d \in D_{C}^{+}} (Q (ω_{l}, θ_{d}) F_{d} + Q (- ω_{l}, θ_{d}) F_{d})

+ \frac{1}{2} \sum_{d^{'} \in D_{C}^{-}} (Q (ω_{l}, θ_{d^{'}}) F_{d^{'}} + Q (- ω_{l}, θ_{d^{'}}) F_{d^{'}})

Now, we pair the terms from $D_{C}^{+}$ and $D_{C}^{-}$ . For each $d \in D_{C}^{+}$ , there is a corresponding $d^{'} \in D_{C}^{-}$ such that $θ_{d^{'}} = - θ_{d}$ and $F_{d^{'}} = \bar{F_{d}}$ . Using the properties $Q (- ω_{l}, θ_{d}) F_{d} = \bar{Q (ω_{l}, θ_{d^{'}}) F_{d^{'}}}$ and $Q (- ω_{l}, θ_{d^{'}}) F_{d^{'}} = \bar{Q (ω_{l}, θ_{d}) F_{d}}$ , the sum of the $D_{C}^{+}$ and $D_{C}^{-}$ terms becomes:

\frac{1}{2} \sum_{d \in D_{C}^{+}} (Q (ω_{l}, θ_{d}) F_{d} + \bar{Q (ω_{l}, θ d^{'}) F_{d^{'}}} + Q (ω_{l}, θ d^{'}) F_{d^{'}} + \bar{Q (ω_{l}, θ_{d}) F_{d}})

= \frac{1}{2} \sum_{d \in D_{C}^{+}} ((Q_{d} F_{d} + \bar{Q_{d} F_{d}}) + (Q d^{'} F_{d^{'}} + \bar{Q d^{'} F_{d^{'}}}))

= \sum_{d \in D_{C}^{+}} Re [Q (ω_{l}, θ_{d}) F_{d}] + \sum_{d^{'} \in D_{C}^{-}} Re [Q (ω_{l}, θ_{d^{'}}) F_{d^{'}}]

Combining this with the $D_{R}$ sum, which is $\sum_{d \in D_{R}} Re [Q (ω_{l}, 0) F_{d}]$ , we arrive at:

\sum_{d = 1}^{n} Re [Q (ω_{l}, θ_{d}) F_{d}]

This is the real part of the product of the terms $Q (ω_{l}, θ_{d})$ and $F_{d}$ defined in Equations A1.28; A1.29. Therefore, the second term of the $x (t)$ can be expressed as:

\sum_{k = 0}^{t - 1} A^{t - 1 - k} B \sum_{l = 1}^{L} \cos (ω_{l} k) u_{0}

= \sum_{l = 1}^{L} \sum_{d = 1}^{n} Re [\frac{r_{d}^{t + 1} e^{i θ_{d} (t - 1)} - r_{d}^{t} e^{i (θ_{d} t - ω_{l})} - r_{d} e^{i (ω_{l} t - θ_{d})} + e^{i ω_{l} (t - 1)}}{r_{d}^{2} - 2 r_{d} \cos (θ_{d} - ω_{l}) + 1} v_{d} w_{d}^{⊤} B u_{0}],

Therefore, the state vector $x (t)$ is decomposed into a passive and a difference term as follows,

x (t) = x_{passive} (t) + x_{diff} (t),

x_{passive} (t) = A^{t} x (0) + \sum_{k = 0}^{t - 1} A^{t - 1 - k} ξ (k)

x_{diff} (t) = \sum_{l = 1}^{L} \sum_{d = 1}^{n} Re [\frac{r_{d}^{t + 1} e^{i θ_{d} (t - 1)} - r_{d}^{t} e^{i (θ_{d} t - ω_{l})} - r_{d} e^{i (ω_{l} t - θ_{d})} + e^{i ω_{l} (t - 1)}}{r_{d}^{2} - 2 r_{d} \cos (θ_{d} - ω_{l}) + 1} v_{d} w_{d}^{⊤} B u_{0}]

where $x_{passive} (t)$ represents the state x that would be realized if no perturbation were applied, i.e., the state that would evolve according to Equation 10 with $u = 0$ in the main text. $x_{diff} (t)$ represents the change in the state trajectory induced by the specific perturbation input $u (k) = \sum_{l = 1}^{L} \cos (ω_{l} k) u_{0}$ . As defined in the previous section, $λ_{d} = r_{d} e^{i θ_{d}}$ is the polar form of the d-th eigenvalue, the vector $v_{d}$ denotes the right eigenvector of A associated with $λ_{d}$ , and $w_{d}^{⊤}$ denotes the corresponding left eigenvector. The numerator of the fraction in Equation A1.46 consists of four distinct oscillatory terms. The resonant behavior of this expression is primarily governed by the denominator, $D (ω_{l}, θ_{d}) = r_{d}^{2} - 2 r_{d} \cos (θ_{d} - ω_{l}) + 1$ . The gain term $1 / D (ω_{l}, θ_{d})$ has a single extremum (a maximum) with respect to $ω_{l}$ , which occurs when $\cos (θ_{d} - ω_{l}) = 1$ , i.e., when the input frequency $ω_{l}$ matches the eigenfrequency (argument) $θ_{d}$ of the system.

The restriction to a common spatial pattern $u_{0}$ in the input $u (k) = \sum_{l = 1}^{L} \cos (ω_{l} k) u_{0}$ is a simplifying assumption made to present the theory of this paper clearly, and it can be naturally replaced by frequency-specific patterns $u_{l}$ . The simplification keeps the formulation consistent with the impulse and step inputs, and it isolates the effect of the input frequencies from that of the stimulation location, thereby keeping the presentation in the main text intuitive. Substituting $u_{0} \to u_{l}$ gives a more general expression, but it leaves the denominator $D (ω_{l}, θ_{d})$ untouched, so the conclusion that the input frequencies should match the mode frequencies is unchanged. Optimization that also includes the stimulation location can be carried out by numerically minimizing Equation A1.14, so the practical estimation procedure and design strategy remain unchanged.

Appendix 1.2.2 Impulse and step input

Impulse and step inputs follow simple rules. Let α represent the scaling factor for impulse input. Then, input $u (k)$ is defined as:

u (k) = {\begin{cases} α u_{0} & if k = 0, \\ 0 & if k > 0. \end{cases}

Here, $u_{0} \in R^{m}$ represents the fundamental input direction applied to the system, determining which input channel receives the impulse. We obtain the decomposition:

x_{passive} (t) = A^{t} x (0) + \sum_{k = 0}^{t - 1} A^{t - 1 - k} ξ (k), x_{diff} (t) = α A^{t - 1} B u_{0} .

For a scaled step input, $u (k)$ is defined as

u (k) = β u_{0}, \forall k \geq 0,

where β is the input strength. We then get the decomposition:

x_{passive} (t) = A^{t} x (0) + \sum_{k = 0}^{t - 1} A^{t - 1 - k} ξ (k), x_{diff} (t) = β \sum_{k = 0}^{t - 1} A^{t - 1 - k} B u_{0} .

From these solutions, we can infer that the input-driven terms ( $x_{diff} (t)$ ) are linear in the input strengths α and β. Consequently, the contribution of these inputs to the state covariance matrix (i.e. the outer product of the state vector) will scale quadratically with $α^{2}$ and $β^{2}$ . This analysis, which suggests that such inputs can increase the rank or magnitude of the covariance matrix eigenvalues, is explored further in the main text.

Appendix 2

Experimental conditions and results

2.1 Simulation details

The specific conditions used for generating neural dynamics in each section are summarized in Appendix 2—table 1. These include factors, such as sampling frequency, trial length, and perturbation input strength, which are crucial for designing reproducible simulations and evaluations. These parameters are chosen to illustrate the theoretical insights as clearly as possible, rather than to characterise the framework’s behaviour across generic regimes.

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⁻³²	1
Figures 4 and 5 (Sinusoidal inputs)	100	50	5000	10,000	10⁻²	1
Figure 6 (Location tuning)	1000	1	1000	200	10⁻²	10
Figure 7 (Neural state classification)	10	2	20	100	10⁻²	10²
Figure 9 (Optimal control)	50	1	50	100	10⁻²	10³
Figure 10 (Practical designing procedure)	100	25	2500	500	10⁻²	0.1
Appendix 2—figure 1 (Nonlinear dynamics)	1000	1	1000	30	10⁻³²	0.15–0.3
Appendix 2—figure 2(B-Matrix uncertainty)	1000	1	1000	30	10⁻³²	0.15–0.3
Appendix 2—figure 3 (Impulse/step)	100	50	5000	1	10⁻²	0.5–5
Appendix 2—figure 4 (White noise/designed)	100	100	10,000	1	10⁻²	0.1

2.2 Robustness to nonlinear state dynamics

We introduce nonlinearity into the state update by replacing the linear update with a nonlinear function of the current state,

x (t + 1) = \frac{A}{β} \tanh (β x (t)) + B u (t) + ξ (t),

on a 4-dimensional two-mode network, which is the extended setup of Figure 2. The detailed parameters are given in Appendix 2—table 1. The $1 / β$ prefactor normalises the dynamics. The $\tanh$ function is a common choice for introducing nonlinearity (Khalil, 2002) while maintaining bounded outputs, and the parameter β controls the degree of nonlinearity (Appendix 2—figure 1a). Because the true map is now nonlinear in x, a direct comparison of $\hat{A}$ against the true A is no longer well-defined; we, therefore, evaluate the misspecified linear OLS estimator by the relative one-step prediction RMSE and the relative free-run (rollout) RMSE.

Under both metrics, the impulse condition tracks the true dynamics more accurately than passive observation across the swept range, with the gap shrinking when β grows large enough that the nonlinearity itself re-excites the 40 Hz hidden mode strongly (Appendix 2—figure 1b and c). Appendix 2—figure 1d and e show the time series for the passive and impulse conditions across the range of β, illustrating that the impulse condition strongly excites the mode dynamics. Appendix 2—figure 1f and g compare the true and predicted time series; under both error metrics, the impulse condition tracks the true dynamics more accurately than passive observation. Appendix 2—figure 1h shows the top three principal components of the data. Even when nonlinearity is present in the dynamics, the impulse condition still re-excites the hidden mode (third PC) and yields better predictions than passive observation.

Appendix 2—figure 1

Download asset Open asset

Robustness of the perturbation framework to nonlinear state dynamics.

(a) The $\tanh$ nonlinearity with the min–max envelope of $β x$ within the ordinary least squares (OLS) window shaded for $β \in {0.3, 1, 5, 10}$ . (b) Relative 1-step prediction RMSE vs. β. (c) Relative free-run (rollout) RMSE vs. β. (d) Training trajectories within the OLS estimation window for $β \in {0.3, 1, 5, 10}$ , 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 $\hat{A}, \hat{B}$ . (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).

2.3 Joint estimation of A and B

The naive estimator that assumes $B = I$ absorbs the unmodelled term $(B_{t r u e} - I) u$ into $\hat{A}$ , so its error in A grows with the mismatch level $ε = ‖ B_{t r u e} - I ‖_{F}$ . A natural remedy is to estimate A and B jointly by regressing the next state on the stacked $[X; U]$ block,

[\hat{A}, \hat{B}] = X_{+} [X_{-}; U]^{†},

which is unbiased whenever the input excites every state coordinate. We compare seven conditions on a 4-dimensional, 2-mode discrete-time system: the ground truth; passive observation; joint estimation under four perturbation patterns, with input on (i) node 1 alone, (ii) node 3 alone, (iii) nodes 1 and 3, and (iv) all four nodes; and a baseline that uses the all-nodes input but estimates A under the naive assumption $\hat{B} = I$ . The true matrix uses the same setup as Figure 2. The detailed parameters are given in Appendix 2—table 1. Every stimulated trial uses 4 impulses whose amplitudes on the active channels are drawn independently from $N (0, I)$ per pulse, so that the rank of the U-block of the regressor equals $min (n_{a c t i v e}, 4)$ . We parametrise the mismatch as $B_{t r u e} = I + ε R$ with a fixed random matrix R normalised so that $‖ R ‖_{F} = 1$ , which makes ε the scalar control parameter for the sweep.

The results show that the perturbation can still improve the estimation of A even when B is unknown, provided that A and B are recovered jointly by OLS. In the joint-estimation results of Appendix 2—figure 2a, only the columns of B corresponding to the stimulated channels are recovered: the (1, 3) pattern fills two columns, and the all-nodes pattern recovers B to the accuracy attainable when B is known. Nevertheless, A can be estimated correctly, including its hidden dynamical mode. The third column of the same panel shows that the naive estimator under the assumption $\hat{B} = I$ recovers A only approximately. Appendix 2—figure 2b shows the corresponding state trajectories. Appendix 2—figure 2c shows $‖ A - \hat{A} ‖_{F}$ as ε is swept on a log scale. Passive observation is pinned at the hidden-mode bias, and the naive $\hat{B} = I$ baseline grows with ε as the unmodelled $(B_{t r u e} - I) u_{t}$ term inflates the residual variance. Among the joint estimators, stimulating node 1 alone yields an intermediate error because node 1 is not directly coupled to the hidden mode; stimulating node 3 yields a markedly lower error because node 3 excites the hidden mode and thereby renders it identifiable; and stimulating all nodes is better still, matching the error achievable when B is known. The (1, 3) pattern lies between the node 3 and all-nodes cases. The error $‖ B - \hat{B} ‖_{F}$ grows linearly with ε for the rank-deficient stimulation patterns and the naive baseline, and matches the known-B baseline only for the all-nodes pattern under joint OLS. Although these results provide only a preliminary demonstration of the robustness of the proposed framework, joint estimation of A and B is a promising direction for future work to relax the assumption of known B.

Appendix 2—figure 2

Download asset Open asset

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 $A$ and $B$ at $ε = 20$ from a single trial. Seven conditions are compared: ground truth, passive observation, a naive baseline that uses the all-nodes input but assumes $\hat{B} = I$ , 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) $‖ A - \hat{A} ‖_{F}$ as a function of ε (log scale). (d) $‖ B - \hat{B} ‖_{F}$ as a function of ε.

2.4 Demonstration of impulse and step input

Impulse- and step-like inputs, such as those applied in TMS or tDCS, can substantially improve system identification accuracy when applied at high intensities. By perturbing the system state further from equilibrium, the resulting decay dynamics provide richer trajectory information. This richer information increases the eigenvalues of the state covariance matrix and, in turn, reduces estimation error.

The relationship between input intensity and system identification accuracy was examined through simulations using matrix A, as illustrated in Appendix 2—figure 3a. For the dynamics generated by matrix A, system identification was performed using impulse and step inputs of varying intensities. As shown in Appendix 2—figure 3b and c, the identification error decreases exponentially as the input intensity increases. This error was observed to be inversely correlated with the sum of the inverses of the eigenvalues of the state covariance matrix. In practice, the strength of perturbation input must be constrained within the limits allowed by experimental conditions. However, our findings indicate that perturbation inputs should be designed to be as strong as possible within biologically acceptable ranges.

Appendix 2—figure 3

Download asset Open asset

Variation of eigenvalues and system matrix error with impulse and step input strengths.

(a) Assumed system matrix $A$ and its eigenvalues. (b) Changes in the sum of eigenvalues, the sum of their reciprocals, and the estimation error of $A$ with respect to impulse input strength. (c) Changes in the same quantities with respect to step input strength.

2.5 White noise vs. concatenated frequency input

Appendix 2—figure 4

Download asset Open asset

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 ( $\sum_{i} (1 / μ_{i})$ ) 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 1^st–3^rd and 30^th–32^nd eigenvalues.

We demonstrate both the effectiveness and suboptimality of white noise input for neural system identification, a standard approach widely used in control engineering (Goodwin and Payne, 1977; Ljung, 1999). White noise contains a broad and nearly uniform range of frequencies, making it a convenient and practical choice for identifying system models when no prior information is available. However, as discussed in the main text, the input should ideally be designed to match the dynamical modes of the system. Thus, we consider white noise input as a practical yet potentially suboptimal method for neural system identification.

To assess the feasibility of the white noise perturbation input, we assume that a neural perturbation system is capable of generating random noise-like signals and that these input signals can be recorded simultaneously. We generated 32-channel time-series data using the matrix A shown in Appendix 2—figure 4a. Data were generated at sampling rate $f_{s} = 100$ Hz, with total length $T = 10, 000$ samples. The system matrix $A \in R^{32 \times 32}$ was constructed as the block diagonal of 16 discrete-time 2×2 rotation blocks whose eigenvalue radii were drawn from $U n i f (0.9, 0.95)$ and modal frequencies from $U n i f (1, 100)$ Hz. The system dynamics exhibit multiple modes, as shown in Appendix 2—figure 4b. The input matrix was $B = I_{32}$ with single-node stimulation applied through its first column. Process noise was i.i.d. Gaussian with per-component standard deviation $σ_{ξ} = 10^{- 2}$ . The designed input was a sinusoidal input with frequencies matching the modal frequencies of the system, which were extracted from the eigenvalues of A. Starting from white noise, the modal frequencies ${ω_{k}}$ extracted from the eigenvalues of $\hat{A}$ were used as the component frequencies of the next input (design-1), and this redesign was iterated twice more (design-2, design-3). To ensure a fair comparison, both the white noise and the designed inputs were normalized to have the same total input power. This setup allows for a controlled comparison between white noise input and frequency-optimized input.

The system identification error decreased when using frequency-designed inputs rather than white noise. Appendix 2—figure 4c shows the sum of reciprocal eigenvalues ( $\sum_{i} (1 / μ_{i})$ ) of the covariance matrix of the state vectors. This value decreased when using frequency-designed inputs rather than white noise, indicating improved accuracy in system identification. Moreover, iterative redesigns progressively reduced the sum of reciprocal eigenvalues (designs 2 and 3), leading to more accurate identification. As shown in Appendix 2—figure 4d, both white noise and designed concatenated frequency inputs excited all dynamical modes without collapsing the smaller axes. The ellipsoids obtained with the designed input had smaller volumes, indicating higher identification accuracy. Therefore, this simulation demonstrates that iterative optimal design can yield more accurate system identification than conventional white noise stimulation.

In summary, white noise input offers a simple and practical strategy for neural system identification when no prior dynamical information is available. Nevertheless, it remains suboptimal compared with concatenated frequency input, which provides higher accuracy and stability through iterative refinement. These findings also highlight the importance of input optimization in achieving precise system identification in complex neural systems.

Data availability

The current manuscript is a computational study, so no biological data have been generated. The code for the simulations and analyses presented in this paper is openly accessible at https://github.com/mikito-ogino/NeuroPerturbID (copy archived at Ogino, 2026) released under a dual MIT or BSD 3-Clause licence.

References

Conference
1. Ahmed S
2. Nozari E
(2022) On the linearizing effect of spatial averaging in large-scale populations of homogeneous nonlinear systems
2022 IEEE 61st Conference on Decision and Control (CDC).

https://doi.org/10.1109/CDC51059.2022.9993260
- Google Scholar
(2017) Effect of tDCS stimulation of motor cortex and cerebellum on EEG classification of motor imagery and sensorimotor band power
Journal of Neuroengineering and Rehabilitation 14:31.

https://doi.org/10.1186/s12984-017-0242-1
- PubMed
- Google Scholar
1. Antal A
2. Alekseichuk I
3. Bikson M
4. Brockmöller J
5. Brunoni AR
6. Chen R
7. Cohen LG
8. Dowthwaite G
9. Ellrich J
10. Flöel A
11. Fregni F
12. George MS
13. Hamilton R
14. Haueisen J
15. Herrmann CS
16. Hummel FC
17. Lefaucheur JP
18. Liebetanz D
19. Loo CK
20. McCaig CD
21. Miniussi C
22. Miranda PC
23. Moliadze V
24. Nitsche MA
25. Nowak R
26. Padberg F
27. Pascual-Leone A
28. Poppendieck W
29. Priori A
30. Rossi S
31. Rossini PM
32. Rothwell J
33. Rueger MA
34. Ruffini G
35. Schellhorn K
36. Siebner HR
37. Ugawa Y
38. Wexler A
39. Ziemann U
40. Hallett M
41. Paulus W
(2017) Low intensity transcranial electric stimulation: Safety, ethical, legal regulatory and application guidelines
Clinical Neurophysiology 128:1774–1809.

https://doi.org/10.1016/j.clinph.2017.06.001
- Google Scholar
1. Arbabi H
2. Mezić I
(2017) Ergodic theory, dynamic mode decomposition, and computation of spectral properties of the koopman operator
SIAM Journal on Applied Dynamical Systems 16:2096–2126.

https://doi.org/10.1137/17M1125236
- Google Scholar
1. Bassett DS
2. Sporns O
(2017) Network neuroscience
Nature Neuroscience 20:353–364.

https://doi.org/10.1038/nn.4502
- PubMed
- Google Scholar
(2021) Concurrent TMS-fMRI for causal network perturbation and proof of target engagement
NeuroImage 237:118093.

https://doi.org/10.1016/j.neuroimage.2021.118093
- PubMed
- Google Scholar
1. Bernal-Casas D
2. Lee HJ
3. Weitz AJ
4. Lee JH
(2017) Studying brain circuit function with dynamic causal modeling for optogenetic fMRI
Neuron 93:522–532.

https://doi.org/10.1016/j.neuron.2016.12.035
- PubMed
- Google Scholar
(2016) Optimally controlling the human connectome: the role of network topology
Scientific Reports 6:30770.

https://doi.org/10.1038/srep30770
- PubMed
- Google Scholar
1. Bikson M
2. Grossman P
3. Thomas C
4. Zannou AL
5. Jiang J
6. Adnan T
7. Mourdoukoutas AP
8. Kronberg G
9. Truong D
10. Boggio P
11. Brunoni AR
12. Charvet L
13. Fregni F
14. Fritsch B
15. Gillick B
16. Hamilton RH
17. Hampstead BM
18. Jankord R
19. Kirton A
20. Knotkova H
21. Liebetanz D
22. Liu A
23. Loo C
24. Nitsche MA
25. Reis J
26. Richardson JD
27. Rotenberg A
28. Turkeltaub PE
29. Woods AJ
(2016) Safety of transcranial direct current stimulation: evidence based update 2016
Brain Stimulation 9:641–661.

https://doi.org/10.1016/j.brs.2016.06.004
- Google Scholar
1. Bogéa Ribeiro L
2. da Silva Filho M
(2023) Systematic review on EEG analysis to diagnose and treat autism by evaluating functional connectivity and spectral power
Neuropsychiatric Disease and Treatment 19:415–424.

https://doi.org/10.2147/NDT.S394363
- PubMed
- Google Scholar
1. Boly M
2. Massimini M
3. Garrido MI
4. Gosseries O
5. Noirhomme Q
6. Laureys S
7. Soddu A
(2012) Brain connectivity in disorders of consciousness
Brain Connectivity 2:1–10.

https://doi.org/10.1089/brain.2011.0049
- PubMed
- Google Scholar
1. Braun U
2. Harneit A
3. Pergola G
4. Menara T
5. Schäfer A
6. Betzel RF
7. Zang Z
8. Schweiger JI
9. Zhang X
10. Schwarz K
11. Chen J
12. Blasi G
13. Bertolino A
14. Durstewitz D
15. Pasqualetti F
16. Schwarz E
17. Meyer-Lindenberg A
18. Bassett DS
19. Tost H
(2021) Brain network dynamics during working memory are modulated by dopamine and diminished in schizophrenia
Nature Communications 12:3478.

https://doi.org/10.1038/s41467-021-23694-9
- PubMed
- Google Scholar
1. Brunton SL
2. Brunton BW
3. Proctor JL
4. Kaiser E
5. Kutz JN
(2017) Chaos as an intermittently forced linear system
Nature Communications 8:19.

https://doi.org/10.1038/s41467-017-00030-8
- PubMed
- Google Scholar
1. Casali AG
2. Gosseries O
3. Rosanova M
4. Boly M
5. Sarasso S
6. Casali KR
7. Casarotto S
8. Bruno M-A
9. Laureys S
10. Tononi G
11. Massimini M
(2013) A theoretically based index of consciousness independent of sensory processing and behavior
Science Translational Medicine 5:198ra105.

https://doi.org/10.1126/scitranslmed.3006294
- PubMed
- Google Scholar
(2018) Time, frequency, and time-varying Granger-causality measures in neuroscience
Statistics in Medicine 37:1910–1931.

https://doi.org/10.1002/sim.7621
- PubMed
- Google Scholar
(2012) Multivariate autoregressive models with exogenous inputs for intracerebral responses to direct electrical stimulation of the human brain
Frontiers in Human Neuroscience 6:317.

https://doi.org/10.3389/fnhum.2012.00317
- Google Scholar
Book
1. Chow C
2. Dan T
3. Styner M
4. Wu G
(2024) Understanding brain dynamics through neural koopman operator with structure-function coupling: 27th international conference, marrakesh, morocco, october 6–10, 2024, proceedings, part II
In: Linguraru MG, editors. Medical Image Computing and Computer Assisted Intervention – MICCAI 2024, Vol. 15002 of Lecture Notes in Computer Science. Springer Nature. pp. 509–518.

https://doi.org/10.1007/978-3-031-72069-7_48
- Google Scholar
1. Deisseroth K
(2011) Optogenetics
Nature Methods 8:26–29.

https://doi.org/10.1038/nmeth.f.324
- PubMed
- Google Scholar
1. Deng S
2. Li J
3. Thomas Yeo BT
4. Gu S
(2022) Control theory illustrates the energy efficiency in the dynamic reconfiguration of functional connectivity
Communications Biology 5:295.

https://doi.org/10.1038/s42003-022-03196-0
- PubMed
- Google Scholar
Book
(2006) Handbook of Time Series Analysis: Recent Theoretical Developments and Applications
Wiley-VCH Verlag.

https://doi.org/10.1002/9783527609970
- Google Scholar
1. Fernández Galán R
(2008) On how network architecture determines the dominant patterns of spontaneous neural activity
PLOS ONE 3:e2148.

https://doi.org/10.1371/journal.pone.0002148
- PubMed
- Google Scholar
1. Friston KJ
(2011) Functional and effective connectivity: a review
Brain Connectivity 1:13–36.

https://doi.org/10.1089/brain.2011.0008
- PubMed
- Google Scholar
(2013) Analysing connectivity with Granger causality and dynamic causal modelling
Current Opinion in Neurobiology 23:172–178.

https://doi.org/10.1016/j.conb.2012.11.010
- PubMed
- Google Scholar
1. George MS
(2019) Whither TMS: a one-trick pony or the beginning of a neuroscientific revolution?
The American Journal of Psychiatry 176:904–910.

https://doi.org/10.1176/appi.ajp.2019.19090957
- PubMed
- Google Scholar
(2019) The roles of supervised machine learning in systems neuroscience
Progress in Neurobiology 175:126–137.

https://doi.org/10.1016/j.pneurobio.2019.01.008
- PubMed
- Google Scholar
(2003) Investigating directed cortical interactions in time-resolved fMRI data using vector autoregressive modeling and Granger causality mapping
Magnetic Resonance Imaging 21:1251–1261.

https://doi.org/10.1016/j.mri.2003.08.026
- Google Scholar
Book
1. Goodwin GC
2. Payne RL
(1977)
Dynamic System Identification. Experiment design and Data Analysis

Academic Press.
- Google Scholar
1. Green M
2. Moore JB
(1986) Persistence of excitation in linear systems
Systems & Control Letters 7:351–360.

https://doi.org/10.1016/0167-6911(86)90052-6
- Google Scholar
1. Greicius MD
2. Flores BH
3. Menon V
4. Glover GH
5. Solvason HB
6. Kenna H
7. Reiss AL
8. Schatzberg AF
(2007) Resting-state functional connectivity in major depression: abnormally increased contributions from subgenual cingulate cortex and thalamus
Biological Psychiatry 62:429–437.

https://doi.org/10.1016/j.biopsych.2006.09.020
- PubMed
- Google Scholar
1. Grimm C
2. Duss SN
3. Privitera M
4. Munn BR
5. Karalis N
6. Frässle S
7. Wilhelm M
8. Patriarchi T
9. Razansky D
10. Wenderoth N
11. Shine JM
12. Bohacek J
13. Zerbi V
(2024) Tonic and burst-like locus coeruleus stimulation distinctly shift network activity across the cortical hierarchy
Nature Neuroscience 27:2167–2177.

https://doi.org/10.1038/s41593-024-01755-8
- PubMed
- Google Scholar
(2023) A meta-analysis suggests that tACS improves cognition in healthy, aging, and psychiatric populations
Science Translational Medicine 15:eabo2044.

https://doi.org/10.1126/scitranslmed.abo2044
- PubMed
- Google Scholar
1. Gu S
2. Pasqualetti F
3. Cieslak M
4. Telesford QK
5. Yu AB
6. Kahn AE
7. Medaglia JD
8. Vettel JM
9. Miller MB
10. Grafton ST
11. Bassett DS
(2015) Controllability of structural brain networks
Nature Communications 6:8414.

https://doi.org/10.1038/ncomms9414
- PubMed
- Google Scholar
1. Hallett M
2. Di Iorio R
3. Rossini PM
4. Park JE
5. Chen R
6. Celnik P
7. Strafella AP
8. Matsumoto H
9. Ugawa Y
(2017) Contribution of transcranial magnetic stimulation to assessment of brain connectivity and networks
Clinical Neurophysiology 128:2125–2139.

https://doi.org/10.1016/j.clinph.2017.08.007
- Google Scholar
Book
1. Hamilton JD
(1994) Time Series Analysis
Princeton University Press.

https://doi.org/10.2307/j.ctv14jx6sm
- Google Scholar
1. Heinrichs JH
(2012) The promises and perils of non-invasive brain stimulation
International Journal of Law and Psychiatry 35:121–129.

https://doi.org/10.1016/j.ijlp.2011.12.006
- PubMed
- Google Scholar
1. Honey CJ
2. Sporns O
3. Cammoun L
4. Gigandet X
5. Thiran JP
6. Meuli R
7. Hagmann P
(2009) Predicting human resting-state functional connectivity from structural connectivity
PNAS 106:2035–2040.

https://doi.org/10.1073/pnas.0811168106
- PubMed
- Google Scholar
Preprint
1. Huang A
2. Ostrow M
3. Singh SH
4. Kozachkov L
5. Fiete I
6. Rajan K
(2025) InputDSA: demixing then comparing recurrent and externally driven dynamics
arXiv.

https://doi.org/10.48550/arXiv.2510.25943
- Google Scholar
(2018) Convergence properties of adaptive systems and the definition of exponential stability
SIAM Journal on Control and Optimization 56:2463–2484.

https://doi.org/10.1137/15M1047805
- PubMed
- Google Scholar
1. Kamiya S
2. Kawakita G
3. Sasai S
4. Kitazono J
5. Oizumi M
(2023) Optimal control costs of brain state transitions in linear stochastic systems
The Journal of Neuroscience 43:270–281.

https://doi.org/10.1523/JNEUROSCI.1053-22.2022
- PubMed
- Google Scholar
1. Karrer TM
2. Kim JZ
3. Stiso J
4. Kahn AE
5. Pasqualetti F
6. Habel U
7. Bassett DS
(2020) A practical guide to methodological considerations in the controllability of structural brain networks
Journal of Neural Engineering 17:026031.

https://doi.org/10.1088/1741-2552/ab6e8b
- PubMed
- Google Scholar
1. Kawakita G
2. Kamiya S
3. Sasai S
4. Kitazono J
5. Oizumi M
(2022) Quantifying brain state transition cost via Schrödinger Bridge
Network Neuroscience 6:118–134.

https://doi.org/10.1162/netn_a_00213
- PubMed
- Google Scholar
Book
1. Khalil HK
(2002)
Nonlinear Systems

Prentice-Hall.
- Google Scholar
1. Kim JZ
2. Soffer JM
3. Kahn AE
4. Vettel JM
5. Pasqualetti F
6. Bassett DS
(2018) Role of graph architecture in controlling dynamical networks with applications to neural systems
Nature Physics 14:91–98.

https://doi.org/10.1038/nphys4268
- PubMed
- Google Scholar
1. Koopman BO
(1931) Hamiltonian systems and transformation in Hilbert space
PNAS 17:315–318.

https://doi.org/10.1073/pnas.17.5.315
- Google Scholar
1. Kucyi A
2. Schrouff J
3. Bickel S
4. Foster BL
5. Shine JM
6. Parvizi J
(2018) Intracranial electrophysiology reveals reproducible intrinsic functional connectivity within human brain networks
The Journal of Neuroscience 38:4230–4242.

https://doi.org/10.1523/JNEUROSCI.0217-18.2018
- PubMed
- Google Scholar
Preprint
1. LaFosse PK
2. Zhou Z
3. O’Rawe JF
4. Friedman NG
5. Scott VM
6. Deng Y
7. Histed MH
(2023) Single-cell optogenetics reveals attenuation-by-suppression in visual cortical neurons
bioRxiv.

https://doi.org/10.1101/2023.09.13.557650
- Google Scholar
1. Lee M
2. Yoon JG
3. Lee SW
(2020) Predicting motor imagery performance from resting-state EEG using dynamic causal modeling
Frontiers in Human Neuroscience 14:321.

https://doi.org/10.3389/fnhum.2020.00321
- PubMed
- Google Scholar
(2023) Inferring causal connectivity from pairwise recordings and optogenetics
PLOS Computational Biology 19:e1011574.

https://doi.org/10.1371/journal.pcbi.1011574
- PubMed
- Google Scholar
Book
1. Ljung L
(1999)
System Identification: Theory for the User

Prentice Hall PTR.
- Google Scholar
1. Lozano AM
2. Lipsman N
3. Bergman H
4. Brown P
5. Chabardes S
6. Chang JW
7. Matthews K
8. McIntyre CC
9. Schlaepfer TE
10. Schulder M
11. Temel Y
12. Volkmann J
13. Krauss JK
(2019) Deep brain stimulation: current challenges and future directions
Nature Reviews. Neurology 15:148–160.

https://doi.org/10.1038/s41582-018-0128-2
- PubMed
- Google Scholar
1. Lu X
2. Cannon M
(2023) Robust adaptive model predictive control with persistent excitation conditions
Automatica 152:110959.

https://doi.org/10.1016/j.automatica.2023.110959
- Google Scholar
1. Mager T
2. Lopez de la Morena D
3. Senn V
4. Schlotte J
5. D´Errico A
6. Feldbauer K
7. Wrobel C
8. Jung S
9. Bodensiek K
10. Rankovic V
11. Browne L
12. Huet A
13. Jüttner J
14. Wood PG
15. Letzkus JJ
16. Moser T
17. Bamberg E
(2018) High frequency neural spiking and auditory signaling by ultrafast red-shifted optogenetics
Nature Communications 9:1750.

https://doi.org/10.1038/s41467-018-04146-3
- Google Scholar
1. Massimini M
2. Ferrarelli F
3. Esser SK
4. Riedner BA
5. Huber R
6. Murphy M
7. Peterson MJ
8. Tononi G
(2007) Triggering sleep slow waves by transcranial magnetic stimulation
PNAS 104:8496–8501.

https://doi.org/10.1073/pnas.0702495104
- Google Scholar
Conference
(2024) MiSO: Optimizing brain stimulation to create neural activity states
Advances in Neural Information Processing Systems 37.

https://doi.org/10.52202/079017-0760
- Google Scholar
(2023) TMS-evoked responses are driven by recurrent large-scale network dynamics
eLife 12:e83232.

https://doi.org/10.7554/eLife.83232
- PubMed
- Google Scholar
(2024) Controllability of functional and structural brain networks
Complexity 2024:e7402894.

https://doi.org/10.1155/2024/7402894
- Google Scholar
1. Nentwich M
2. Ai L
3. Madsen J
4. Telesford QK
5. Haufe S
6. Milham MP
7. Parra LC
(2020) Functional connectivity of EEG is subject-specific, associated with phenotype, and different from fMRI
NeuroImage 218:117001.

https://doi.org/10.1016/j.neuroimage.2020.117001
- PubMed
- Google Scholar
Book
1. Nise NS
(2020)
Control Systems Engineering

John Wiley & Sons.
- Google Scholar
1. Nozari E
2. Bertolero MA
3. Stiso J
4. Caciagli L
5. Cornblath EJ
6. He X
7. Mahadevan AS
8. Pappas GJ
9. Bassett DS
(2024) Macroscopic resting-state brain dynamics are best described by linear models
Nature Biomedical Engineering 8:68–84.

https://doi.org/10.1038/s41551-023-01117-y
- PubMed
- Google Scholar
Book
1. Ogata K
(2010)
Modern Control Engineering

Prentice Hall.
- Google Scholar
Software
1. Ogino M
(2026) NeuroPerturbID, version swh:1:rev:1dd29a2f5fbe40c8a7e22b117d5ef5da55cf0890
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:e6f83307be8aa194c7055d0780995b826723e3c2;origin=https://github.com/mikito-ogino/NeuroPerturbID;visit=swh:1:snp:414922a8e1da9e17cc136d09f2ff9a54f4ccde37;anchor=swh:1:rev:1dd29a2f5fbe40c8a7e22b117d5ef5da55cf0890
Preprint
1. Ostrow M
2. Eisen A
3. Fiete I
(2024) Delay embedding theory of neural sequence models
arXiv.

https://doi.org/10.48550/arXiv.2406.11993
- Google Scholar
1. Park AH
2. Lee SH
3. Lee C
4. Kim J
5. Lee HE
6. Paik S-B
7. Lee KJ
8. Kim D
(2016) Optogenetic mapping of functional connectivity in freely moving mice via insertable wrapping electrode array beneath the skull
ACS Nano 10:2791–2802.

https://doi.org/10.1021/acsnano.5b07889
- PubMed
- Google Scholar
1. Paulus W
(2011) Transcranial electrical stimulation (tES - tDCS; tRNS, tACS) methods
Neuropsychological Rehabilitation 21:602–617.

https://doi.org/10.1080/09602011.2011.557292
- PubMed
- Google Scholar
(2016) Application of transcranial electric stimulation (tdcs, tacs, trns): from motor-evoked potentials towards modulation of behaviour
European Psychologist 21:4–14.

https://doi.org/10.1027/1016-9040/a000242
- Google Scholar
1. Petersen SE
2. Sporns O
(2015) Brain networks and cognitive architectures
Neuron 88:207–219.

https://doi.org/10.1016/j.neuron.2015.09.027
- Google Scholar
1. Reid AT
2. Headley DB
3. Mill RD
4. Sanchez-Romero R
5. Uddin LQ
6. Marinazzo D
7. Lurie DJ
8. Valdés-Sosa PA
9. Hanson SJ
10. Biswal BB
11. Calhoun V
12. Poldrack RA
13. Cole MW
(2019) Advancing functional connectivity research from association to causation
Nature Neuroscience 22:1751–1760.

https://doi.org/10.1038/s41593-019-0510-4
- PubMed
- Google Scholar
1. Rocchi F
2. Canella C
3. Noei S
4. Gutierrez-Barragan D
5. Coletta L
6. Galbusera A
7. Stuefer A
8. Vassanelli S
9. Pasqualetti M
10. Iurilli G
11. Panzeri S
12. Gozzi A
(2022) Increased fMRI connectivity upon chemogenetic inhibition of the mouse prefrontal cortex
Nature Communications 13:1056.

https://doi.org/10.1038/s41467-022-28591-3
- PubMed
- Google Scholar
(2009) Safety, ethical considerations, and application guidelines for the use of transcranial magnetic stimulation in clinical practice and research
Clinical Neurophysiology 120:2008–2039.

https://doi.org/10.1016/j.clinph.2009.08.016
- PubMed
- Google Scholar
1. Rubinov M
2. Sporns O
(2010) Complex network measures of brain connectivity: uses and interpretations
NeuroImage 52:1059–1069.

https://doi.org/10.1016/j.neuroimage.2009.10.003
- PubMed
- Google Scholar
(2014) Real-time fMRI brain computer interfaces: self-regulation of single brain regions to networks
Biological Psychology 95:4–20.

https://doi.org/10.1016/j.biopsycho.2013.04.010
- PubMed
- Google Scholar
1. Seguin C
2. Jedynak M
3. David O
4. Mansour S
5. Sporns O
6. Zalesky A
(2023) Communication dynamics in the human connectome shape the cortex-wide propagation of direct electrical stimulation
Neuron 111:1391–1401.

https://doi.org/10.1016/j.neuron.2023.01.027
- PubMed
- Google Scholar
(2015) Granger causality analysis in neuroscience and neuroimaging
The Journal of Neuroscience 35:3293–3297.

https://doi.org/10.1523/JNEUROSCI.4399-14.2015
- PubMed
- Google Scholar
1. Shakeel A
2. Onojima T
3. Tanaka T
4. Kitajo K
(2021) Real-time implementation of EEG oscillatory phase-informed visual stimulation using a least mean square-based AR model
Journal of Personalized Medicine 11:38.

https://doi.org/10.3390/jpm11010038
- PubMed
- Google Scholar
Preprint
(2025) Quantifying state-dependent control properties of brain dynamics from perturbation responses
bioRxiv.

https://doi.org/10.1101/2025.02.18.638784
- Google Scholar
1. Shimkin N
2. Feuer A
(1987) Persistency of excitation in continuous-time systems
Systems & Control Letters 9:225–233.

https://doi.org/10.1016/0167-6911(87)90044-2
- Google Scholar
(2022) Causal mapping of human brain function
Nature Reviews Neuroscience 23:361–375.

https://doi.org/10.1038/s41583-022-00583-8
- Google Scholar
1. Siegle JH
2. Jia X
3. Durand S
4. Gale S
5. Bennett C
6. Graddis N
7. Heller G
8. Ramirez TK
9. Choi H
10. Luviano JA
11. Groblewski PA
12. Ahmed R
13. Arkhipov A
14. Bernard A
15. Billeh YN
16. Brown D
17. Buice MA
18. Cain N
19. Caldejon S
20. Casal L
21. Cho A
22. Chvilicek M
23. Cox TC
24. Dai K
25. Denman DJ
26. de Vries SEJ
27. Dietzman R
28. Esposito L
29. Farrell C
30. Feng D
31. Galbraith J
32. Garrett M
33. Gelfand EC
34. Hancock N
35. Harris JA
36. Howard R
37. Hu B
38. Hytnen R
39. Iyer R
40. Jessett E
41. Johnson K
42. Kato I
43. Kiggins J
44. Lambert S
45. Lecoq J
46. Ledochowitsch P
47. Lee JH
48. Leon A
49. Li Y
50. Liang E
51. Long F
52. Mace K
53. Melchior J
54. Millman D
55. Mollenkopf T
56. Nayan C
57. Ng L
58. Ngo K
59. Nguyen T
60. Nicovich PR
61. North K
62. Ocker GK
63. Ollerenshaw D
64. Oliver M
65. Pachitariu M
66. Perkins J
67. Reding M
68. Reid D
69. Robertson M
70. Ronellenfitch K
71. Seid S
72. Slaughterbeck C
73. Stoecklin M
74. Sullivan D
75. Sutton B
76. Swapp J
77. Thompson C
78. Turner K
79. Wakeman W
80. Whitesell JD
81. Williams D
82. Williford A
83. Young R
84. Zeng H
85. Naylor S
86. Phillips JW
87. Reid RC
88. Mihalas S
89. Olsen SR
90. Koch C
(2021) Survey of spiking in the mouse visual system reveals functional hierarchy
Nature 592:86–92.

https://doi.org/10.1038/s41586-020-03171-x
- PubMed
- Google Scholar
1. Sporns O
(2014) Contributions and challenges for network models in cognitive neuroscience
Nature Neuroscience 17:652–660.

https://doi.org/10.1038/nn.3690
- PubMed
- Google Scholar
Book
1. Sporns O
(2016) Networks of the Brain
The MIT Press.

https://doi.org/10.7551/mitpress/8476.001.0001
- Google Scholar
1. Steinberg EE
2. Janak PH
(2013) Establishing causality for dopamine in neural function and behavior with optogenetics
Brain Research 1511:46–64.

https://doi.org/10.1016/j.brainres.2012.09.036
- PubMed
- Google Scholar
(2020) Inferring causal networks of dynamical systems through transient dynamics and perturbation
Physical Review. E 102:042309.

https://doi.org/10.1103/PhysRevE.102.042309
- PubMed
- Google Scholar
1. Stevens MC
(2009) The developmental cognitive neuroscience of functional connectivity
Brain and Cognition 70:1–12.

https://doi.org/10.1016/j.bandc.2008.12.009
- PubMed
- Google Scholar
1. Suk HI
2. Wee CY
3. Lee SW
4. Shen D
(2016) State-space model with deep learning for functional dynamics estimation in resting-state fMRI
NeuroImage 129:292–307.

https://doi.org/10.1016/j.neuroimage.2016.01.005
- Google Scholar
1. Suweis S
2. Tu C
3. Rocha RP
4. Zampieri S
5. Zorzi M
6. Corbetta M
(2019) Brain controllability: Not a slam dunk yet
NeuroImage 200:552–555.

https://doi.org/10.1016/j.neuroimage.2019.07.012
- PubMed
- Google Scholar
Book
1. Takens F
(1980) Detecting strange attractors in turbulence
In: Rand D, Young LS, editors. Dynamical Systems and Turbulence. Springer. pp. 366–381.

https://doi.org/10.1007/BFb0091924
- Google Scholar
(2025) What are the optimal transcranial direct current stimulation parameters and design elements to modulate corticospinal excitability? A systematic review and longitudinal meta-analysis
Neurological Research and Practice 7:86.

https://doi.org/10.1186/s42466-025-00449-1
- PubMed
- Google Scholar
1. Tik M
2. Woletz M
3. Schuler A-L
4. Vasileiadi M
5. Cash RFH
6. Zalesky A
7. Lamm C
8. Windischberger C
(2023) Acute TMS/fMRI response explains offline TMS network effects – An interleaved TMS-fMRI study
NeuroImage 267:119833.

https://doi.org/10.1016/j.neuroimage.2022.119833
- Google Scholar
(2015) Is first-order vector autoregressive model optimal for fMRI Data?
Neural Computation 27:1857–1871.

https://doi.org/10.1162/NECO_a_00765
- PubMed
- Google Scholar
1. Tong L
2. Han S
3. Xue Y
4. Chen M
5. Chen F
6. Ke W
7. Shu Y
8. Ding N
9. Bewersdorf J
10. Zhou ZJ
11. Yuan P
12. Grutzendler J
(2023) Single cell in vivo optogenetic stimulation by two-photon excitation fluorescence transfer
iScience 26:107857.

https://doi.org/10.1016/j.isci.2023.107857
- PubMed
- Google Scholar
1. Tseng SY
2. Chen RC
3. Chong FC
4. Kuo TS
(1995) Evaluation of parametric methods in EEG signal analysis
Medical Engineering & Physics 17:71–78.

https://doi.org/10.1016/1350-4533(95)90380-T
- Google Scholar
1. Tu C
2. Rocha RP
3. Corbetta M
4. Zampieri S
5. Zorzi M
6. Suweis S
(2018) Warnings and caveats in brain controllability
NeuroImage 176:83–91.

https://doi.org/10.1016/j.neuroimage.2018.04.010
- Google Scholar
1. Wagenmaker A
2. Mi L
3. Rozsa M
4. Bull MS
5. Svoboda K
6. Daie K
7. Golub MD
8. Jamieson K
(2024)
Active learning of neural population dynamics using two-photon holographic optogenetics

Advances in Neural Information Processing Systems 37:31659–31687.
- PubMed
- Google Scholar
(2006) Functional connectivity networks are disrupted in left temporal lobe epilepsy
Annals of Neurology 59:335–343.

https://doi.org/10.1002/ana.20733
- PubMed
- Google Scholar
1. Yan Y
2. Murphy TH
(2024) Decoding state-dependent cortical-cerebellar cellular functional connectivity in the mouse brain
Cell Reports 43:114348.

https://doi.org/10.1016/j.celrep.2024.114348
- PubMed
- Google Scholar

Article and author information

Author details

Mikito Ogino

Graduate School of Arts and Sciences, The University of Tokyo, Meguro-ku, Japan

Contribution
Conceptualization, Software, Formal analysis, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing

For correspondence
ogino-m@g.ecc.u-tokyo.ac.jp

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-1089-4469
Daiki Sekizawa
1. Graduate School of Arts and Sciences, The University of Tokyo, Meguro-ku, Japan
2. RIKEN Center for Brain Science, RIKEN, Wako, Japan
Contribution
Conceptualization, Formal analysis, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0009-0004-0196-2612
Jun Kitazono

Graduate School of Data Science, Yokohama City University, Yokohama, Japan

Contribution
Conceptualization, Validation, Methodology, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-6701-3947
Masafumi Oizumi

Graduate School of Arts and Sciences, The University of Tokyo, Meguro-ku, Japan

Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Writing – original draft, Project administration, Writing – review and editing

For correspondence
c-oizumi@g.ecc.u-tokyo.ac.jp

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-8802-2607

Funding

Japan Science and Technology Agency (Moonshot R&D Grant Number JPMJMS2012)

Masafumi Oizumi

Japan Society for the Promotion of Science (KAKENHI Grant Number 24K20462)

Mikito Ogino

Japan Society for the Promotion of Science (KAKENHI Grant Number 23KJ0799)

Daiki Sekizawa

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We are grateful to Yumi Shikauchi, Shunsuke Kamiya, and Daiki Kiyooka for their insightful feedback and valuable discussions, which significantly contributed to the development of this research. This work was supported by JST Moonshot R&D Grant Number JPMJMS2012, JSPS KAKENHI Grant Number 24K20462, and JSPS KAKENHI Grant Number 23KJ0799.

Version history

Preprint posted: November 19, 2025
Sent for peer review: December 15, 2025
Reviewed Preprint version 1: March 9, 2026
Reviewed Preprint version 2: September 9, 2026
Version of Record published: October 13, 2026

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.110030. This DOI represents all versions, and will always resolve to the latest one.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.