A new perspective on brain stimulation interventions: Optimal stochastic tracking control of brain network dynamics

In this section, we first demonstrate the brain dynamics considering neuronal noise compared to the network model without considering noise, as illustrated in Fig. 6. We use the data of a representative participant with ID 2121 as an example, where Fig. 6.a and Fig. 6.b correspond to the results for 100-ROI and 400-ROI systems, respectively. It is evident that the simulated dynamics with neuronal noise more closely approximate the actual dynamics measured by fMRI, which can be also supported by the statistical result of KL divergence. In contrast, the dynamics without considering noise tends to converge to zero over time, due to the eigenvalues of the coupling matrix $\mathbf{A}$ being all less than zero, leading to the convergence of the state $\mathbf{x}$.

Although the results of optimal state approaching control and optimal stochastic tracking control cannot be directly compared, we can quantify their relationships with the intrinsic characteristics of the brain system to explore any correlated features. For optimal state approaching control, we used clustering methods to select four representative discrete brain states as target states for both the 100-ROI and 400-ROI brain systems. We then used the peak states of each individual participant as the original states and pushed them to the four target states. Fig. 7.a and d provide visualizations of the spatial distributions of the four target states, while Fig. 7.b and e quantify the average distribution of these states within each resting-state system. It can be observed that the four states exhibit distinct spatial distribution characteristics, and the four states in the 100-ROI and 400-ROI brain systems can be matched to each other. Specifically, state 1 of 100-ROI system corresponds to state 4 of 400-ROI system, state 2 of 100-ROI system corresponds to state 1 of 400-ROI system, state 3 of 100-ROI system corresponds to state 2 of 400-ROI system, and state 4 of 100-ROI system corresponds to state 3 of 400-ROI system. We retained the original order of the clustering results for visualization, which does not substantially affect the results.

Fig. 7.c and f illustrate the relationship between the energy of optimal state approaching control and other features (tracking energy, average controllability, modal controllability, and the value of state itself) across five different groups. Fig. 7.c presents the results for the 100-ROI brain network system, showing that the control energy required to approach each target state (1 to 4) is predominantly related to the state itself. Similarly, for the 400-ROI brain network system (Fig. 7.f), despite identifying more statistically significant correlations between the control energy and three features beyond and state value, most significant correlations are concentrated between the energy of optimal state approaching control and the target state itself. These results indicate that the energy of optimal state approaching control is highly dependent on the target state itself, rather than the intrinsic characteristics of the brain network system.

Using fMRI data from stroke and post-stroke aphasia cases as examples, we attempted to drive the disease dynamics towards healthy dynamics. We here used two open source fMRI datasets (Aphasia Recovery Cohort (ARC) Dataset, and Resting-state for Older Adults Dataset), which are described below.

Dataset 1: Resting-state fMRI data were obtained from the University of South Carolina and the Medical University of South Carolina [61]. The Institutional Review Board (IRB) at the University of South Carolina determined that the study was exempt from further ethical review. The final dataset comprised 148 post-stroke aphasia (including four aphasia types: Anomic, Broca, Conduction, Global) patients and 26 normal stroke patients, after excluding those with head motion greater than 3mm. Participants were instructed to remain still, keep their eyes open, and stay awake during the scan. All fMRI images were collected using a Siemens Trio 3.0T scanner with a 16-channel head coil. For each participant, the fMRI scan closest to the time of stroke onset was selected for further analysis.

Dataset 2: For finding target states/dynamics for patients, a target group of 28 healthy old participants underwent whole-brain imaging at the Gateway MRI Center at UNCG, utilizing a Siemens 3.0T Tim Trio MRI scanner with a 16-channel head coil, was selected for analysis. The study received approval from the IRB at the University of North Carolina at Greensboro [62]. Resting-state fMRI data were acquired following anatomical scans, with 300 measurements collected during the 10-minute scanning sessions. Participants were instructed to remain still, awake, and with their eyes open throughout the procedure. Functional imaging was performed using echo-planar imaging sequences sensitive to blood oxygen level-dependent (BOLD) contrast, with a slice thickness of 4.0 mm, 32 slices, a TE of 30 ms, a TR of 2000 ms, and a flip angle of 70°. The scans covered the entire cerebral cortex and portions of the cerebellum. The primary objective of this dataset was to explore how connectome-based modeling of intrinsic functional connectivity within the default mode network could predict memory discrimination differences between young and older adults. Statistical information of participants for population number, age, and Western Aphasia Battery (WAB) score is presented in Table 1.

