Linear Structure Function Coupling

---
name: linear-structure-function-coupling
description: 脑结构-功能耦合线性生成框架。从扩散加权成像(DWI)结构连接预测静息态fMRI功能连接，揭示整合枢纽和中介枢纽的机制，支持虚拟损伤实验。适用于脑连接组学、神经影像分析、结构-功能关系研究。触发词：结构功能耦合、脑连接组、结构连接、功能连接、整合枢纽、structure-function coupling、structural connectivity、functional connectivity、integrator hub、mediator hub。
user-invocable: true
---

# 脑结构-功能耦合线性生成框架

**来源论文：** arXiv:2507.06136 - A Linear Generative Framework for Structure-Function Coupling in the Human Brain

## 核心方法论

### 1. 结构-功能耦合问题

**结构连接 (SC)：** 白质通路网络（DWI 测量）

**功能连接 (FC)：** 脑区神经活动相关性（rs-fMRI 测量）

**核心问题：** 结构架构如何塑造功能模式？

### 2. 线性生成模型

**核心方程：**

\[
FC = \sum_{k} \alpha_k \cdot M_k(SC)
\]

其中：
- \(M_k\) - 结构动机（motif）变换
- \(\alpha_k\) - 权重系数
- 通过直接/间接通路预测 FC

### 3. 枢纽分类

| 类型 | 功能 | 作用 |
|------|------|------|
| **整合枢纽** | 结构支点 | 促进同步 |
| **中介枢纽** | 结构杠杆 | 协调竞争动力学 |

### 4. 虚拟损伤实验

预测连接破坏如何级联影响功能网络

## Python 实现

