By name: brain-network-controllability description: 结构脑网络可控性分析实用指南。基于网络控制理论计算平均可控性、模态可控性、最小控制能量等指标，研究脑网络状态转移能力。适用于脑网络分析、神经调控研究、认知控制研究。触发词：脑网络可控性、控制理论、结构连接、控制能量、平均可控性、模态可控性、controllability、control theory、structural brain network。 user-invocable: true
结构脑网络可控性分析方法论
来源论文： arXiv:1908.03514 - A practical guide to methodological considerations in the controllability of structural brain networks
核心方法论
基于网络控制理论的脑网络可控性框架：
1. 线性动力学模型
[\mathbf{x}(t+1) = \mathbf{A}\mathbf{x}(t) + \mathbf{B}\mathbf{u}(t)]
其中：
- (\mathbf{x}(t))：脑区状态向量
- (\mathbf{A})：结构连接矩阵（状态转移矩阵）
- (\mathbf{B})：控制输入矩阵
- (\mathbf{u}(t))：控制信号
2. 可控性指标
平均可控性（Average Controllability）
- 驱动系统到附近状态的能力
- 与网络整合和扩散相关
模态可控性（Modal Controllability）
- 驱动系统到远距离状态的能力
- 与网络分离和专门化相关
控制能量（Control Energy）
- 实现状态转移所需的最小能量
- 与认知努力相关
Python 实现
import numpy as np
from scipy.linalg import eig, inv, solve, svd
from scipy.integrate import solve_ivp
from typing import Dict, Tuple, List, Optional
from dataclasses import dataclass
import warnings
@dataclass
class ControllabilityConfig:
"""可控性分析配置"""
normalize: bool = True # 是否归一化可控性指标
energy_horizon: float = 1.0 # 控制能量计算时间范围
timesteps: int = 100 # 状态转移时间步数
regularization: float = 1e-6 # 正则化参数
eigenvalue_threshold: float = 0.99 # 特征值阈值（稳定性）
class BrainNetworkControllability:
"""脑网络可控性分析器"""
def __init__(self, A: np.ndarray, config: Optional[ControllabilityConfig] = None):
"""
Args:
A: 结构连接矩阵 (n x n)
config: 分析配置
"""
self.A = A.copy().astype(float)
self.n = A.shape[0]
self.config = config or ControllabilityConfig()
# 稳定性处理
self._stabilize_matrix()
def _stabilize_matrix(self):
"""确保系统稳定（谱半径 < 1）"""
eigenvalues = np.linalg.eigvals(self.A)
spectral_radius = np.max(np.abs(eigenvalues))
if spectral_radius >= 1:
# 缩放矩阵
self.A = self.A / spectral_radius * self.config.eigenvalue_threshold
warnings.warn(f"Matrix scaled for stability. Original spectral radius: {spectral_radius:.4f}")
def gramian(self, T: int = None) -> np.ndarray:
"""计算可控性格拉姆矩阵
W_c = sum_{t=0}^{T-1} A^t @ B @ B^T @ (A^t)^T
Args:
T: 时间范围
Returns:
W_c: 可控性格拉姆矩阵 (n x n)
"""
if T is None:
T = self.config.timesteps
# 假设B = I（所有节点都是控制节点）
B = np.eye(self.n)
W_c = np.zeros((self.n, self.n))
A_power = np.eye(self.n)
for t in range(T):
W_c += A_power @ B @ B.T @ A_power.T
A_power = A_power @ self.A
return W_c
def average_controllability(self, node_indices: List[int] = None) -> np.ndarray:
"""计算平均可控性
对于控制节点i，平均可控性定义为：
AC_i = trace(W_c^{(i)})
其中 W_c^{(i)} 是以节点i为控制输入的格拉姆矩阵
Args:
node_indices: 要分析的节点索引，None则分析所有节点
Returns:
ac: 每个节点的平均可控性 (n,)
"""
if node_indices is None:
node_indices = list(range(self.n))
ac = np.zeros(len(node_indices))
T = self.config.timesteps
for idx, i in enumerate(node_indices):
# B矩阵：只有第i列非零
B = np.zeros((self.n, 1))
B[i, 0] = 1
# 计算格拉姆矩阵
W_c = np.zeros((self.n, self.n))
A_power = np.eye(self.n)
for t in range(T):
W_c += A_power @ B @ B.T @ A_power.T
A_power = A_power @ self.A
# 平均可控性 = trace(W_c)
ac[idx] = np.trace(W_c)
if self.config.normalize:
ac = ac / ac.max() if ac.max() > 0 else ac
return ac
def modal_controllability(self, node_indices: List[int] = None) -> np.ndarray:
"""计算模态可控性
基于特征向量计算：
phi_i = 1 - sum_j (A_ij)^2 * lambda_j^2 / ||v_j||^2
其中 lambda_j 是特征值，v_j 是特征向量
Args:
node_indices: 要分析的节点索引
Returns:
mc: 每个节点的模态可控性 (n,)
"""
if node_indices is None:
node_indices = list(range(self.n))
# 特征分解
eigenvalues, eigenvectors = eig(self.A)
eigenvalues = np.real(eigenvalues)
eigenvectors = np.real(eigenvectors)
# 计算每个节点的模态可控性
mc = np.zeros(len(node_indices))
for idx, i in enumerate(node_indices):
# 节点i的模态可控性
# phi_i = sum_j (1 - lambda_j^2) * v_{i,j}^2
mc[idx] = np.sum(
(1 - eigenvalues**2) * eigenvectors[i, :]**2
)
if self.config.normalize:
mc = mc / mc.max() if mc.max() > 0 else mc
return mc
def minimum_control_energy(self, x0: np.ndarray, xf: np.ndarray,
T: float = 1.0,
controlled_nodes: List[int] = None) -> float:
"""计算最小控制能量
E = integral_0^T u(t)^T u(t) dt
使用最优控制理论求解
Args:
x0: 初始状态 (n,)
xf: 目标状态 (n,)
T: 控制时间
controlled_nodes: 控制节点索引
Returns:
E: 最小控制能量
"""
if controlled_nodes is None:
controlled_nodes = list(range(self.n))
# 控制矩阵B
B = np.zeros((self.n, len(controlled_nodes)))
for j, i in enumerate(controlled_nodes):
B[i, j] = 1
# 计算控制能量
# 离散时间控制
timesteps = int(T * 100)
dt = T / timesteps
# 状态转移
def dynamics(t, x, u_func):
u = u_func(t)
return self.A @ x + B @ u
# 最优控制：u*(t) = B^T @ exp(A^T(T-t)) @ W^{-1} @ (xf - exp(AT) @ x0)
# 其中 W = integral_0^T exp(At) @ B @ B^T @ exp(A^T t) dt
# 简化计算：使用伪逆
delta_x = xf - np.linalg.matrix_power(self.A, timesteps) @ x0
# 计算可控性格拉姆矩阵
W = np.zeros((self.n, self.n))
A_power = np.eye(self.n)
for t in range(timesteps):
W += A_power @ B @ B.T @ A_power.T * dt
A_power = A_power @ self.A
# 正则化
W += self.config.regularization * np.eye(self.n)
# 最优控制输入
try:
v = solve(W, delta_x, assume_a='pos')
except np.linalg.LinAlgError:
v = np.linalg.lstsq(W, delta_x, rcond=None)[0]
# 计算能量
E = delta_x.T @ v
return max(0, E)
def optimal_control_trajectory(self, x0: np.ndarray, xf: np.ndarray,
T: float = 1.0,
controlled_nodes: List[int] = None) -> Dict:
"""计算最优控制轨迹
Args:
x0: 初始状态
xf: 目标状态
T: 控制时间
controlled_nodes: 控制节点
Returns:
trajectory: 包含状态和控制轨迹的字典
"""
if controlled_nodes is None:
controlled_nodes = list(range(self.n))
n_controls = len(controlled_nodes)
timesteps = int(T * 100)
dt = T / timesteps
# B矩阵
B = np.zeros((self.n, n_controls))
for j, i in enumerate(controlled_nodes):
B[i, j] = 1
# 计算格拉姆矩阵
W = np.zeros((self.n, self.n))
A_power = np.eye(self.n)
for t in range(timesteps):
W += A_power @ B @ B.T @ A_power.T * dt
A_power = A_power @ self.A
W += self.config.regularization * np.eye(self.n)
# 状态轨迹
t = np.linspace(0, T, timesteps + 1)
x = np.zeros((timesteps + 1, self.n))
x[0] = x0
# 控制轨迹
u = np.zeros((timesteps, n_controls))
# 最优控制计算
delta_x = xf - np.linalg.matrix_power(self.A, timesteps) @ x0
try:
v = solve(W, delta_x, assume_a='pos')
except np.linalg.LinAlgError:
v = np.linalg.lstsq(W, delta_x, rcond=None)[0]
# 模拟轨迹
A_power = np.eye(self.n)
for i in range(timesteps):
# 控制输入
exp_At = np.linalg.matrix_power(self.A, timesteps - i)
u[i] = B.T @ exp_At.T @ v
# 更新状态
x[i+1] = self.A @ x[i] + B @ u[i]
return {
'time': t,
'states': x,
'controls': u,
'energy': np.sum(u**2) * dt
}
def radial_propagation_connectivity(self,
distances: np.ndarray) -> np.ndarray:
"""计算径向传播连接性
考虑活动通过相邻组织的径向传播
Args:
distances: 脑区之间的物理距离矩阵 (n x n)
Returns:
A_radial: 径向传播连接矩阵
"""
# 径向传播因子
# 连接强度随距离衰减
sigma = np.mean(distances[distances > 0]) / 2 # 特征长度
A_radial = np.exp(-distances**2 / (2 * sigma**2))
# 结合结构连接
A_combined = self.A * A_radial
return A_combined
def energy_landscape_complexity(self,
n_samples: int = 100) -> Dict:
"""计算能量景观复杂度
随机采样状态转移，分析能量分布
Args:
n_samples: 采样数量
Returns:
metrics: 能量景观度量
"""
energies = []
for _ in range(n_samples):
# 随机初始和目标状态
x0 = np.random.randn(self.n)
xf = np.random.randn(self.n)
# 计算最小控制能量
E = self.minimum_control_energy(x0, xf)
energies.append(E)
energies = np.array(energies)
return {
'mean_energy': np.mean(energies),
'std_energy': np.std(energies),
'min_energy': np.min(energies),
'max_energy': np.max(energies),
'energy_range': np.max(energies) - np.min(energies),
'energy_skewness': self._skewness(energies),
'energy_kurtosis': self._kurtosis(energies)
}
@staticmethod
def _skewness(x: np.ndarray) -> float:
"""计算偏度"""
return np.mean((x - np.mean(x))**3) / np.std(x)**3
@staticmethod
def _kurtosis(x: np.ndarray) -> float:
"""计算峰度"""
return np.mean((x - np.mean(x))**4) / np.std(x)**4 - 3
def full_analysis(self,
x0: np.ndarray = None,
xf: np.ndarray = None) -> Dict:
"""完整可控性分析
Returns:
results: 所有可控性指标
"""
results = {}
# 1. 平均可控性
results['average_controllability'] = self.average_controllability()
# 2. 模态可控性
results['modal_controllability'] = self.modal_controllability()
# 3. 可控性权衡
results['controllability_tradeoff'] = -np.corrcoef(
results['average_controllability'],
results['modal_controllability']
)[0, 1]
# 4. 控制能量（如果提供了状态）
if x0 is not None and xf is not None:
results['minimum_energy'] = self.minimum_control_energy(x0, xf)
results['optimal_trajectory'] = self.optimal_control_trajectory(x0, xf)
# 5. 能量景观复杂度
results['energy_landscape'] = self.energy_landscape_complexity(n_samples=50)
# 6. 网络统计
results['network_stats'] = {
'spectral_radius': np.max(np.abs(np.linalg.eigvals(self.A))),
'density': np.mean(self.A > 0),
'mean_weight': np.mean(self.A[self.A > 0]) if np.any(self.A > 0) else 0
}
return results
def generate_modular_network(n_nodes: int, n_modules: int = 4,
p_intra: float = 0.3, p_inter: float = 0.05,
weight_mean: float = 0.5, weight_std: float = 0.2) -> np.ndarray:
"""生成模块化脑网络
Args:
n_nodes: 节点数量
n_modules: 模块数量
p_intra: 模块内连接概率
p_inter: 模块间连接概率
weight_mean: 连接权重均值
weight_std: 连接权重标准差
Returns:
A: 连接矩阵
"""
A = np.zeros((n_nodes, n_nodes))
module_size = n_nodes // n_modules
rng = np.random.default_rng(42)
for i in range(n_nodes):
for j in range(n_nodes):
if i == j:
continue
# 确定是否在同一模块
i_module = i // module_size
j_module = j // module_size
if i_module == j_module:
p = p_intra
else:
p = p_inter
# 建立连接
if rng.random() < p:
A[i, j] = np.clip(
rng.normal(weight_mean, weight_std),
0.1, 1.0
)
# 对称化
A = (A + A.T) / 2
return A
# 使用示例
def example_analysis():
"""示例：脑网络可控性分析"""
# 生成模拟连接矩阵
n_nodes = 50
A = generate_modular_network(n_nodes, n_modules=5)
# 创建分析器
analyzer = BrainNetworkControllability(A)
# 完整分析
results = analyzer.full_analysis()
print("="*60)
print("脑网络可控性分析结果")
print("="*60)
print(f"\n网络统计:")
print(f" 谱半径: {results['network_stats']['spectral_radius']:.4f}")
print(f" 连接密度: {results['network_stats']['density']:.4f}")
print(f" 平均权重: {results['network_stats']['mean_weight']:.4f}")
print(f"\n可控性指标:")
ac = results['average_controllability']
mc = results['modal_controllability']
print(f" 平均可控性: mean={ac.mean():.4f}, std={ac.std():.4f}")
print(f" 模态可控性: mean={mc.mean():.4f}, std={mc.std():.4f}")
print(f" 可控性权衡 (负相关): r={results['controllability_tradeoff']:.4f}")
print(f"\n能量景观:")
el = results['energy_landscape']
print(f" 平均能量: {el['mean_energy']:.4f}")
print(f" 能量标准差: {el['std_energy']:.4f}")
print(f" 能量范围: {el['energy_range']:.4f}")
# 识别高可控性节点
print(f"\n高平均可控性节点 (前5):")
top_ac = np.argsort(ac)[-5:][::-1]
for i, node in enumerate(top_ac):
print(f" {i+1}. 节点 {node}: AC={ac[node]:.4f}")
print(f"\n高模态可控性节点 (前5):")
top_mc = np.argsort(mc)[-5:][::-1]
for i, node in enumerate(top_mc):
print(f" {i+1}. 节点 {node}: MC={mc[node]:.4f}")
return analyzer, results
def example_state_transition():
"""示例：状态转移控制"""
n_nodes = 20
A = generate_modular_network(n_nodes)
analyzer = BrainNetworkControllability(A)
# 定义初始和目标状态
x0 = np.zeros(n_nodes)
x0[0:5] = 1.0 # 前5个节点激活
xf = np.zeros(n_nodes)
xf[15:20] = 1.0 # 后5个节点激活
# 计算最优轨迹
trajectory = analyzer.optimal_control_trajectory(x0, xf, T=1.0)
print("\n状态转移控制:")
print(f" 控制能量: {trajectory['energy']:.4f}")
print(f" 时间步数: {len(trajectory['time'])}")
print(f" 初始状态范数: {np.linalg.norm(x0):.4f}")
print(f" 最终状态范数: {np.linalg.norm(trajectory['states'][-1]):.4f}")
return trajectory
## Activation Keywords
- 脑网络可控性
- 控制理论
- 结构连接
- 控制能量
- 平均可控性
- 模态可控性
- controllability
- control theory
- structural brain network
- 能量景观
## Tools Used
- numpy
- scipy
- networkx
- matplotlib
## Instructions for Agents
1. 理解线性动力学模型：x(t+1) = Ax(t) + Bu(t)
2. 计算平均可控性：驱动系统到附近状态的能力
3. 计算模态可控性：驱动系统到远距离状态的能力
4. 计算控制能量：实现状态转移所需的最小能量
5. 分析可控性权衡：平均可控性与模态可控性的负相关关系
## Examples
```python
# 使用示例
from brain_network_controllability import BrainNetworkControllability, generate_modular_network
# 1. 生成模块化脑网络
A = generate_modular_network(n_nodes=50, n_modules=5)
# 2. 创建可控性分析器
analyzer = BrainNetworkControllability(A)
# 3. 完整分析
results = analyzer.full_analysis()
print(f"平均可控性: mean={results['average_controllability'].mean():.4f}")
print(f"模态可控性: mean={results['modal_controllability'].mean():.4f}")
# 4. 状态转移控制
x0 = np.zeros(50)
xf = np.zeros(50)
xf[15:20] = 1.0
energy = analyzer.minimum_control_energy(x0, xf)
if name == "main": analyzer, results = example_analysis() print("\n" + "="*60) trajectory = example_state_transition()