For the fMRI raw scan data in the above two datasets, we performed preprocessing operations. The first ten volumes of each scan were excluded to ensure magnetic stability and to optimize the generation of a consistent BOLD activity signal. The preprocessing pipeline was implemented utilizing the CONN toolbox [2]. Subsequently, the fMRI data underwent the following procedures:

1. Correction for slice timing;

2. Realignment to correct for head motion;

3. Normalization to a 3 mm MNI space using an EPI template provided by the CONN software package;

4. Application of linear temporal detrending to remove the linear drift of the BOLD signal;

5. Covariate regression, including the white matter signal, cerebrospinal fluid (CSF) signal, and head motion signal, utilising the Fristion-24 Parameters strategy;

6. Temporal filtering with a bandwidth of 0.01-0.1 Hz.

We utilized the 100-ROI and 400-ROI parcellation templates from Schaefer2018 parcellation [21] at different scales to parcel the cortex, and for further estimation of brain network systems. We employed a regression-out approach (one of the most used harmonization techniques [63, 64]) for preprocessed data to eliminate any potential factors that could influence the results, including sites, equipments, acquisition parameters, and the age and gender of individuals.

NCT requires the establishment of brain dynamic model, along with equations describing the connectivity between brain regions. In this context, we employ the stochastic linear time-invariant model:

where $\mathbf{x}(t)$ is an $N\times 1$ vector that represents the brain state (i.e. the BOLD time series for fMRI measurement) at time $t$, and $N$ is the number of brain regions. In the network adjacency matrix (or coupled matrix) $\mathbf{A}$, each $a_{ij}$-th element gives the quantitative anisotropy between region $i$ and region $j$. $\mathbf{u}(t)$ is the input which is constant in time.

The coupled matrix $\mathbf{A}$ in Eq. 1 determines the brain network dynamics, which was described by a multivariate Ornstein-Uhlenbeck process [65, 66, 67]:

where the time constant $\tau_{\mathbf{x}}$ abstracts the intrinsic dynamics of each ROI, each activity ${x}_{i}(t)$ of node $i$ decays exponentially according to the time constant $\tau_{\mathbf{x}}$ and evolves depending on the activity of other populations. The network effective connectivity (EC) between different ROIs, embodied by the matrix $\mathbf{C}$, whose topological skeleton is determined by structural data. The fluctuating inputs ${w}_{i}(t)$ are independent and correspond to a diagonal variance matrix $\mathbf{\Sigma}$. In matrix notation, it is:

All ${x}_{i}(t)$ have zero mean, of which the spatiotemporal covariances are denoted by ${Q}_{\tau}^{ij}$, where $\tau\in\{0,1\}$ is the time lag, and they satisfy the following consistency equations:

where $e$ denotes the matrix exponential, superscript $\top$ denotes the matrix transpose. $\mathbf{J}$ is the Jacobian of the dynamical system, can be considered as a appropriate estimation for $\mathbf{A}$, which depends on the time constant $\tau_{x}$ and the network EC:

where $\delta_{ji}$ is the Kronecker delta.

Most studies treat the networks extracted from structural imaging as the $\mathbf{A}$ matrix in network control, however, it may not necessarily correspond to the measured brain network dynamics. For obtaining a fine estimation of $\mathbf{J}$ in Eq. 6 for brain network dynamics, i.e. the estimation of $\mathbf{A}$, we tune the above model in Sec. 4.2.1 to make its covariance matrices $\mathbf{Q}_{0}$ and $\mathbf{Q}_{1}$ reproduce the empirical $\mathbf{\hat{Q}}_{0}$ and $\mathbf{\hat{Q}}_{1}$ according to previous studies [66, 67], which can be calculated as follows: For a session of $T$ time points, we denote the BOLD time series by ${s}_{i}(t)$ for each region $1\leq i\leq N$ with time indexed by $1\leq t\leq T$. The mean signal is $\mathbf{\bar{s}}_{i}=\frac{1}{T}\sum_{t}{s}_{i}(t)$ for all $i$. Following previous studies [66, 67], two BOLD covariance matrices without and with time lag are:

Pearson correlation ${K_{ij}}$ can be calculated from the covariances in Eq. 7 and 8:

The iterative optimization procedure for $\mathbf{C}$ in Eq. 6 is related to a concept of natural gradient descent that accounts for the non-linearity of the mapping between $\mathbf{J}$ and the matrix pair $\mathbf{Q}_{0}$ and $\mathbf{Q}_{1}$:

where $\tau_{TR}$ is the time resolution (TR) value of fMRI in seconds.

The difference matrices $\Delta\mathbf{Q}_{0}=\mathbf{\hat{Q}}_{0}-\mathbf{Q}_{0}$ and $\Delta\mathbf{Q}_{1}=\mathbf{\hat{Q}}_{1}-\mathbf{Q}_{1}$ determine the model error:

The Jacobian update is applied to decrease the model error $E$ at each optimization step:

A update that is more robust to empirical noise is:

Optimal state approaching control requires the coupled matrix $\mathbf{A}$ to be symmetric. To ensure the symmetry of the coupled matrix, we modified the Jacobian update $\Delta\mathbf{J}$ as follows:

From the Jacobian update $\Delta\mathbf{J}$, the connectivity update was obtained:

Non-negativity is imposed on the EC values during the optimization. To take properly the effect of cross-correlated inputs into account, the $\mathbf{\Sigma}$ update is:

Finally, the time constant $\tau_{x}$ can also be tuned as:

where $\lambda_{\text{max}}$ is the maximum negative real part of the eigenvalues of $\mathbf{J}$. Repeating the parameter updates, the best fit of matrix $\mathbf{J}$ corresponds to the minimum model error $E$.

Our long-term goal is to use the model described above to predict optimal stimulation parameters, i.e., network control. In the network dynamics control model, we first consider two important factors: one is the stochastic fluctuation in the brain dynamics, which may represent spontaneous activity or noise in brain networks; the other is the optimal control from one network dynamics to a target dynamics, which can be seen as a tracking problem in the control domain, rather than pushing one discrete brain state to another. Given the linear stochastic system state equation like Eq. 1 with input:

where the input matrix $\mathbf{B}$ was defined to allow input dominated by the stimulation region of interest (ROI). The output equation is:

where $\mathbf{z}(t)$ is a $N$-dimensional output vector. In brain network systems, the matrix $\mathbf{D}$ can be considered as an identity matrix, implying that $\mathbf{z}(t)=\mathbf{x}(t)$.

Let the $N$-dimensional target vector $\mathbf{z}_{r}(t)$ be the output of the target dynamic model under the influence of stochastic noise. The state equation of the target model is known to be:

where $\mathbf{w}_{r}(t)$ is stochastic noise with zero mean and variance matrix $\mathbf{Q}_{r}$. The output equation of the target model is:

where $\mathbf{D}_{r}$ is an identity matrix, implying that $\mathbf{z}_{r}(t)=\mathbf{x}_{r}(t)$.

To take an initial step toward the optimal control goal, we seek to estimate the optimal error and energy required to reach the target brain dynamics, and therefore define the following objective: Determine the linear optimal feedback control $\mathbf{u}^{*}(t)$ such that the system output $\mathbf{z}(t)$ tracks the target dynamics $\mathbf{z}_{r}(t)$ and minimizes the following stochastic performance index:

where the final time $t_{f}$ is fixed, and $\mathbf{Q}$ and $\mathbf{R}$ are symmetric non-negative definite and symmetric positive definite matrices, respectively.

For this stochastic output feedback tracking system problem, if the system $\{\mathbf{A},\mathbf{D}\}$ is completely observable and the target model $\{\mathbf{A}_{r},\mathbf{D}_{r}\}$ is completely observable, then the optimal control is:

where

where $\mathbf{P}_{11}(t)$ and $\mathbf{P}_{12}(t)$ satisfy the following matrix differential equations and their boundary conditions:

The proof and details can be seen in Supplementary Materials.