```python
import numpy as np
from typing import Dict, List, Tuple, Optional
from dataclasses import dataclass
from collections import defaultdict
from scipy import stats
from sklearn.linear_model import Ridge
from sklearn.model_selection import cross_val_score


@dataclass
class StructureFunctionConfig:
    """结构-功能耦合配置"""
    n_regions: int = 100           # 脑区数量
    n_motifs: int = 5              # 结构动机数量
    
    # 模型参数
    regularization: float = 1.0    # 正则化参数
    
    # 枢纽识别
    hub_threshold: float = 0.8     # 枢纽阈值（百分位）


class StructuralMotifExtractor:
    """结构动机提取器"""
    
    def __init__(self, config: StructureFunctionConfig):
        self.config = config
        
    def extract_direct_connections(self, SC: np.ndarray) -> np.ndarray:
        """提取直接连接动机
        
        Args:
            SC: 结构连接矩阵
            
        Returns:
            direct: 直接连接动机
        """
        return SC.copy()
    
    def extract_indirect_connections(self, 
                                      SC: np.ndarray,
                                      max_path_length: int = 3) -> List[np.ndarray]:
        """提取间接连接动机
        
        Args:
            SC: 结构连接矩阵
            max_path_length: 最大路径长度
            
        Returns:
            indirect_motifs: 间接连接动机列表
        """
        motifs = []
        
        # 2步路径
        path_2 = SC @ SC
        motifs.append(path_2 / (path_2.max() + 1e-10))
        
        # 3步路径
        if max_path_length >= 3:
            path_3 = SC @ path_2
            motifs.append(path_3 / (path_3.max() + 1e-10))
            
        return motifs
    
    def extract_triadic_motifs(self, SC: np.ndarray) -> np.ndarray:
        """提取三元动机
        
        Args:
            SC: 结构连接矩阵
            
        Returns:
            triadic: 三元动机
        """
        n = SC.shape[0]
        triadic = np.zeros((n, n))
        
        for i in range(n):
            for j in range(n):
                # 计算 i→j 的三元动机强度
                for k in range(n):
                    if k != i and k != j:
                        triadic[i, j] += SC[i, k] * SC[k, j]
                        
        return triadic / (triadic.max() + 1e-10)
    
    def extract_all_motifs(self, SC: np.ndarray) -> np.ndarray:
        """提取所有结构动机
        
        Args:
            SC: 结构连接矩阵
            
        Returns:
            motifs: 动机矩阵 (n_motifs, n_regions, n_regions)
        """
        motifs = []
        
        # 直接连接
        motifs.append(self.extract_direct_connections(SC))
        
        # 间接连接
        indirect = self.extract_indirect_connections(SC)
        motifs.extend(indirect)
        
        # 三元动机
        motifs.append(self.extract_triadic_motifs(SC))
        
        return np.array(motifs)


class LinearGenerativeModel:
    """线性生成模型"""
    
    def __init__(self, config: StructureFunctionConfig):
        self.config = config
        self.motif_extractor = StructuralMotifExtractor(config)
        
        # 模型权重
        self.weights = None
        self.motifs = None
        
    def fit(self, SC: np.ndarray, FC: np.ndarray) -> Dict:
        """拟合模型
        
        Args:
            SC: 结构连接矩阵 (n_regions, n_regions)
            FC: 功能连接矩阵 (n_regions, n_regions)
            
        Returns:
            training_info: 训练信息
        """
        # 提取动机
        self.motifs = self.motif_extractor.extract_all_motifs(SC)
        n_motifs = len(self.motifs)
        
        # 准备训练数据
        n_regions = SC.shape[0]
        
        # 展平动机
        X = self.motifs.reshape(n_motifs, -1).T  # (n_regions^2, n_motifs)
        y = FC.flatten()  # (n_regions^2,)
        
        # 岭回归
        model = Ridge(alpha=self.config.regularization)
        model.fit(X, y)
        
        self.weights = model.coef_
        
        # 交叉验证
        cv_scores = cross_val_score(model, X, y, cv=5)
        
        return {
            'weights': self.weights,
            'cv_score_mean': cv_scores.mean(),
            'cv_score_std': cv_scores.std(),
            'n_motifs': n_motifs
        }
    
    def predict(self, SC: np.ndarray) -> np.ndarray:
        """预测功能连接
        
        Args:
            SC: 结构连接矩阵
            
        Returns:
            FC_pred: 预测的功能连接
        """
        if self.weights is None:
            raise ValueError("Model not fitted. Call fit() first.")
            
        # 提取动机
        motifs = self.motif_extractor.extract_all_motifs(SC)
        
        # 加权组合
        FC_pred = np.zeros((SC.shape[0], SC.shape[0]))
        
        for i, w in enumerate(self.weights):
            if i < len(motifs):
                FC_pred += w * motifs[i]
                
        # 对称化
        FC_pred = (FC_pred + FC_pred.T) / 2
        
        return FC_pred
    
    def compute_prediction_accuracy(self,
                                     FC_true: np.ndarray,
                                     FC_pred: np.ndarray) -> Dict:
        """计算预测准确度
        
        Args:
            FC_true: 真实 FC
            FC_pred: 预测 FC
            
        Returns:
            accuracy: 准确度指标
        """
        # 展平
        true = FC_true.flatten()
        pred = FC_pred.flatten()
        
        # 移除对角线
        n = FC_true.shape[0]
        mask = ~np.eye(n, dtype=bool)
        
        true_masked = true[mask.flatten()]
        pred_masked = pred[mask.flatten()]
        
        # 相关系数
        correlation = np.corrcoef(true_masked, pred_masked)[0, 1]
        
        # R²
        r2 = 1 - np.sum((true_masked - pred_masked)**2) / np.sum((true_masked - true_masked.mean())**2)
        
        return {
            'correlation': correlation,
            'r2': r2,
            'rmse': np.sqrt(np.mean((true_masked - pred_masked)**2))
        }


class HubClassifier:
    """枢纽分类器"""
    
    def __init__(self, config: StructureFunctionConfig):
        self.config = config
        
    def classify_hubs(self, 
                      SC: np.ndarray,
                      FC: np.ndarray) -> Dict:
        """分类枢纽
        
        Args:
            SC: 结构连接矩阵
            FC: 功能连接矩阵
            
        Returns:
            classification: 分类结果
        """
        n = SC.shape[0]
        
        # 计算节点强度
        sc_strength = SC.sum(axis=1) + SC.sum(axis=0)
        fc_strength = FC.sum(axis=1) + FC.sum(axis=0)
        
        # 计算参与系数
        participation = self._compute_participation_coefficient(SC)
        
        # 识别枢纽
        sc_threshold = np.percentile(sc_strength, self.config.hub_threshold * 100)
        hub_mask = sc_strength >= sc_threshold
        
        # 分类
        integrator_hubs = []
        mediator_hubs = []
        
        for i in range(n):
            if hub_mask[i]:
                # 高参与系数 → 整合枢纽
                # 低参与系数但高强度 → 中介枢纽
                if participation[i] > np.median(participation[hub_mask]):
                    integrator_hubs.append(i)
                else:
                    mediator_hubs.append(i)
                    
        return {
            'integrator_hubs': integrator_hubs,
            'mediator_hubs': mediator_hubs,
            'sc_strength': sc_strength,
            'fc_strength': fc_strength,
            'participation': participation
        }
    
    def _compute_participation_coefficient(self, SC: np.ndarray) -> np.ndarray:
        """计算参与系数
        
        Args:
            SC: 结构连接矩阵
            
        Returns:
            participation: 参与系数
        """
        n = SC.shape[0]
        
        # 简化的社区检测（使用模块度优化）
        # 这里用度排序作为代理
        degree = SC.sum(axis=1)
        
        # 分成两个社区
        sorted_indices = np.argsort(degree)
        community = np.zeros(n, dtype=int)
        community[sorted_indices[n//2:]] = 1
        
        # 计算参与系数
        participation = np.zeros(n)
        
        for i in range(n):
            s = SC[i, :].sum()
            if s == 0:
                continue
                
            # 每个社区的连接比例
            for c in [0, 1]:
                s_c = SC[i, community == c].sum()
                p_c = s_c / s
                participation[i] += p_c**2
                
            participation[i] = 1 - participation[i]
            
        return participation


class VirtualLesionExperiment:
    """虚拟损伤实验"""
    
    def __init__(self, model: LinearGenerativeModel):
        self.model = model
        
    def simulate_lesion(self,
                        SC: np.ndarray,
                        lesion_region: int) -> np.ndarray:
        """模拟损伤
        
        Args:
            SC: 结构连接矩阵
            lesion_region: 损伤脑区
            
        Returns:
            SC_lesioned: 损伤后的结构连接
        """
        SC_lesioned = SC.copy()
        
        # 移除该区域的所有连接
        SC_lesioned[lesion_region, :] = 0
        SC_lesioned[:, lesion_region] = 0
        
        return SC_lesioned
    
    def analyze_cascade_effect(self,
                                SC: np.ndarray,
                                FC_original: np.ndarray,
                                lesion_regions: List[int]) -> Dict:
        """分析级联效应
        
        Args:
            SC: 结构连接矩阵
            FC_original: 原始功能连接
            lesion_regions: 损伤脑区列表
            
        Returns:
            cascade: 级联效应分析
        """
        results = []
        
        for region in lesion_regions:
            # 模拟损伤
            SC_lesioned = self.simulate_lesion(SC, region)
            
            # 预测新的 FC
            FC_pred = self.model.predict(SC_lesioned)
            
            # 计算变化
            fc_change = FC_pred - FC_original
            global_change = np.mean(np.abs(fc_change))
            
            # 局部变化
            local_change = np.mean(np.abs(fc_change[region, :]))
            
            results.append({
                'lesion_region': region,
                'global_fc_change': global_change,
                'local_fc_change': local_change,
                'fc_predicted': FC_pred
            })
            
        return {
            'lesion_results': results,
            'most_disruptive_region': max(results, 
                                         key=lambda x: x['global_fc_change'])['lesion_region']
        }


def visualize_structure_function_coupling(SC: np.ndarray,
                                           FC: np.ndarray,
                                           FC_pred: np.ndarray):
    """可视化结构-功能耦合"""
    import matplotlib.pyplot as plt
    
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 1. 结构连接
    ax = axes[0, 0]
    im = ax.imshow(SC, cmap='hot')
    ax.set_title('Structural Connectivity (SC)')
    ax.set_xlabel('Brain Region')
    ax.set_ylabel('Brain Region')
    plt.colorbar(im, ax=ax)
    
    # 2. 真实功能连接
    ax = axes[0, 1]
    im = ax.imshow(FC, cmap='coolwarm', vmin=-1, vmax=1)
    ax.set_title('Functional Connectivity (FC) - True')
    ax.set_xlabel('Brain Region')
    ax.set_ylabel('Brain Region')
    plt.colorbar(im, ax=ax)
    
    # 3. 预测功能连接
    ax = axes[1, 0]
    im = ax.imshow(FC_pred, cmap='coolwarm', vmin=-1, vmax=1)
    ax.set_title('FC - Predicted')
    ax.set_xlabel('Brain Region')
    ax.set_ylabel('Brain Region')
    plt.colorbar(im, ax=ax)
    
    # 4. 预测 vs 真实
    ax = axes[1, 1]
    n = SC.shape[0]
    mask = ~np.eye(n, dtype=bool)
    ax.scatter(FC[mask], FC_pred[mask], alpha=0.3, s=10)
    ax.plot([-1, 1], [-1, 1], 'r--', label='Identity')
    ax.set_xlabel('True FC')
    ax.set_ylabel('Predicted FC')
    ax.set_title(f'Prediction Accuracy: r={np.corrcoef(FC[mask], FC_pred[mask])[0,1]:.3f}')
    ax.legend()
    ax.set_xlim(-1, 1)
    ax.set_ylim(-1, 1)
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('structure_function_coupling.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    return 'structure_function_coupling.png'


# 使用示例
def example_structure_function_coupling():
    """示例：结构-功能耦合分析"""
    print("="*60)
    print("脑结构-功能耦合线性生成框架")
    print("="*60)
    
    config = StructureFunctionConfig(n_regions=100)
    
    # 生成模拟数据
    n = config.n_regions
    
    # 结构连接（稀疏）
    SC = np.random.rand(n, n) * 0.3
    SC = (SC + SC.T) / 2
    SC[np.random.rand(n, n) > 0.1] = 0  # 稀疏化
    np.fill_diagonal(SC, 0)
    
    # 功能连接（基于 SC 加噪声）
    FC = SC @ SC + np.random.randn(n, n) * 0.1
    FC = (FC + FC.T) / 2
    np.fill_diagonal(FC, 0)
    FC = FC / (np.abs(FC).max() + 1e-10)
    
    # 创建模型
    model = LinearGenerativeModel(config)
    
    # 拟合
    print("\n拟合模型...")
    training_info = model.fit(SC, FC)
    print(f"  动机数量: {training_info['n_motifs']}")
    print(f"  CV 分数: {training_info['cv_score_mean']:.3f} ± {training_info['cv_score_std']:.3f}")
    
    # 预测
    print("\n预测功能连接...")
    FC_pred = model.predict(SC)
    
    # 评估
    accuracy = model.compute_prediction_accuracy(FC, FC_pred)
    print(f"  相关系数: {accuracy['correlation']:.3f}")
    print(f"  R²: {accuracy['r2']:.3f}")
    print(f"  RMSE: {accuracy['rmse']:.3f}")
    
    # 枢纽分类
    print("\n枢纽分类...")
    classifier = HubClassifier(config)
    classification = classifier.classify_hubs(SC, FC)
    print(f"  整合枢纽: {len(classification['integrator_hubs'])} 个")
    print(f"  中介枢纽: {len(classification['mediator_hubs'])} 个")
    
    # 虚拟损伤
    print("\n虚拟损伤实验...")
    lesion_exp = VirtualLesionExperiment(model)
    cascade = lesion_exp.analyze_cascade_effect(
        SC, FC,
        classification['integrator_hubs'][:3]
    )
    print(f"  最具破坏性区域: {cascade['most_disruptive_region']}")
    
    print(f"\n关键发现:")
    print(f"  ✅ 线性模型预测 FC")
    print(f"  ✅ 识别整合/中介枢纽")
    print(f"  ✅ 虚拟损伤预测级联效应")
    
    return model


## Activation Keywords
- 结构功能耦合
- 脑连接组
- 结构连接
- 功能连接
- 整合枢纽
- structure-function coupling
- structural connectivity
- functional connectivity
- integrator hub
- mediator hub

## Tools Used
- numpy
- scipy
- sklearn

## Instructions for Agents
1. 从 DWI 提取结构连接矩阵
2. 从 rs-fMRI 提取功能连接矩阵
3. 提取结构动机（直接/间接/三元）
4. 拟合线性生成模型
5. 识别整合枢纽和中介枢纽
6. 进行虚拟损伤实验预测级联效应

## Examples
```python
# 结构-功能耦合分析示例
from linear_structure_function_coupling import (
    LinearGenerativeModel, HubClassifier,
    VirtualLesionExperiment, StructureFunctionConfig
)

# 1. 配置
config = StructureFunctionConfig(n_regions=100)

# 2. 创建模型
model = LinearGenerativeModel(config)

# 3. 拟合
training_info = model.fit(SC, FC)

# 4. 预测
FC_pred = model.predict(SC)

# 5. 评估
accuracy = model.compute_prediction_accuracy(FC, FC_pred)
print(f"相关系数: {accuracy['correlation']:.3f}")

# 6. 枢纽分类
classifier = HubClassifier(config)
hubs = classifier.classify_hubs(SC, FC)

# 7. 虚拟损伤
lesion_exp = VirtualLesionExperiment(model)
cascade = lesion_exp.analyze_cascade_effect(SC, FC, hubs['integrator_hubs'])
```

if __name__ == "__main__":
    example_structure_function_coupling()
```

## Related Skills

- `brain-network-controllability` - 脑网络可控性
- `generative-brain-dynamics-models` - 脑动力学生成模型
- `time-varying-brain-connectivity` - 时变脑连接

## References

- arXiv:2507.06136 - A Linear Generative Framework for Structure-Function Coupling in the Human Brain
- Topics: Neurons and Cognition (q-bio.NC), Quantitative Methods (q-bio.QM)