Reverse engineering the brain input: Network control theory to identify cognitive task-related control nodes
Zhichao Liang Yinuo Zhang Jushen Wu Quanying Liu


*[1]
    0pt[0pt][0pt]















































Reverse engineering the brain input: Network control theory to identify cognitive task-related control nodes





This work was funded in part by the National Key R&D Program of China (2021YFF1200804), Shenzhen Science and Technology Innovation Committee (2022410129, KCXFZ2020122117340001).


    

Zhichao Liang1,
Yinuo Zhang1,
Jushen Wu1, and
Quanying Liu 1^*
1Department of Biomedical Engineering, Southern University of Science and Technology, Shenzhen, 518055, China

* Corresponding Author: Quanying Liu   Email: liuqy@sustech.edu.cn
    
=====================================================================================================================================================================================================================================================================================================



The human brain receives complex inputs when performing cognitive tasks, which range from external inputs via the senses to internal inputs from other brain regions. However, the explicit inputs to the brain during a cognitive task remain unclear.


Here, we present an input identification framework for reverse engineering the control nodes and the corresponding inputs to the brain. The framework is verified with synthetic data generated by a predefined linear system, indicating it can robustly reconstruct data and recover the inputs.
Then we apply the framework to the real motor-task fMRI data from 200 human subjects.
Our results show that the model with sparse inputs can reconstruct neural dynamics in motor tasks (EV=0.779) and the identified 28 control nodes largely overlap with the motor system. 
Underpinned by network control theory, our framework offers a general tool for understanding brain inputs.



Control nodes identification, Inputs identification, Brain inputs, Brain network control, fMRI




§ INTRODUCTION

When performing a complex cognitive task, multiple brain regions are engaged in sequence, for gathering sensory inputs, reasoning, making decisions, and executing actions<cit.>.
Our brain receives complex task inputs and hierarchically coordinates these task-related regions, from low-level sensory regions to high-level associative regions<cit.>. 
Uncovering which brain areas receive task-relevant inputs and how these inputs are sequentially transmitted to other brain areas are core scientific questions in neuroscience. 
These motivations for exploring task inputs to the brain system pose new technical challenges.


In neuroscience, task-related brain regions are typically identified through statistical analysis of functional magnetic resonance imaging (fMRI)<cit.>, including univariate analysis (e.g., general linear model (GLM) and t-test) and multivariate analysis (e.g., repeated measures ANOVA). Specifically, GLM regresses each voxel of fMRI data with task-related variables (e.g., task onset, reaction time, or behavioral performance). The regions that demonstrate statistical significance are identified as task-related brain areas. GLM has helped to identify the working memory-related brain regions <cit.> and the Chinese language-related brain region <cit.>. 
Recently, machine learning has been applied to uncover the brain regions and activity patterns across tasks and stimuli, such as multi-voxel pattern analysis (MVPA) <cit.> and graph convolutional networks <cit.>. Ye et al. proposed a Spatial Temporal-pyramid Graph Convolutional Network (STpGCN) to decode 23 cognitive tasks and a model-agnostic explanation tool (BrainNetX) to identify task-related brain regions <cit.>. 
Although these machine learning methods succeed in identifying task-relevant brain regions, they fall short in extracting the explicit time courses of inputs to the brain.
The ultimate goal is the reverse engineering of brain inputs, which involves uncovering the relationship between endogenous or exogenous stimuli and neural responses in complex cognitive tasks. 

Network control theory offers versatile tools for investigating the brain system. For instance, it can model input-output relationships in the brain with precision <cit.> and elucidate the brain state transitions driven by minimal energy<cit.>.
Previous studies have shown that the brain network can be controlled by specific regions<cit.>. This finding suggests that altering brain states does not require stimulation inputs to all nodes. Nevertheless, in practical scenarios, considering the large amount of energy input required, changing the brain states by input to a single brain node is unrealistic. It is widely recognized that multiple brain regions receive inputs during cognitive tasks. Adopting the terminology from control theory, we refer to these brain regions as control nodes. Identifying control nodes and determining the input sequences of these nodes represent two key questions in this work.



Here, we present an input identification framework for the reverse engineering of brain inputs during cognitive tasks. We formulate this as a joint optimization problem: i) identifying critical control nodes and ii) inferring their corresponding control sequences. 
The objective function aims to minimize the reconstruction error of the neural dynamics using a sparse control set and smooth inputs. To solve this objective function, we propose a gradient-based augmented Lagrangian method, taking into account the consideration of the relaxation of boolean variables. 
The results from synthetic data validate the feasibility of inferring the control nodes and the input sequence. We apply our framework to the task fMRI data, identifying the motor-related brain regions and the explicit input sequence during motor tasks. 
Our framework sheds light on the underlying principles linking control nodes and the topological properties of the brain. 