Here, we introduce several definitions, where Kullback-Leibler (KL) divergence and optimal control energy are employed to assess control effect, while average controllability and modal controllability are used to quantify the intrinsic controllability of brain network systems.

Due to the involvement of stochastic processes in optimal stochastic tracking control, the control effect cannot be measured using Euclidean distances between discrete brain states. Instead, it should be assessed using the KL divergence to quantify the similarity between the statistic of dynamics before/after control and target dynamics. In the context of measuring the similarity between random variables, the KL divergence is a widely used metric in information theory and statistics, and its definition is given as follows:

The KL divergence is always non-negative and is zero if and only if $P$ and $Q$ are identical almost everywhere. However, it is important to note that the KL divergence is not a true metric as it is not symmetric, i.e., $D_{\text{KL}}(P\|Q)\neq D_{\text{KL}}(Q\|P)$, and it does not satisfy the triangle inequality. For each node in brain networks, the KL divergence from the dynamics before control to the target dynamics, and the KL divergence from the dynamics after control to the target dynamics were estimated and compared. We disregarded the estimation of the control effectiveness for optimal state approaching control (which can be directly estimated using methods such as Euclidean distance), because its results cannot be directly compared with those of optimal stochastic tracking control.

In the context of optimal stochastic tracking control, we employed a time step of 1 $s$ and 1000 time points for the simulation and calculation of control energy to approximate the resolution used in fMRI sampling. In optimal state approaching control, we similarly used 1000 time points but with a time step of 0.001 $s$ to define the control horizon, as this optimal control method is not suitable for a broad horizon. Consequently, the magnitudes of control energy obtained from the two methods cannot be directly compared, and only their spatial distributions can be produced to correlation analysis.

It is noteworthy that optimal state approaching control aims to drive one brain state to closely approach another target brain state within a limited short time frame. In contrast, optimal stochastic tracking control aims to track an entire brain dynamics towards another target brain dynamics. This highlights the advantage of optimal stochastic tracking control not only in considering the neuronal noise, but also in achieving its intended goal and aligns with the long-term rehabilitation objectives for conditions such as stroke and aphasia.

We quantified the relationship between control energy and the intrinsic controllability of brain network systems. The two types of controllability [43] include:

For optimal stochastic tracking control, the target dynamics is the averaged fMRI time series of healthy participants, and the original dynamics is the fMRI time series of each stroke and aphasia patients. For optimal state approaching control, since the state represents a state of a specific moment in time, we employed a clustering method for state selection, following previous research [47].

Target brain states of healthy participants: To define the brain states of healthy participants as target brain states, we referenced the selection method for brain states described in [47]. Initially, we detected the peaks of the root mean square (RMS) of whole-brain activity amplitude (N = 100 or 400 regions/nodes defined by Schaefer2018 parcellation). We defined a peak as a frame where the RMS exceeds the values of the preceding and following frames. Subsequently, we aggregated the corresponding patterns from participants and calculated the Lin’s concordance among all patterns. We then applied a variant of modularity maximization to cluster the $Npeak\times Npeak$ matrix, which resulted in 10 clusters. Most clusters were relatively small and/or contained peaks from only a few participants. We defined a brain state as the cluster center that represents peak activations from at least 50% of the participants. Ultimately, we identified four distinct brain states for healthy participants. This approach clusters BOLD time series into recurrent states and is similar to established methods for extracting co-activation patterns (CAPs) from BOLD data.

Original brain states of stoke and post-stroke aphasia patients: To define the brain states of patients as original brain states, and to ensure precision and alignment with personalized treatment applications, we calculated an individual brain state for each stroke and aphasia patient. We selected the peaks of the root mean square (RMS) of whole-brain activity amplitude as the brain state, specifically utilizing the z-scored mean of resting-state fMRI (rs-fMRI) signal amplitudes over time for each brain region [69].

For the traditional optimal brain network control in which the target is to approach the original state to a targeted state (refer to as optimal state approaching control), given the linear deterministic model:

the cost function is:

subject to

The optimal input can be estimated as follows:

where $\mathbf{u}_{k}^{*}$ can be obtained through solving corresponding differential equations (see details in Supplementary Materials).