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 [1, 2, 3, 4, 5]. This interaction—called neuronal connectivity [6, 7]—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 [8, 9, 10, 11]. The functional connectivity is estimated through recording techniques such as neural spike recording techniques [12, 13], electrocorticography (ECoG) [14, 15], functional magnetic resonance imaging (fMRI) [16, 17] and electroencephalography (EEG) [18, 19]. While various methods have been proposed and used to estimate these connections (see for example [2, 20] for a comprehensive review), one typical method is based on dynamic modeling, among which is a simple but widely used linear auto-regressive model (Fig. 1a) [21, 22, 23, 24, 25]. Connectivity is statistically estimated from the time-series data of neural activity (Fig. 1b) by fitting the model parameters (Fig. 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 [8, 26, 27, 28, 29, 30].

Fig. 1:

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 Fig. 1a, the true connectivity matrix governs the generation of time-series signals (Fig. 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 (Fig. 1c), because some modes quickly decay and vanish from the observable data (see also Fig. 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.

Fig. 2:

To overcome these limitations, the application of perturbations has emerged as a powerful approach for improving the estimation of connectivity in neuroscience [31, 32, 29, 33]. 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 [32]. Recent studies using stimulation techniques such as optogenetics [34, 35] and transcranial magnetic stimulation (TMS) [36, 29, 37] 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 [31, 29]. 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 [38, 39, 29, 40, 41]. 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 [42, 43, 44, 45, 46, 47], and treat external perturbations as control inputs to design properties of neural stimulation (Fig. 1d). If the optimal perturbation input can be systematically designed, it becomes possible to steer the neural system toward states that are maximally informative (Fig. 1e), thereby enhancing the accuracy of the inferred connectivity (Fig. 1f). While recent studies have begun to develop algorithmic approaches to active stimulation design for specific experimental platforms [48, 49], 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 [31, 36, 50, 47] and control of neural states [42, 51, 44, 45]. 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 dynamics with external inputs and stochastic system noise [52, 42, 53, 54, 43]. We assume that the brain state evolves according to the following linear dynamics:

where x(t) ∈ ℝⁿ 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 ∈ ℝ^n×n characterizes intrinsic interactions reflecting functional connectivity. The input matrix B ∈ ℝ^n×m specifies how external inputs u(t) ∈ ℝ^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) ℝⁿ represents stochastic fluctuations (e.g., synaptic noise or unobserved inputs), modeled as i.i.d. zero-mean Gaussian noise across time with covariance matrix Σ_ξ = 𝔼 [ξ(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 Refs. [55, 56].

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

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):

where ∥ · ∥_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 Supplementary Material B.2 and B.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:

Here, A and ξ(t) are identical to those defined in Eq. 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:

where

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 [57], 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 [58], 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 [59], 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 [60, 61].

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 [62], 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 [63], 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 [64, 65], 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 inputs and impulse as typically observed in tACS, tDCS, and TMS—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 Fig. 2a. Neural data are generated based on the connectivity matrix, yielding the time series shown in Fig. 2b. These data are not influenced by external input (passive state), and thus the oscillations of x₃ and x₄ 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 10Hz dynamics, while the third principal component, associated with the 40Hz mode, contributes negligibly to the total variance (Fig. 2c). Consequently, the matrix estimated by OLS deviates substantially from the true connectivity matrix (Fig. 2d), particularly showing errors in the lower-right components, as highlighted in Fig. 2e. This failure arises because the true matrix A contains two distinct dynamical modes: one is a continuous oscillatory mode at 10Hz, and the other is a damped mode at 40Hz. 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. Figures 2f–2i show the application of an impulse input u to the system. This perturbation input primarily affects the nodes corresponding to x₃ and x₄, 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 Figs. 2j–2m. 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 [66, 67, 68, 69]. 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 Supplementary Material B.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 → ∞), the estimation error of A asymptotically converges to the following expression [70]:

where with 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 the Supplementary Material A.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 , 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 Eq. 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

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 Fig. 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 (Fig. 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 (Fig. 3b). As shown in Fig. 3c, when the neural dynamics is dominated by a single eigenvalue µ₁, the other eigenvalues µ₂ and µ₃ remain small, resulting in an elongated reciprocal ellipsoid (Fig. 3d) and a large estimation error. When perturbation inputs that excite the directions associated with µ₂ and µ₃ are applied, the ellipsoid expands and becomes closer to a sphere (Fig. 3e). As a result, the estimation error decreases in all directions, as illustrated in Fig. 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.

Fig. 3:

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 Eq. 1 can be obtained by iterating the state transition matrix. Specifically,

We assume a cosine input , 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.

where each eigenvalue of A is written in polar form as , 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 Eq. 5. The vector v_d denotes the right eigenvector of A associated with λ_d, and denotes the corresponding left eigenvector. The detailed derivation is provided in the Supplementary Material A.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 Eq. 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:

For step input, the state vector is described by,

where α and β are intensity of the impulse and step inputs. The detailed derivation is provided in the Supplementary Material A.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, α² and β². However, estimation accuracy is not determined by input intensity alone. Eqs.14 and 15 show that the resulting dynamics x_diff(t) are critically dependent on the term Bu₀. This term represents how the input vector u₀ (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₀ is the key design parameter that determines which dynamical modes are excited. An improperly chosen u₀ 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₀ 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 the Supplementary Material B.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 ∑µ_i, the sum of reciprocal eigenvalues 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 Eq. 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 Figs. 4c and 4d. Both the major and minor axes of the ellipsoids in Fig. 4c are extended by the sinusoidal inputs, especially by the 10 Hz input. In contrast, the major and minor axes of the ellipsoids in Fig. 4d are substantially shortened by the 10 Hz sinusoidal input.

Fig. 4:

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 Fig. 5a, along with simulated time series and the corresponding matrix estimation. Because the system matrix contains two distinct pairs of eigenvalues (10 Hz 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 Hz and 20 Hz components; inputs with a single frequency yield larger values (Fig. 5b). To further evaluate the effect of input design, we compared the two-frequency input with a flat-spectrum input (Fig. 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 Figs. 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 Fig. 5e, leading to larger estimation errors. Small eigenvalues disproportionately increase the sum of reciprocal eigenvalues, as seen in Figs. 5f and 5g. 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).

Fig. 5:

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 Figs. 6a and 6b. 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. Fig. 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.

Fig. 6:

Stimulating Node 7 (hub node) minimizes the total estimation error of A (Fig. 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 µ₅–µ₇ are small and their reciprocals 1/µ₅–1/µ₇ are large, producing the elongated ellipsoid seen in Fig. 6d. When Node 1 receives an impulse (Fig. 6e), the eigenvalue associated with the 15 Hz mode grows and the 1/µ₅ and 1/µ₆ axes shrink; however, the 1/µ₇ axis—which correspond to directions that Node 1 alone does not excite— remains large. Stimulating Node 7 (Fig. 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. Figures 6h and 6i 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 (described in a later section) and joint estimation of A and B (Supplementary Material B.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 [31, 36, 50, 47]. 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 (Fig. 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 [50, 47]. 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 . The resulting time-series data are shown in Fig. 7b. Using this data, we estimated the underlying dynamical system via a control-based identification approach presented in Eq. 4, which corresponds to an estimation of functional connectivity.

Fig. 7:

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 (Fig. 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 Fig. 7d. ROC curve analysis further confirms that perturbation substantially improves discriminability across all tasks and outperforms the random perturbation condition (Fig. 7e). These results can be explained by Eq. 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 Fig. 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 [42, 51]. 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 (Fig. 8a). Moreover, use of the perturbation approach to derive the model allows the controllability Gramian [42, 47], as well as the estimation of control costs for each network area, to be more reliably assessed (Fig. 8b and 8c). 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.

Fig. 8:

We designed a network as shown in Fig. 9a to clearly show the differences in system identification between passive and perturbations 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³). 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 to 25 while keeping all other nodes at 0 without any movement.

Fig. 9:

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 (Fig. 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 (Fig. 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 (Fig. 9d). Here, we adopt average controllability as an example measure of control cost [42, 46]. The key point here is that, since the controllability Gramian involves multiple products of A (Eq. 21), even small errors of 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 [71, 72]. As shown in Figs. 5 and B.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 (Fig. 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 Â, and the target node of the perturbation is determined through numerical optimization that minimizes . Table B.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.

Fig. 10:

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

1. Record spontaneous activity as the passive state and estimate A_passive.

2. Based on the estimated matrix, the optimal frequency and target node of the perturbation input are determined through numerical optimization following Eq. 9.

3. The designed perturbation u_design1 is applied to the system, and a more precise matrix A_design1 is identified.

4. 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 Fig. 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 . The error from the second iteration (2nd design) is smaller than that from the first (1st 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 outperformed it 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 [9, 73, 42]. Linear models are particularly effective for providing accurate approximations of nonlinear systems within a specific operating range [74]. 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 [75, 76] and system identification methods leveraging machine learning techniques [77, 78]. 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 [79, 80, 81]. 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 [82, 17]. For stimulus intensity, our proposed framework suggested that a stronger stimulus improves system identification; however, an excessively strong stimulus poses safety risks [83]. 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 [83, 84]. Furthermore, the selection of stimulation sites is often dictated by experimental accessibility or ethical considerations [85, 86]. 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 [87], to reconstruct hidden dynamics from partial observations [88, 89, 90, 91]. 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 an essential in 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 [47]. 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 Fig. 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) (Eq. 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 the Supplementary Material B.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

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 [42, 51] with some limitations [92, 93]. 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₀ to the terminal state x_T:

under

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

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

where W is the controllability Gramian, defined as:

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 [42, 47].

Average Controllability for Control Cost

While the optimal control problem (Eq. 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 [42, 46]. It is defined as the trace of the controllability Gramian:

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.

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.

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.

Additional files

Supplementary Material

Additional information

Funding

MEXT | Japan Science and Technology Agency (JST)

https://doi.org/10.52926/jpmjms2012

• Masafumi Oizumi

Japan Society for the Promotion of Science (JSPS) (24K20462)

• Mikito Ogino

Japan Society for the Promotion of Science (JSPS) (23KJ0799)

• Daiki Sekizawa