§ THE INPUT IDENTIFICATION FRAMEWORK


We present a control theory-driven framework designed to identify brain inputs during cognitive tasks. The framework consists of a dynamic model of the brain system, an objective function with associated constraints, and an optimization algorithm.
It is important to note that, although we implement a linear dynamic model in this study, this framework can be extended to any input-output models (e.g., SINDy <cit.> and Koopman <cit.>).




 §.§ Brain network dynamical model



Here, we model the brain system as a discrete-time linear time-invariant network with input (Fig. <ref>). The brain network model is formulated as Eq (<ref>).

    𝐱(t+1)=𝐀𝐱(t)+𝐁𝐮(t),

where vector 𝐱(t) ∈ℝ^N denotes the state of the brain network with N regions at time t. Within our brain network model, the symmetric matrix  𝐀∈ℝ^N × N is the state transition matrix describing the interactions between brain nodes; the diagonal matrix 𝐁∈ℝ^N × N is the input matrix with the value in diagonal either 1, indicating the brain region is selected as a control node, or 0, indicating the region is not a control node; the vector 𝐮(t) ∈ℝ^N represents the control strategy to the brain network at time t. 


In this study, our goal is to infer the input signals of the brain system during specific tasks according to the observed fMRI data. We define 𝐀 as the structural connectivity of the brain regions, which can be obtained by DTI imaging. The system dynamics 𝐱(t) are observed, reflected in fMRI data. The parameters 𝐁 and 𝐮(t) remain unknown.
Technically, the key challenge lies in identifying the binary input matrix 𝐁—representing critical control nodes in the brain network. To this end, we propose a neural perturbation-based approach to identify 𝐁, and employ a data-driven approach to infer the input signal 𝐮(t) from the task fMRI data. 







 §.§ Objective and optimization of input identification




  §.§.§ Objective function and contraints

Since we model the brain system with explicit expression defined by Eq. (<ref>), the input identification problem can be conceptualized as a parameter estimation problem by fitting data.
Inspired by machine learning and neuroscience prior knowledge, we define the objective function with three terms, including a data reconstruction loss to minimize the reconstruction error and two input regularization terms to ensure the smoothness of brain inputs. The brain system is constrained by the linear network dynamic model and input matrix 𝐁. The objective function with system constraints is shown in Eq. (<ref>),



    min _B, U_1:T   ∑_t=1^T𝐱(t)-𝐱(t)_2^2+λ_1 ∑_t=1^T𝐮(t)_2^2 
       +λ_2 ∑_t=1^T-1𝐮(t+1)-𝐮(t)_2^2
     s.t. :     𝐱(t+1)=𝐀𝐱(t)+𝐁𝐮(t) 
    
       ∑_i=1^N 𝐁_i,i =α,     𝐁_i,i∈{0,1},

where 𝐱(t) ∈ℝ^N is the state of the brain network at time step t; 𝐱(t) is the generative state in the brain network control framework; 𝐀 is an N × N adjacency matrix that describes the structural connections among brain regions; 𝐁 is the control matrix where 𝐁_i,i∈{0, 1} for i={1,2,…, N}, indicating whether the node i is a critical node in the specific task. There are in total three hyperparameters, i.e., λ_1, λ_2, and α. We use the hyperparameters λ_1 and λ_2 to constrain the energy and the slope of input signals, respectively. The hyperparameter α determines the number of identified critical nodes.




  §.§.§ Optimization Algorithm

Note that Eq (<ref>) is a nonconvex optimization problem due to the Boolean constraints 𝐁_i,i∈{0, 1}, making it a challenge to solve directly. Here, we relax the Boolean constraint using continuous epsilon relaxation with |𝐁_i,i (1-𝐁_i,i)| ≤ϵ,
which yields the following optimization problem,

    min _𝐁, 𝐮_1:T   ∑_t=1^T𝐱(t)-𝐱(t)_2^2+λ_1 ∑_t=1^T𝐮(t)_2^2 
       +λ_2 ∑_t=1^T-1𝐮(t+1)-𝐮(t)_2^2
     s.t. :     𝐱(t+1)=𝐀𝐱(t)+𝐁 u(t) 
       ∑_i=1^N 𝐁_i,i =α,    |𝐁_i,i (1-𝐁_i,i)| ≤ϵ_k,

where ϵ_k is a violation parameter that smoothes the integer constraints and it will converge to a sufficiently small number during the optimization. 

To solve the optimization problem in Eq. (<ref>), 

we utilize the Augmented Lagrangian Method (ALM) with gradient descent for iterative optimization, progressively converging towards the optimal solution. 
Importantly, similar to training recurrent neural networks, the gradient of ℒ to B and U can be computed automatically by the AUTOGRAD in Pytorch and updated efficiently by backpropagation through time (BPTT).











            



  

























§ RESULTS

We apply the framework to both synthetic data generated by a pre-defined linear dynamic model and the motor task-fMRI data to validate the efficacy of our approach.



 §.§ Validation with simulated data

We first validate the performance of the proposed algorithm by addressing the reverse problem of deducing control nodes and their corresponding inputs from synthetic data generated by a pre-defined linear model. Specifically, we randomly initialize a state-transition matrix 𝐀 with 6 nodes, and two input sequences sin(10π t) and cos(10π t) to the node 2 & 3, respectively. Following the linear model in Eq. (<ref>), we synthesize the evolving dynamics (Fig. <ref>a). 
Given the synthetic data, we apply our method to iteratively update 𝐁 and 𝐮, and ultimately it converges to a steady solution (Fig. <ref>b). The quantitative results demonstrate the effectiveness of our method on both clean and noisy data (Fig. <ref>c).






 §.§ Identifying the inputs using fMRI data under motor task




  §.§.§ Motor task and fMRI recording

To validate the proposed framework, we utilized the task-fMRI dataset from the Human Connectome Project (HCP 1200 release). We randomly selected 200 out of 1200 subjects. All selected subjects include T1 MRI, dMRI, and fMRI recording during the motor tasks. The motor tasks comprised of 5 sub-tasks, including left/right foot, left/right hand, and tongue. Each motor task includes 17 Repetition Times (TRs) with a TR=0.72 second.




We parcellate the brain into 100 non-overlapping regions of interest (ROIs) using Yeo atlas<cit.>. 
The 100×100 structural connectivity is obtained from DTI images, and 100 × 17 × trials time courses are extracted from motor task-fMRI data (trials = 200 subjects × 5 sub-tasks). 



  §.§.§ Identified control nodes and inputs

We extend our framework to identify the critical control nodes with empirical task-fMRI data. For the group-level analysis, we first apply two metrics to quantify the performance of fitting accuracy across max control nodes and regions of the whole brain, including the Explained-Variance (EV), and the  Bayesian Information Criteria (BIC). EV quantifies the fitting accuracy of the task-related functional dynamics of the whole brain across the number of control nodes (From 25 to 45) (Fig. <ref>a). Fig. <ref>a (top) shows the mean and variance of global whole-brain EV across tasks with different numbers of control nodes in each subject. The curves increase with the increase of α until they reach a value with high predictive performance. Meanwhile, we apply BIC to select the optimal number of control nodes. Fig. <ref>a (bottom) shows the mean and variance of BIC across tasks with different numbers of control nodes in each subject. It shows that a bell-shaped curve with the increase of α, and the optimal number of control nodes is α=28 of motor tasks. 

To interpret the identified nodes, we compare them with the average-level activation map of motor tasks in the HCP dataset (Fig. <ref>b&c). We found that most of the inferred control nodes are well in line with the activated map, largely overlapping in motor and somatosensory areas. The inferred input sequences of ROIs are illustrated in Fig. <ref>d.

















§ DISCUSSION AND CONCLUSION




  §.§.§ From methodology perspective

Our study presented a general framework based on network control theory for identifying brain inputs during cognitive tasks. We proposed the Augmented Lagrangian Method with gradient descent to iteratively search the steady state of the objective function, which guarantees finding the optimal with the sparse control nodes and the smooth input sequences. The framework was validated with both synthetic data and real task-based fMRI data from the HCP dataset. 
Different from statistical analyses <cit.> and machine learning methods <cit.>, we brought the network control theory to node identification. Rather than focusing on partially observed data, we use the whole-brain fMRI to infer the system parameters 𝐁 and 𝐮(t).




  §.§.§ From neuroscience perspective





Revealing how external stimuli are processed by the brain and dynamically transmitted in different brain regions is the core question in neuroscience. Our work delves into the brain's input-output response to task stimuli, focusing on two key aspects: identifying control nodes and inferring input sequences.
Identifying control nodes enables us to pinpoint core brain regions relevant to cognitive functions. Also, it is important in clinical applications, such as selecting electrical stimulation sites for precise brain function modulation<cit.>.
Moreover, characterizing the input sequences is key to understanding how information flows between brain regions <cit.>. The temporal dynamics elucidate the information process across brain regions and the integration and coordination of information during complex cognitive tasks.






In summary, we proposed a new framework from control theory perspective for understanding brain inputs and system dynamics during cognitive tasks. It is a potential tool for unraveling the input-output relations in the brain under task stimulation and electrical stimulation.


















































































































































































































IEEEtran

