Introduction
Deep brain stimulation (DBS) has emerged as a promising therapy for refractory psychiatric disorders like obsessive-compulsive disorder (OCD) and major depressive disorder (MDD), proving effective where conventional treatments fall short (Dougherty et al., Reference Dougherty, Brennan, Stewart, Wilhelm, Widge and Rauch2018; Lozano et al., Reference Lozano, Lipsman, Bergman, Brown, Chabardes, Chang, Matthews, McIntyre, Schlaepfer, Schulder, Temel, Volkmann and Krauss2019; Scangos et al., Reference Scangos, State, Miller, Baker and Williams2023). Yet substantial variability in treatment outcomes (Denys et al., Reference Denys, Graat, Mocking, de Koning, Vulink, Figee, Ooms, Mantione, van den Munckhof and Schuurman2020; Nora Runia et al., Reference Runia, Yücel, Lok, de Jong, Denys, van Wingen and Bergfeld2022; Xue et al., Reference Xue, Qiu, Li, Xie, Wang, Chen, Xie, Bai, Yang, Meng, Wang and Zhang2026) reveals the limitations of current trial-and-error protocols, highlighting the need for mechanistically informed, predictive strategies. Modern network neuroscience demonstrates that DBS acts not only as a focal intervention, but as a perturbation that modulates distributed brain dynamics underlying both pathology and treatment response (Cha et al., Reference Cha, Choi, Rajendra, McGrath, Riva-Posse, Holtzheimer, Figee, Kopell and Mayberg2024; Hollunder et al., Reference Hollunder, Ostrem, Sahin, Rajamani, Oxenford, Butenko, Neudorfer, Reinhardt, Zvarova, Polosan, Akram, Vissani, Zhang, Sun, Navratil, Reich, Volkmann, Yeh, Baldermann and Horn2024; Lynch et al., Reference Lynch, Elbau, Ng, Ayaz, Zhu, Wolk, Manfredi, Johnson, Chang, Chou, Summerville, Ho, Luecke, Bukhari, Buchanan, Victoria, Solomonov, Goldwaser, Moia and Liston2024; Meyer et al., Reference Meyer, Hollunder, Li, Butenko, Dembek, Hart, Nombela, Mosley, Akram, Acevedo, Borron, Chou, Castaño Montoya, Strange, Barcia, Tyagi, Castle, Smith, Choi and Horn2024), challenging its classical view as a focal intervention. This paradigm shift demands computational models that can simulate brain-wide network engagement (Cagnan, Denison, McIntyre, & Brown, Reference Cagnan, Denison, McIntyre and Brown2019; Kulkarni & Bassett, Reference Kulkarni and Bassett2025), while enabling not only patient-specific classification but also prediction of symptom trajectories before treatment.
A critical challenge in moving toward personalized DBS involves understanding how stimulation reshapes widely distributed neural dynamics. Past research has used electrophysiological recordings (e.g. local field potentials) to predict multiregional responses to time-varying stimulation in primate models (Schwock et al., Reference Schwock, Bloch, Khateeb, Zhou, Atlas and Yazdan-Shahmorad2025; Yang et al., Reference Yang, Qiao, Sani, Sedillo, Ferrentino, Pesaran and Shanechi2021). Diffusion MRI-derived structural connectomes have long served as fundamental scaffolds for these dynamic processes, informing biophysical models of pathological activity (e.g. epilepsy and neurodegeneration) and therapeutic targeting (Horn et al., Reference Horn, Reich, Ewert, Li, Al-Fatly, Lange, Roothans, Oxenford, Horn, Paschen, Runge, Wodarg, Witt, Nickl, Wittstock, Schneider, Mahlknecht, Poewe, Eisner and Kühn2022; Hutchings et al., Reference Hutchings, Han, Keller, Weber, Taylor and Kaiser2015; Lang et al., Reference Lang, Momi, Vetkas, Santyr, Yang, Kalia, Griffiths and Lozano2024; Ríos et al., Reference Ríos, Oxenford, Neudorfer, Butenko, Li, Rajamani, Boutet, Elias, Germann, Loh, Deeb, Wang, Setsompop, Salvato, de Almeid, Foote, Amaral, Rosenberg, Tang-Wai and Horn2022; Sanz Perl et al., Reference Sanz Perl, Fittipaldi, Gonzalez Campo, Moguilner, Cruzat, Fraile-Vazquez, Herzog, Kringelbach, Deco, Prado, Ibanez and Tagliazucchi2023). Network control theory combined with mean-field modeling has further clarified how anatomical connectivity governs emergent dynamics (Deco, Tononi, Boly, & Kringelbach, Reference Deco, Tononi, Boly and Kringelbach2015; Ju & Bassett, Reference Ju and Bassett2020; Meier et al., Reference Meier, Perdikis, Blickensdörfer, Stefanovski, Liu, Maith, Dinkelbach, Baladron, Hamker and Ritter2022; Parkes et al., Reference Parkes, Kim, Stiso, Brynildsen, Cieslak, Covitz, Gur, Gur, Pasqualetti, Shinohara, Zhou, Satterthwaite and Bassett2024). Notably, the relationship between brain connectivity and network-wide stimulation effects depends on the pre-stimulation dynamical state (Papadopoulos, Lynn, Battaglia, & Bassett, Reference Papadopoulos, Lynn, Battaglia and Bassett2020; Perl et al., Reference Perl, Zamora-Lopez, Montbrió, Monge-Asensio, Vohryzek, Fittipaldi, Campo, Moguilner, Ibañez, Tagliazucchi, Yeo, Kringelbach and Deco2023). Although directed and perturbation-based effective-connectivity approaches can provide information about causal interactions (Maran, Muller, & Fulcher, Reference Maran, Muller and Fulcher2025), the present study instead uses network deconvolution to estimate putative direct functional dependencies from undirected resting-state FC. Surrogate brain models using resting-state fMRI can now map connectivity changes under virtual perturbations (Luo et al., Reference Luo, Peng, Liang, Cai, Xu, Li, Hu, Zhou and Liu2025), and digital twin models integrate multimodal in vivo human data to simulate global dynamics and pinpoint rehab-critical regions (Zhang et al., Reference Zhang, Du, Li, Yang, Cao, Luo, Yang, Tang, Chu, Xiao, Li, Jiang, Wang, Du, Shi, Ma, Xiong, Song, Zhang and Jiang2025). Together, these approaches synthesize principles from statistical physics, information theory, and network neuroscience to bridge local circuits and whole-brain function (Kulkarni & Bassett, Reference Kulkarni and Bassett2025; Zhang et al., Reference Zhang, Liu, Liang, Cheng, Lou, Duan, Gao, Hu and Liu2026).
Here, we modeled the propagation of focal stimulation over a deconvolution-derived direct functional network estimated from whole-brain resting-state FC (Chen et al., Reference Chen, Zhang, Li, Huang, Lv, Yu, Chen, Sun and Wang2018). Under the assumed graph-diffusion model, network deconvolution attenuates higher-order indirect statistical dependencies and estimates a first-order network representation of putative direct functional dependencies (Feizi, Marbach, Médard, & Kellis, Reference Feizi, Marbach, Médard and Kellis2013; Lin, Pan, & Pan, Reference Lin, Pan and Pan2024). Moreover, we formulated the brain under regional stimulation as a linear dynamical system driven by time-varying oscillatory inputs, enabling systematic in silico perturbations across all regions and analysis of brain-wide responses. We retrospectively applied this framework to resting-state fMRI data from OCD and MDD patients (Figure 1), exploring an abstract perturbation parameter space defined by target, frequency, and amplitude. The best-fitting model parameter combination was selected by comparing simulated FC states with empirical DBS-ON fMRI-derived FC patterns. Importantly, we integrated this dynamic network modeling with machine learning to translate model-derived changes in functional connectivity into predictors of individual clinical outcomes. In addition, as a proof-of-concept for personalized therapeutic planning, we evaluated our framework’s ability to stratify individualized circuit-level treatment strategies in an independent cohort of OCD patients who have not yet received DBS therapy. By integrating whole-brain network modeling with clinically actionable insights, the present study contributes a mechanistic foundation for presurgical precision neuromodulation, with potential implications for informing therapeutic efficacy improvements in a wide range of brain disorders.
Conceptual schematic of the virtual neural modulation computational flowchart. (a) Network-based computational framework for virtual neuromodulation. Resting-state FC networks are first transformed using network deconvolution to obtain a deconvolution-derived first-order functional network. Linear dynamical modeling is then used to simulate how abstract virtual perturbations propagate through this network. The resulting simulated FC states are compared across candidate targets and model perturbation parameters to identify best-fitting network response patterns. (b) Experimental design. Patients first underwent DBS surgery and clinical parameter optimization. Electrodes were implanted in the nucleus accumbens (NAc) in OCD patients. In MDD patients, electrodes were implanted along the vALIC/NAc axis, with the deepest contact extending into the NAc and the three upper contacts located in the ventral anterior limb of the internal capsule (vALIC). Following optimization, patients entered a randomized crossover phase, in which the stimulator was set to DBS-ON or DBS-OFF. In each phase, patients underwent MRI scanning and concurrent clinical symptom assessments. (c) Individualized network features and clinical symptom scores were used to train models linking network weights to symptom variation, enabling patient-specific exploration of DBS-related clinical associations. (d) The model simplifies whole-brain resting-state dynamics as a first-order linear dynamical system. In this system, the signal observed at each node is treated as the result of diffusion through a deconvolution-derived functional network (d). Based on this framework, network deconvolution is used to estimate a first-order network representation of putative direct functional dependencies from the observed FC network by attenuating higher-order indirect statistical dependencies. After estimating this deconvolution-derived network, abstract periodic perturbation inputs, parameterized by target, frequency, and dimensionless model amplitude, are applied to selected target nodes to simulate DBS-related network responses within the linear dynamical model. (e) A series of simulated DBS-ON FC states are generated under different model perturbation parameter combinations. The network distance between each simulated DBS-ON FC state and the empirical DBS-ON FC state is then calculated to evaluate model fit. The best-fitting model parameter combination is selected based on the smallest network distance. Finally, this procedure is applied across multiple candidate stimulation targets, allowing model-based target ranking within the virtual perturbation framework.
Figure 1. Long description
The image is divided into five panels labeled A through E.
Panel A, Network-based computational framework: A flowchart starts with Network Topology showing a brain with nodes n 1 through n i and a connectivity matrix. This flows into a central box containing Network Deconvolution and Linear Dynamics. Arrows lead to a Virtual Perturbation stage where node n 1 is highlighted with a lightning bolt. The process concludes with a Generative Model involving individuals, different targets, and different stimulations, resulting in a Simulated Network with its own connectivity matrix.
Panel B, Experimental design: A timeline labeled Stage shows D B S Optimization with an anatomical cross-section of the N Ac. This is followed by a Crossover Phase with two branches: Active and Sham. Each branch includes M R I and Clinical Symptom assessment icons.
Panel C, Machine Learning: A brain diagram with personal features w 1, w 2, through w n is analyzed as a Biomarker. This leads to a Predict Symptom graph showing a positive linear correlation between individual data points and a dashed trend line.
Panel D, Linear dynamical system: This panel details the math behind the network. It shows an equation where observed connection F minus I equals direct connection D plus indirect connections represented as D squared (2 hops), D cubed (3 hops), and D to the 4th power (4 hops). This results in an inferred direct network matrix D. An oscillatory input with parameters for targets, frequencies, and amplitudes is applied to a linear dynamic system to produce the final Simulated Network.
Panel E, Model evaluation: On the left, an Empirical Network matrix is vectorized. On the right, a stack of Simulated Network matrices is vectorized across axes of frequencies and amplitudes. A central 3 D surface plot shows Euclidean Distance on the vertical axis against frequencies and amplitudes on the horizontal axes, marking a red dot for minimal distance. At the bottom, a series of brain maps for stimulation targets 1 through N are shown, ending in a bar chart sorting targets from low to high distance.
Empirical and modeled DBS effects in OCD dataset. (a) Whole-brain FC (standardized z-scores) matrices for DBS OCD and HC (stage off vs. stage on), excluding artifacts and noise regions. (b) FC z-value matrices (lower triangular) and effect size of change after ANOVA analysis (upper triangular) in the fronto-subcortical network. The left panel is the modeled DBS effect, and the right panel is the empirical DBS effects. Y-axis (vertical): Odd-numbered rows from bottom to top represent right hemisphere regions, while even-numbered rows correspond to left hemisphere regions; X-axis (horizontal): Odd-numbered columns from left to right denote left hemisphere areas, with even-numbered columns indicating right hemisphere areas. (c) Connectivity within the fronto-subcortical network showing positive (red) and negative (blue) group-by-treatment interaction effects at an uncorrected p < 0.05, with thicker lines indicating larger effect sizes. The left panel shows simulated DBS effects, and the right panel shows empirical DBS effects. Green rectangles represent frontal regions and purple rectangles represent subcortical regions. The stacked bar plots below show the proportions of these exploratory FC changes involving frontal (green) and subcortical (purple) regions. (d) Proportion of increased (red) and decreased (blue) FC changes among the connections shown in panel (c). The top panel shows cortico-subcortical FCs, and the bottom panel shows cortico-cortical FCs. Within each panel, the left column displays empirical DBS results, and the right column displays simulated DBS results. Statistical differences in the distributions were assessed using Fisher’s exact test.
Figure 2. Long description
Panel A: Four heatmaps arranged in a two-by-two grid show whole-brain F C standardized z-scores. The top row shows D B S O C D patients in Off and On stages. The bottom row shows H C in Off and On stages. Axes range from 1 to 64 regions. A color bar on the right indicates z-scores from negative 0.2 (brown) to 1.0 (orange).
Panel B: Two large matrices comparing simulated (left) and empirical (right) D B S effects in the fronto-subcortical network. The lower triangular portion shows F C z-values, while the upper triangular portion shows effect size (Cohen’s d) of change. Labels include P F C regions and subcortical areas like Cau, Put, Tha, N A c, and S T N. The simulated panel shows stronger positive effect sizes (red) compared to the empirical panel.
Panel C: Two circular chord diagrams showing connectivity. Frontal regions are marked with green perimeter rectangles and subcortical regions with purple. Red lines indicate positive interaction effects; blue lines indicate negative. The simulated panel (left) shows denser, thicker red connections than the empirical panel (right). Below each circle, a stacked bar plot shows the percentage of F C changes, with frontal regions occupying approximately 75 percent of the bar.
Panel D: Four donut charts comparing increased (red) and decreased (blue) F C changes. The top row (cortico-subcortical) shows empirical D B S at 10 percent increased versus simulated at 17 percent. The bottom row (cortico-cortical) shows empirical at 0 percent increased versus simulated at 33 percent. Brackets labeled n dot s dot indicate no significant statistical difference between empirical and simulated distributions.
Materials and methods
Study cohorts
This retrospective study analyzed three distinct participant groups: (1) an OCD baseline cohort (104 patients and 60 healthy controls [HCs]); (2) a treatment-refractory OCD cohort receiving deep brain stimulation (DBS; 14 patients and 15 controls); and (3) a treatment-resistant major depressive disorder (MDD) cohort receiving DBS (13 patients and 16 HCs). The baseline OCD data were previously published (Chen et al., Reference Chen, Wang, Lv, Lv, van Wingen, Fridgeirsson, Denys, Voon and Wang2022; Lv et al., Reference Lv, Lv, Yin, Zhang, Sun, Voon and Wang2021; Yin et al., Reference Yin, Zhang, Lv, Chen, Zeljic, Gong, Zhan, Jin, Wang and Sun2018; Zhao et al., Reference Zhao, Lv, Zhang, Zhang, Wang, Zhang, Wang and Wang2024). DBS-treated OCD patients were exclusively recruited from the Academic Medical Center (University of Amsterdam), while MDD patients were additionally enrolled from St. Elisabeth Hospital (Tilburg); all patient data were derived from prior publications (Bergfeld et al., Reference Bergfeld, Mantione, Hoogendoorn, Ruhé, Notten, van Laarhoven, Visser, Figee, de Kwaasteniet, Horst, Schene, van den Munckhof, Beute, Schuurman and Denys2016; Figee et al., Reference Figee, Luigjes, Smolders, Valencia-Alfonso, Van Wingen, De Kwaasteniet, Mantione, Ooms, De Koning, Vulink, Levar, Droge, van den Munckhof, Schuurman, Nederveen, van den Brink, Mazaheri, Vink and Denys2013; Fridgeirsson et al., Reference Fridgeirsson, Figee, Luigjes, van den Munckhof, Schuurman, van Wingen and Denys2020). Clinical assessments included the Yale-Brown Obsessive Compulsive Scale (Y-BOCS) and the 17-item Hamilton Depression Rating Scale (HAM-D). All data were obtained from electronic databases and were fully anonymized before analysis. The study protocol was approved by the Medical Ethics Committee of the Academic Medical Center, University of Amsterdam.
Participants with excessive head motion (>3 mm translation displacement) or MRI abnormalities were excluded, resulting in the removal of 12 baseline OCD patients, 1 baseline HC, 3 DBS-OCD patients, and 1 DBS-OCD control. For the MDD cohort, which required MRI scans during both DBS states, we excluded two patients and two controls. The final analyzed samples were: 92 baseline OCD patients and 59 HCs; 11 DBS-OCD patients and 14 HCs; and 11 DBS-MDD patients and 14 HCs (see Table 1).
Demographic and clinical characteristics of subjects
Table 1. Long description
The first table focuses on O C D and contains three columns: D B S O C D (N = 11), D B S H C (N = 14), and Local O C D (N = 92).
* Age (years): 45.0 plus or minus 10.8, 43.7 plus or minus 9.2, and 30.8 plus or minus 9.6 respectively.
* Gender (males/females): 7/4, 7/7, and 56/36.
* Disease duration (year): 25.8 plus or minus 14.2 for D B S O C D and 9.7 plus or minus 6.2 for Local O C D.
* Y dash B O C S: 29.0 plus or minus 6.8 for Local O C D.
* Y dash B O C S (D B S dash O F F): 28.5 plus or minus 6.3 for D B S O C D.
* Y dash B O C S (D B S dash O N): 21.5 plus or minus 6.7 for D B S O C D.
* D B S frequency: 185 Hz (N = 2) and 130 Hz (N = 9) for the D B S O C D group.
* D B S amplitude: 4.8 plus or minus 0.7 for D B S O C D.
The second table focuses on M D D and contains two columns: D B S M D D (N = 11) and D B S H C (N = 14).
* Age (years): 50.9 plus or minus 7.8 and 53.7 plus or minus 9.5.
* Gender (males/females): 4/7 and 6/8.
* Disease duration (year): 20.7 plus or minus 12.7 for D B S M D D.
* H A M dash D dash 17 (D B S dash O F F): 22.9 plus or minus 4.7 for D B S M D D.
* H A M dash D dash 17 (D B S dash O N): 15.8 plus or minus 7.4 for D B S M D D.
* D B S frequency: 180 Hz (N = 5), 160 Hz (N = 2), and 130 Hz (N = 4) for the D B S M D D group.
* D B S amplitude: 4.6 plus or minus 0.6 for D B S M D D.
Note: HAM-D-17, Hamilton Depression Rating Scale; Y-BOCS, Yale–Brown Obsessive Compulsive Scale. Values are mean ± standard deviation.
DBS intervention protocols
The detailed DBS intervention protocols have previously been published and briefly summarized here. The treatment-refractory OCD cohort underwent nucleus accumbens (NAc) DBS for ≥1 year to achieve stable clinical improvement (DBS-ON phase), followed by a 1-week stimulation discontinuation (DBS-OFF phase). Resting-state fMRI and clinical assessments were acquired under both the stable DBS-ON condition and the DBS-OFF condition after stimulation discontinuation. Stimulation parameters included a mean voltage of 4.8 V (range: 3.5–6.2 V), frequencies of 130 or 185 Hz, and pulse widths of 90 or 150 μs (Figee et al., Reference Figee, Luigjes, Smolders, Valencia-Alfonso, Van Wingen, De Kwaasteniet, Mantione, Ooms, De Koning, Vulink, Levar, Droge, van den Munckhof, Schuurman, Nederveen, van den Brink, Mazaheri, Vink and Denys2013; Fridgeirsson et al., Reference Fridgeirsson, Figee, Luigjes, van den Munckhof, Schuurman, van Wingen and Denys2020).
The treatment-refractory MDD cohort received bilateral DBS implanted with 4-contact electrodes, with three proximal contacts positioned in the ventral anterior limb of the internal capsule (vALIC) and the distal contact extending into the nucleus accumbens (NAc). Following an optimization period (up to 52 weeks) to establish a 4-week stable response, patients participated in a double-blind, randomized crossover trial with 1-week DBS-ON and DBS-OFF phase. Resting-state fMRI and clinical assessments were acquired during both stimulation conditions. Parameters included voltages of 2.5–6.0 V, frequencies of 130, 160, or 180 Hz, and a 90 μs pulse width (Bergfeld et al., Reference Bergfeld, Mantione, Hoogendoorn, Ruhé, Notten, van Laarhoven, Visser, Figee, de Kwaasteniet, Horst, Schene, van den Munckhof, Beute, Schuurman and Denys2016).
MRI data preprocessing
Resting-state fMRI data were preprocessed using SPM12 (http://www.fil.ion.ucl.ac.uk/spm) and the CONN toolbox (https://web.conn-toolbox.org/). The preprocessing pipeline was designed following previously published methodologies to minimize artifacts and noise introduced by DBS devices (Chen et al., Reference Chen, Wang, Lv, Lv, van Wingen, Fridgeirsson, Denys, Voon and Wang2022). Preprocessing steps included discarding the first 10 volumes, slice timing correction, motion correction, spatial co-registration to the MNI space, resampling in a 4 mm isotropic grid, and spatial smoothing (full width at half maximum = 2 mm). Further denoising steps included linear detrending, regressing out Friston-12 motion parameters (original and derivative of the translation and rotation values estimated during motion correction), regressing out principal components extracted from CSF and white matter signals by a CompCor method (Behzadi, Restom, Liau, & Liu, Reference Behzadi, Restom, Liau and Liu2007; Fridgeirsson et al., Reference Fridgeirsson, Figee, Luigjes, van den Munckhof, Schuurman, van Wingen and Denys2020), and temporal filtering (0.01–0.08 Hz).
Functional network construction
We constructed the functional connectivity (FC) network using a customized brain atlas comprising 94 regions, including cortical regions parcellated from the Regional Map atlas (Bezgin et al., Reference Bezgin, Vakorin, van Opstal, McIntosh and Bakker2012; Zhan et al., Reference Zhan, Wei, Liang, Xu, He, Robbins and Wang2021), and subcortical regions parcellated from the FreeSurfer atlas (Fischl et al., Reference Fischl, Salat, Busa, Albert, Dieterich, Haselgrove, Van Der Kouwe, Killiany, Kennedy, Klaveness, Montillo, Makris, Rosen and Dale2002). The subthalamic nucleus (STN) was additionally defined using the ATAG template (Keuken et al., Reference Keuken, Bazin, Crown, Hootsmans, Laufer, Müller-Axt, Sier, van der Putten, Schäfer, Turner and Forstmann2014). The mean time series for each region of interest was obtained by averaging the BOLD signal across all constituent voxels. We computed FC between all regional pairs using Pearson correlation coefficients of these time series. All correlation coefficients underwent Fisher’s z-transformation to normalize their distribution for subsequent analyses. Following our established methodology (Chen et al., Reference Chen, Wang, Lv, Lv, van Wingen, Fridgeirsson, Denys, Voon and Wang2022), we focused on 32 predefined regions (22 prefrontal and 10 subcortical) within cortico-striato-thalamo-cortical (CSTC) circuits, as detailed in Supplementary Table S1. This approach yielded 496 functional connections while simultaneously addressing several technical challenges, including incomplete scanning coverage and DBS-induced metal artifacts.
For the MDD dataset, we applied similar preprocessing procedures with two modifications. We excluded 20 regions using identical quality control criteria as for the OCD data, and then expanded our analysis focus to 30 regions relevant to the core circuitry underlying MDD neuropathology, including the default mode network (DMN), ventral attention network (VAN), and subcortical regions (listed in Supplementary Table S2) (Drysdale et al., Reference Drysdale, Grosenick, Downar, Dunlop, Mansouri, Meng, Fetcho, Zebley, Oathes, Etkin, Schatzberg, Sudheimer, Keller, Mayberg, Gunning, Alexopoulos, Fox, Pascual-Leone, Voss and Liston2017; Lynch et al., Reference Lynch, Elbau, Ng, Ayaz, Zhu, Wolk, Manfredi, Johnson, Chang, Chou, Summerville, Ho, Luecke, Bukhari, Buchanan, Victoria, Solomonov, Goldwaser, Moia and Liston2024), which resulted in 435 functional connections. Each DMN and VAN region’s functional network assignment was determined using the Yeo7 parcellation scheme (Yeo et al., Reference Yeo, Krienen, Sepulcre, Sabuncu, Lashkari, Hollinshead, Roffman, Smoller, Zöllei, Polimeni, Fischl, Liu and Buckner2011), based on the proportional overlap of voxels with established network boundaries. Following prior studies (Bergfeld et al., Reference Bergfeld, Mantione, Hoogendoorn, Ruhé, Notten, van Laarhoven, Visser, Figee, de Kwaasteniet, Horst, Schene, van den Munckhof, Beute, Schuurman and Denys2016; Fridgeirsson et al., Reference Fridgeirsson, Bergfeld, de Kwaasteniet, Luigjes, van Laarhoven, Notten, Beute, van den Munckhof, Schuurman, Denys and van Wingen2025), we used the NAc as the DBS-targeted area to perform FC network analyses in both simulated and empirical data. Because the present framework is based on gray-matter resting-state FC and does not explicitly model white-matter stimulation fields or fiber activation, we represented the intervention as a node-level perturbation of the bilateral NAc, corresponding to the ventral striatal gray-matter component of the vALIC/NAc target. This approximation was used for modeling purposes and should not be interpreted as a reconstruction of the full vALIC stimulation volume.
Identifying the empirical network effects induced by DBS
Given the double-blind, crossover design of the DBS patient and matched control datasets, we employed repeated-measures ANOVA to assess DBS-induced effects. The model included two fixed factors: group (patient vs. control) and treatment condition (DBS-OFF vs. DBS-ON). We specifically examined the group × treatment interaction effect to identify differential neural responses to DBS between patients and controls. Analyses were conducted using MATLAB’s ranova function. Edge-wise group-by-treatment interaction effects were initially evaluated using an uncorrected threshold of p < 0.05 for exploratory visualization. To account for multiple comparisons, Benjamini–Hochberg false discovery rate (FDR-BH) correction was subsequently applied separately across all connections in the OCD and MDD data, with statistical significance defined as pFDR < 0.05. Effect sizes for significant FCs were estimated using Cohen’s d. For both the OCD and MDD datasets, we classified significant FCs based on their anatomical distribution – comparing cortico-cortical versus cortico-subcortical changes – and assessed their spatial patterning using Fisher’s exact test (significance threshold: p < 0.05).
Modeling how brain signals spread and respond to DBS
To understand how DBS affects the brain, we first need a clear map of how brain regions “talk” to one another. Standard functional connectivity networks are often “blurry” because they capture both direct conversations between two regions and “echoes” or indirect signals passing through multiple other areas. To estimate direct neural interactions underlying observed FC, we employed a linear dynamic system model. In this framework, resting-state fMRI signals are modeled as a first-order linear process:
where x represents brain regional activity over time, w is Gaussian noise, and D is the direct network capturing signal transitions between regions. Under this model, FC reflects the cumulative effects of both direct and indirect interactions among brain regions. To recover the direct network D from the observed FC matrix F, we adopted a network deconvolution approach, which strips away these “polysynaptic” or indirect echoes, leaving us with an effectively sparser weighted map of the direct “highways” between brain regions. This sparsification is reflected in a more concentrated edge-weight distribution, with weaker connections attenuated relative to stronger connections. Consistently, the Gini coefficient of absolute edge weights increased after deconvolution, indicating greater concentration of connection weights. A representative comparison of FC before and after deconvolution in the OCD-DBS dataset is provided in Supplementary Figure S1. Importantly, this model-derived network is undirected and should not be interpreted as causal effective connectivity, anatomical connectivity, or evidence of monosynaptic interactions between brain regions.
Once the deconvolution-derived functional network was obtained, we simulated DBS-related network responses by applying an abstract periodic perturbation to a selected target node, such as the NAc. This periodic input was used as a mathematical forcing term to approximate sustained external perturbation within the discrete-time dynamical system, rather than to reproduce the physical waveform of clinical DBS. By applying this perturbation to the deconvolution-derived network, we calculated how the modeled input propagated through the network and produced a simulated DBS-ON FC state. This allowed us to explore how an individual’s baseline functional architecture may constrain DBS-related network-state changes under different model perturbation settings, as illustrated in Figure 1d. The full mathematical derivation is provided in the Supplementary Materials.
Individualized search of best-fitting model perturbation parameters
To identify the model parameter combination that best reproduced each individual’s empirical DBS-ON FC state, we searched an abstract perturbation parameter space of the discrete-time dynamical model. The model frequency was sampled across the normalized angular-frequency interval [0, π], and the model amplitude ranged from 0 to 10 in dimensionless arbitrary units. These parameters should not be interpreted as clinically programmable DBS settings and do not have a direct one-to-one mapping to physical stimulation frequency, voltage, current, or pulse width. Accordingly, clinical DBS frequencies, such as 130, 180, or 185 Hz, cannot be directly compared with the model frequency range. For OCD and MDD, the virtual stimulation target was initially set to bilateral NAc. We performed a greedy grid search over normalized model frequency with a step size of 0.02π and dimensionless model amplitude with a step size of 0.1. Each combination of model frequency and amplitude was used to generate a simulated DBS-ON FC state (example shown in Supplementary Figure S2 and Supplementary Movie S1). The best-fitting model parameter combination was defined as the one yielding the minimum Euclidean distance between the simulated FC state and the empirical DBS-ON FC profile (Supplementary Figure S3 for OCD cohorts and Supplementary Figure S4 for MDD cohorts). This procedure was used to explore the model response space and estimate individualized best-fitting perturbation patterns (Figure 1e), rather than to optimize clinical DBS device parameters.
Simulated DBS effects were validated using Euclidean network distance comparisons between (1) simulated and empirical DBS-ON FCs and (2) empirical DBS-OFF versus DBS-ON FCs. Statistical significance was assessed via Wilcoxon signed-rank tests. To evaluate target specificity, we compared network distances for clinical targets (NAc) versus alternative sites. For topological analysis, we generated null models using the GRETNA toolbox, preserving edge-level topology while randomizing baseline FCs. Each of 1,000 randomized FCs underwent DBS simulation, with the resulting null distribution used to quantify model fit.
To reduce feature dimensionality, candidate connections were restricted to the strongest 10% of deconvolution-derived connections within predefined disorder-relevant networks, including the CSTC circuitry for OCD and DMN, VAN, and subcortical networks for MDD. Within these constrained feature spaces, one-component partial least squares regression (PLSR) was used to characterize associations between DBS-ON FC (Fisher z-transformed) and clinical symptom scores. A single PLS component was used to maintain model parsimony given the limited sample size. Multivariate FC signatures were identified based on cross-validated predictive performance. Predictive performance was evaluated using leave-one-out cross-validation (LOOCV), with feature standardization and model fitting performed using training subjects only. Performance was quantified by the Pearson correlation between predicted and observed symptom scores. To assess whether simulated FC expressed the same symptom-related pattern, PLSR was refitted to the empirical DBS-ON data using the same FC signature and applied to each individual’s simulated FC. The association between model-estimated and observed clinical scores was then assessed using Pearson correlation (Figure 3d and Supplementary Figure S5).
Model simulation performance and prediction of DBS treatment outcome. (a) Cortical changes based on the Regional Map atlas in FC within the fronto-subcortical network (seeded at the left and right NAc, which are the actual DBS stimulation sites). The left side of the colorbar (Cohens’ d, normalized to −0.5 and 0.5) indicates the modeled DBS results, while the right side corresponds to the empirical results. (b) The null distribution of the test statistic using a randomized DBS-OFF network as an input to the model. The vertical dashed line indicates the observed value based on the actual DBS-OFF network. (c) Simulated DBS-ON to empirical DBS-ON network distances by stimulating different targets, averaged over all patients. The location of NAc is highlighted in red, and the dot represents the mean distance at the region group level, including the subcortical, frontal, temporal, insula, occipital, and parietal regions. (d) Pipeline to build the symptom indicator based on partial least squares regression (PLSR) and to predict the DBS treatment outcome. (e) Association between observed DBS-ON Y-BOCS scores and model-estimated Y-BOCS scores derived from simulated FC (Pearson’s r = 0.465, 95% CI [−0.187, 0.833], p = 0.149).
Figure 3. Long description
Panel A shows brain surface maps comparing Simulated D B S and Empirical D B S results for N Ac L and N Ac R sites. A color bar for Cohens’ d ranges from negative 0.5 in blue to 0.5 in red. Panel B is a distribution plot of Euclidean distance (Model minus O N) on the x-axis from 5.0 to 7.0. A dashed red line at approximately 5.05 marks the observed value, while a null distribution of orange dots and a density curve are centered around 6.8. Panel C is a bar chart showing Euclidean distance on the y-axis (5.00 to 5.25) across various brain regions grouped by color: purple for subcortical, green for frontal, pink for temporal, light blue for insula, brown for occipital, and grey for parietal. The N Ac region is highlighted in red text. Panel D illustrates a Partial Least Squared Regression pipeline. It starts with a box containing Empirical D B S O N Y dash B O C S (y sub 1 to y sub n) and Empirical D B S O N F C (z) matrices. An arrow points to a beta coefficient vector, which is applied to a brain model with green frontal nodes and purple subcortical nodes connected by red positive and blue negative lines. This leads to a Simulated D B S O N F C (z) matrix. Panel E is a scatter plot showing the association between Empirical D B S O N Y dash B O C S on the x-axis and Simulated D B S O N Y dash B O C S on the y-axis. Data points follow a positive linear trend with a Pearson’s r of 0.465, indicated by a dashed line and a shaded 95 percent confidence interval.
Model validation of empirical and modeled DBS effects in MDD dataset. (a) FC z-value matrices (lower triangular) and ANOVA-derived effect sizes of DBS-induced changes (upper triangular) for the ventral attention network, the default mode network, and subcortical regions. Network boundaries are separated by dashed lines. The left panel represents the simulated DBS effects, while the right panel depicts the empirical DBS effects. (b) Correlation of effect size values between empirical and simulated outcomes within the ventral attention network, the default mode network, and subcortical regions, and across each network. ***p < 0.001. (c) Connectivity within the large-scale network showing positive (red) and negative (blue) group-by-treatment interaction effects at an uncorrected p < 0.05, with thicker lines indicating larger effect sizes. The left panel shows simulated DBS effects and the right panel shows empirical DBS effects. Red rectangles represent default mode network regions, violet rectangles represent ventral attention network regions, and purple rectangles represent subcortical regions. The stacked bar plots below show the proportions of these exploratory FC changes across the corresponding network categories. (d) Proportion of increased (red) or decreased (blue) FCs for the connections shown in panel (c), grouped into cortico-subcortical (upper) or cortico-cortical (bottom) connectivity, with the left panel (empirical DBS results) and the right panel (simulated DBS results). Statistical differences in the distributions were assessed using Fisher’s exact test.
Figure 4. Long description
Panel A: Two heatmaps showing F C z-values in the lower triangular and ANOVA effect sizes in the upper triangular. The left heatmap represents simulated D B S effects and the right represents empirical D B S effects. Regions include I a dot R, I p dot R, C C a dot R, P H C dot R, H C dot R, C C p dot R, P F C p o l dot R, P C m dot R, P F C d m dot R, P F C d l dot R, C a u dot R, P u t dot R, T h a dot R, N A c dot R, and G P dot R.
Panel B: A correlation matrix showing effect size values between empirical and simulated outcomes. Three orange-shaded squares with triple asterisks indicate p < 0.001. Color-coded bars represent Subcortical Regions (light blue), Default Mode Network (red), and Ventral Attention Network (violet).
Panel C: Two circular chord diagrams showing connectivity changes. The left circle shows simulated D B S effects and the right shows empirical D B S effects. Red lines indicate positive interaction effects and blue lines indicate negative effects. Thicker lines represent larger effect sizes. Below each circle is a stacked bar plot showing the percentage of F C changes across the three network categories.
Panel D: Four donut charts comparing empirical D B S (left) and simulated D B S (right). The top row shows Cortico-subcortical connectivity: empirical is 50 percent increased (red) and 50 percent decreased (blue); simulated is 33 percent increased and 67 percent decreased. The bottom row shows Cortico-cortical connectivity: empirical is 73 percent increased and 27 percent decreased; simulated is 100 percent increased and 0 percent decreased. Both rows are marked n dot s dot for non-significant differences.
Model simulation performance and prediction of DBS treatment outcome in MDD dataset. (a) Cortical changes based on the Regional Map atlas in FC within the large-scale networks (seeded at the left and right NAc, which is the actual DBS stimulation sites). The left side of the colorbar (Cohens’ d, normalized to −0.5 and 0.5) indicates the simulated DBS results, while the right side corresponds to the empirical results. (b) The null distribution of the test statistic using a randomized DBS-OFF network as an input to the model. The vertical dashed line indicates the observed value based on the actual DBS-OFF network. (c) Simulated DBS-ON to empirical DBS-ON network distances by stimulating different targets, averaged over all patients. The location of NAc is highlighted in red, and the dot represents the mean distance at the region group level, including the subcortical, frontal, temporal, occipital, parietal, and insula regions. (d) FC with beta value derived from PLSR analysis. The connectivity thickness and transparency correspond to the magnitude of the beta values, indicating the strength of FC’s contribution. (e) Association between observed DBS-ON HAM-D scores and model-estimated HAM-D scores derived from simulated FC (Pearson’s r = 0.627, 95% CI [0.043, 0.892], p = 0.039).
Figure 5. Long description
Panel A: A series of brain surface maps comparing Simulated D B S and Empirical D B S results. The maps show cortical changes in F C seeded at the left and right N Ac. A color bar for Cohens’ d ranges from negative 0.5 in blue to 0.5 in red. Simulated results show localized blue clusters, while empirical results show more widespread red and blue activations across the cortex.
Panel B: A distribution plot of Euclidean distance between Model and O N states. An orange bell curve and scatter of points are centered around 6.9. A vertical dashed red line marks the observed value at approximately 4.8, significantly to the left of the null distribution.
Panel C: A bar chart showing Euclidean distance for different brain regions. The y-axis ranges from 4.70 to 5.00. Bars are color-coded by region: Subcortical in purple, Frontal in green, Temporal in pink, Occipital in maroon, Parietal in grey, and Insula in light blue. N Ac is highlighted in red text. A dotted line connects mean dots for each group, showing an upward trend from subcortical to insula regions.
Panel D: Two transparent brain models showing functional connectivity. Red lines indicate positive beta values and blue lines indicate negative values. Nodes are categorized as Default Mode Network nodes in red, Ventral Attention Network nodes in purple, and Subcortical nodes in light blue.
Panel E: A scatter plot correlating Empirical D B S-O N H A M-D scores on the x-axis with Simulated D B S-O N H A M-D scores on the y-axis. Orange dots show a positive linear trend with a Pearson’s r of 0.627. A dashed regression line is surrounded by a light orange shaded confidence interval.
Clinical variables and potential confound
Illness duration, baseline symptom severity, and clinical stimulation settings were summarized descriptively where available (Table 1). Medication status was not consistently available across cohorts. These variables were not entered as covariates in the repeated-measures ANOVA, model-fitting procedure, or PLSR-based symptom association analyses. Given the small DBS sample sizes and incomplete medication-level information, the present study was not sufficiently powered to perform reliable covariate-adjusted or sensitivity analyses for these potential confounds.
External proof-of-concept simulation in an independent OCD dataset
To evaluate generalizability, we applied our framework to an independent local dataset of OCD patients who had not undergone DBS treatment. A resting-state FC network was computed for each patient, and virtual stimulation targeting the NAc was simulated using the same model-based parameter search described earlier. This procedure generated individualized DBS-ON FCs across all combinations of frequency and amplitude. Stimulation parameters were personalized by minimizing the Euclidean distance between each simulated FC and the average empirical DBS-ON FCs derived from 11 DBS-treated OCD patients. The parameter combination yielding the minimum Euclidean distance between the simulated and empirical DBS-ON FC profiles was identified as the individual’s best-fitting model perturbation pattern. Consistent with clinical practice, simulated stimulation was restricted to the bilateral NAc. The resulting FCs were then input into our pretrained PLSR model – previously trained on empirical DBS-ON FCs and corresponding clinical scores – to predict DBS-ON phase Y-BOCS symptom severity (Figure 6a).
Simulation of individualized NAc-targeted DBS effects using different datasets as a proof-of-concept. (a) Schematic overview of the analytical pipeline. Resting-state FC data from a local OCD cohort (n = 92) were first obtained. Using the same parameter search procedure as described in the main framework, individualized FCs under simulated NAc stimulation were generated. These simulated FCs were then integrated with the previously derived OCD PLSR biomarker (β map) to predict DBS-ON phase symptom severity (Y-BOCS scores) for each patient. (b) Permutation analysis comparing the average Euclidean distance of model-simulated NAc stimulation effects in three model-prioritized candidates against distributions derived from randomly selected triplets of patients (n = 1,000 permutations, total N = 92). The observed mean distance for the three candidates (red dashed line) was significantly lower than expected by chance (p = 0.003). (c) For the three model-prioritized candidates, the average simulated FC patterns (lower triangle) and the mean empirical DBS-ON FC from the 11 DBS-OCD patients (upper triangle) are shown. The bar plot on the right displays the corresponding correlation coefficients between network pairs (p < 0.001, Bonferroni-corrected). Network colors: green, Frontal regions; light purple, Subcortical regions. (d) Group-level comparison of simulated DBS-ON to empirical DBS-ON network distances across stimulation targets for the three model-prioritized candidates. Each dot represents the mean Euclidean distance between simulated and empirical DBS-ON networks for a given target, averaged across the three candidates. The NAc is highlighted in red. Regional groupings include subcortical, frontal, temporal, occipital, insular, and parietal areas. (e) Model-estimated symptom outcomes for all 92 patients. Upper: paired baseline and model-estimated DBS-ON Y-BOCS scores per individual. Lower: model-estimated DBS-ON scores versus model-estimated percentage reductions from baseline. Red points indicate the three model-prioritized candidates. The dashed red line denotes the conventional >35% Y-BOCS responder threshold used in prior OCD DBS studies. Because no DBS was performed in this cohort, this threshold is solely a benchmark for model-based candidate prioritization and does not classify observed clinical response.
Figure 6. Long description
Panel A: A flowchart titled Analytical Pipeline. It starts with a Local O C D Cohort of N equals 92 individuals represented by orange icons. Arrows point to a series of heatmaps representing functional connectivity matrices. These are processed through a simulation step involving 3D surface plots of distance across frequencies and amplitudes, leading to an apply step where a brain model with green nodes and colored edges represents a biomarker.
Panel B: A histogram showing the distribution of Distance Difference from negative 0.6 to 1.0. A red dashed line at approximately negative 0.5 indicates an observed p-value of 0.0030.
Panel C: A square functional connectivity matrix heatmap. The upper triangle shows empirical D B S-O N data and the lower triangle shows simulated data. A color scale at the top ranges from brown (negative 0.1) to orange (0.8). Right-side bars indicate correlation values of 0.82 for Frontal regions (green) and 0.48 for Subcortical regions (light purple).
Panel D: A bar chart of Euclidean distance across various brain targets. Targets are grouped by color: Subcortical (purple), Frontal (green), Temporal (pink), Occipital (brown), Insula (blue), and Parietal (gray). The N A c target is highlighted in red within the subcortical group.
Panel E: Two clinical outcome plots. The top plot is a slope chart connecting Baseline Y-B O C S scores (ranging 10 to 40) to Simulated D B S-O N Y-B O C S scores. Most lines show a downward trend. The bottom plot is a scatter plot showing Simulated D B S-O N Y-B O C S on the left Y-axis and Percentage improvement on the right Y-axis. A red dashed line marks the 35 percent improvement threshold. Three prioritized candidates are highlighted as red dots across both plots.
Because these independent OCD patients had not undergone DBS, no observed post-operative clinical outcomes were available. The estimated clinical improvement was therefore derived entirely from the model-estimated DBS-ON Y-BOCS score. Specifically, the model-estimated percentage improvement was calculated relative to each individual’s baseline Y-BOCS score as (Y-BOCS_baseline – Y-BOCS_model-estimated DBS-ON)/Y-BOCS_baseline × 100%. The >35% threshold was adopted as a conventional OCD DBS responder benchmark from prior literature; however, in this cohort, it served only as a reference for model-based candidate prioritization and should not be interpreted as an observed treatment response.
To further examine the model-based specificity of NAc stimulation, we performed a permutation-based ranking analysis (n = 1,000 iterations), identifying individuals whose simulated NAc stimulation consistently ranked among the top candidate targets (individual regional rankings were shown in Supplementary Figure S6). For these selected subjects, the average modeled FC pattern was correlated with the mean empirical DBS-ON FC from the 11 DBS-treated patients to evaluate network modulation similarity. Resulting p-values underwent Bonferroni correction for multiple comparisons. Additionally, predicted Y-BOCS scores were used to quantify expected clinical improvement for these individuals.
Results
Model-simulated and empirical DBS network effects in OCD
For the OCD cohort, we computed 64 × 64 FC matrices for both DBS and HC groups across DBS-OFF and DBS-ON states (Figure 2a). A repeated-measures ANOVA compared model-simulated FC changes (Figure 2b, left) against empirical DBS-induced changes (Figure 2b, right). At an exploratory, uncorrected threshold of p < 0.05, 28 connections in the model-simulated data (Figure 2c, left) and 17 connections in the empirical data (Figure 2c, right) showed group-by-treatment interaction effects (Figure 2c). However, no individual connection survived FDR-BH correction at pFDR < 0.05 in either analysis. Accordingly, the following edge-level findings are presented as descriptive, hypothesis-generating patterns rather than confirmatory evidence of specific DBS-modulated connections. Among the uncorrected empirical effects, 70.59% involved frontal regions and 29.41% involved subcortical regions; the corresponding proportions in the simulated data were 75% and 25%, respectively. The uncorrected effects also showed differing directional patterns across subnetworks: 9 of 10 cortico-subcortical connections decreased, whereas all 7 cortico-cortical connections increased in the empirical data (Figure 2d, “Empirical DBS”). In the simulated data, 5 of 6 cortico-subcortical connections decreased, whereas 12 of 18 cortico-cortical connections increased (Figure 2d, “Simulated DBS”). Fisher’s exact tests did not detect significant differences between the empirical and simulated directional distributions (p > 0.05). Four connections showed overlapping uncorrected effect patterns between the simulated and empirical analyses, including decreased connectivity between NAc.L–PFCm.L, NAc.L–CCs.L, and Cau.L–PFCcl.L, and increased connectivity between PFCcl.R–PFCol.R. Spatial congruence was further supported by strong correlations between simulated and empirical DBS network effects for NAc.L (r = 0.718, p < 0.001) and NAc.R (r = 0.808, p < 0.001) when projecting effect sizes onto cortical surfaces (Figure 3a).
To quantify global alignment, we computed the Euclidean network distance between simulated and empirical DBS-ON states. Simulations closely matched empirical data (Wilcoxon signed-rank test: p < 0.001). To test whether the model leveraged biological network structure, we compared its performance to a topological null model (matched for edge count, amplitudes, and node degrees). The true model significantly outperformed the null (p < 0.001; Figure 3b), confirming that DBS effects emerged from biologically constrained topology. Finally, aligning with clinical benchmarks, stimulation of the NAc – the clinically validated target – yielded significantly better fits than alternative regions (Figure 3c), reinforcing the model’s validity in target-specific prediction.
Virtually predicted and empirical clinical outcomes of DBS in OCD
To assess individual-level clinical outcome associations, we adopted a two-step framework. First, the resulting multivariate FC signature comprised 15 connections, predominantly involving frontal connectivity (Figure 3d). Under LOOCV, predicted and observed DBS-ON Y-BOCS scores were positively correlated (Pearson’s r = 0.625, 95% CI [0.041, 0.891], p = 0.040; R 2 = 0.379). Next, the PLSR model fitted to the empirical DBS-ON data was applied to the corresponding connections in the individually simulated FC derived from DBS-OFF network states. The model-estimated and observed Y-BOCS scores remained positively associated, although the correlation did not reach statistical significance (r = 0.465, 95% CI [−0.187, 0.833], p = 0.149; R 2 = 0.14; Figure 3e). Additional permutation analyses provided further support for the stability of the empirical FC–symptom association relative to randomized clinical labels and randomly sampled edge sets from the same candidate feature space (Supplementary Figure S5a). Together, these findings suggest that empirical DBS-ON FC patterns contained information relevant to individual symptom prediction under LOOCV, whereas predictive performance based on simulated FC was more limited and did not reach statistical significance.
Model-simulated and empirical DBS network effects in MDD
To characterize DBS-induced network changes in MDD, we analyzed 30 × 30 functional connectivity matrices using repeated-measures ANOVA, comparing model-simulated effects with empirical measures (Figure 4a,b). Network-specific effects were evaluated using the Yeo 7-network parcellation, with particular focus on the DMN and VAN. Pearson correlations (Bonferroni-corrected) of Cohen’s d effect sizes demonstrated significant agreement between simulated and empirical ΔFC across most networks (Figure 4b). At an exploratory, uncorrected threshold of p < 0.05, 18 connections showed group-by-treatment interaction effects in both the simulated and empirical analyses (Figure 4c). However, no individual connection survived FDR-BH correction at pFDR < 0.05 in either analysis. Accordingly, the following edge-level findings are presented as descriptive, hypothesis-generating patterns rather than confirmatory evidence of specific DBS-modulated connections. Among the uncorrected empirical effects, 61.1% involved the DMN, 27.8% involved subcortical regions, and 11.1% involved the VAN; the corresponding proportions in the simulated data were 41.7%, 30.6%, and 27.7%, respectively. The uncorrected effects also showed similar directional patterns across connection types: 50% of the cortico-subcortical connections decreased in the empirical data, compared with 67% in the simulated data, whereas 73% of the cortico-cortical connections increased in the empirical data, compared with 100% in the simulated data. Fisher’s exact tests did not detect significant differences between the empirical and simulated directional distributions (p > 0.05; Figure 4d). The overlapping uncorrected exploratory patterns between the simulated and empirical analyses included positive interaction effects involving Cau.R–GP.R, NAc.R–NAc.L, CCa.R–Ia.L, and CCp.R–PFCdm.L, as well as a negative interaction effect involving CCa.L–NAc.R. Seed-based connectivity analysis using NAc as the region of interest further confirmed this spatial agreement, particularly prominent for NAc.R (r = 0.720, p < 0.001; Figure 5a) with more modest effects seen for NAc.L (r = 0.503, p < 0.05).
Virtually predicted and empirical clinical outcomes of DBS in MDD
We quantitatively evaluated our model’s ability to reproduce DBS-induced network changes by computing the Euclidean distance between simulated and empirical DBS-ON networks in MDD. Network similarity analyses revealed significantly greater correspondence between simulated and empirical DBS-ON states compared to DBS-OFF baselines in MDD (Wilcoxon signed-rank test, p < 0.01). Importantly, these results were consistent across both OCD and MDD datasets, demonstrating robust performance in simulating DBS-driven network reorganization. To determine whether these connectivity changes were governed by underlying network topology, we contrasted our model with a topological null model (preserving edge counts, amplitudes, and node degrees). Significantly better fit was observed for the true biological network organization compared to randomized connectivity patterns (p < 0.001; Figure 5b), indicating that DBS effects depend on network architecture rather than chance configurations. Target specificity was further investigated by comparing NAc stimulation to alternative regions. Network distance analysis confirmed that NAc stimulation best recapitulated empirical effects (Figure 5c), supporting NAc as the optimal gray-matter proxy node within the present model.
To assess individual-level clinical outcome prediction in MDD, we applied the same two-step framework. The resulting multivariate FC signature comprised 20 connections, primarily involving DMN-related cortical connectivity together with subcortical and insular connections (Figure 5d). Under LOOCV, predicted and observed DBS-ON HAM-D scores were positively correlated (Pearson’s r = 0.659, 95% CI [0.097, 0.902], p = 0.028; R 2 = 0.433). The PLSR model fitted to the empirical DBS-ON data was then applied to the corresponding connections in the individually simulated FC. Model-estimated and observed HAM-D scores were also positively correlated (r = 0.627, 95% CI [0.043, 0.892], p = 0.039; R 2 = 0.071; Figure 5e). Additional permutation analyses provided further support for the stability of the empirical FC–symptom association relative to randomized clinical labels and randomly sampled edge sets from the same candidate feature space (Supplementary Figure S5b). Together, these findings suggest that the symptom-related multivariate FC pattern identified in the empirical DBS-ON data was also expressed in the simulated FC states.
Computationally guided prioritization of potential DBS candidates
We implemented our computational framework to make personalized response predictions for NAc-DBS across 92 OCD patients from an independent cohort. Patient-specific virtual stimulation effects were simulated from baseline FC patterns and quantitatively compared to empirical DBS-ON FC obtained from 11 DBS-treated patients. Through permutation-based ranking analysis (1,000 iterations), we identified three exceptional patients demonstrating significantly stronger predicted responses to NAc stimulation compared to the full cohort (p = 0.003; Figure 6b). Remarkably, the simulated connectivity patterns of these predicted responders showed strong spatial correspondence with empirical DBS effects observed in implanted patients. Regional analyses revealed particularly robust agreement in frontal (r = 0.82) and subcortical (r = 0.48) circuits (all p < 0.001, Bonferroni-corrected; Figure 6c). Computational target prioritization further confirmed the therapeutic relevance of NAc stimulation for these individuals, consistently ranking it among the top candidate targets (Figure 6d).
As a proof-of-concept candidate-prioritization analysis, the model-prioritized candidates showed elevated baseline severity, with Y-BOCS scores ranging from 35 to 40, and model-estimated reductions exceeding the conventional >35% Y-BOCS responder threshold used in prior OCD DBS studies (Gadot et al., Reference Gadot, Najera, Hirani, Anand, Storch, Goodman, Shofty and Sheth2022). Because these individuals had not undergone DBS, these estimated reductions should be interpreted as hypothetical model-derived outcomes rather than observed clinical improvements. Together, the convergent model-derived features, including network-state similarity to empirical DBS-ON FC patterns, NAc target specificity, and estimated symptom reduction, support the feasibility of using the framework as a hypothesis-generating tool for candidate prioritization in future NAc-DBS studies.
Discussion
Our findings demonstrate that DBS engages large-scale networks in a systematic, disease-specific manner. In OCD, NAc stimulation produced highly convergent alterations within the fronto–subcortical system, recapitulating canonical CSTC (cortico-striato-thalamo-cortical) models of compulsivity (Hollunder et al., Reference Hollunder, Ostrem, Sahin, Rajamani, Oxenford, Butenko, Neudorfer, Reinhardt, Zvarova, Polosan, Akram, Vissani, Zhang, Sun, Navratil, Reich, Volkmann, Yeh, Baldermann and Horn2024). These results suggest that effective intervention reduces pathological striatal gating while enhancing flexible frontal integration. Conversely, NAc-DBS in MDD primarily reorganized the DMN and associated affective circuits. This distribution aligns with contemporary models identifying DMN-centered maladaptive self-referential processing as a core depressive phenotype (Lynch et al., Reference Lynch, Elbau, Ng, Ayaz, Zhu, Wolk, Manfredi, Johnson, Chang, Chou, Summerville, Ho, Luecke, Bukhari, Buchanan, Victoria, Solomonov, Goldwaser, Moia and Liston2024; Tozzi et al., Reference Tozzi, Zhang, Pines, Olmsted, Zhai, Anene, Chesnut, Holt-Gosselin, Chang, Stetz, Ramirez, Hac, Korgaonkar, Wintermark, Gotlib, Ma and Williams2024). Furthermore, modulation within the salience system indicates a rebalancing of affect-biased attention, a recognized hallmark of MDD pathology (Segal et al., Reference Segal, Parkes, Aquino, Kia, Wolfers, Franke, Hoogman, Beckmann, Westlye, Andreassen, Zalesky, Harrison, Davey, Soriano-Mas, Cardoner, Tieg, Yücel, Braganza, Suo and Fornito2023). Collectively, these results suggest a unifying cross-diagnostic principle: DBS does not impose a generic “stimulation signature” but instead interacts with the intrinsic architecture of each disorder to redirect pathological dynamics toward adaptive configurations.
The present framework should be interpreted as a simplified, retrospective model of DBS-related FC-state differences rather than a complete model of therapeutic efficacy. The virtual perturbation captures how an abstract input propagates through a deconvolution-derived functional network under a first-order linear dynamical assumption, whereas the empirical DBS-ON/OFF data reflect resting-state FC differences between stimulation conditions in implanted patients. To evaluate whether the simulated responses were biologically meaningful, we did not rely solely on the resting-state FC used to construct the model; instead, we examined simulated-to-empirical DBS-ON network distances, comparisons with empirical DBS-OFF/DBS-ON differences, topological null models, target-sensitivity analyses, and cross-disorder consistency across OCD and MDD. These analyses support the model’s ability to capture aspects of DBS-related FC-state similarity, but they do not imply that a linear model fully captures brain dynamics. Chronic therapeutic response likely also depends on nonlinear stimulation effects, white-matter activation, neuroplasticity, medication status, and state-dependent clinical changes that are not captured here. Thus, the observed associations with symptom improvement should be considered exploratory, and true prospective prediction will require independent validation using pretreatment imaging and observed post-DBS outcomes.
While group-level analyses capture shared network principles, neuromodulation research is increasingly shifting toward individualized, circuit-informed modeling. The present framework offers a way to explore how patient-specific functional topology may shape simulated responses to candidate stimulation targets. Similar modeling principles have matured in neurological disorders – where larger implanted cohorts and richer longitudinal stimulation data provide important opportunities for prospective validation. In epilepsy, virtual brain twin models have been developed to estimate patient-specific epileptogenic networks and simulate stimulation responses, setting a strong precedent for individualized network-based perturbation modeling (Wang et al., Reference Wang, Triebkorn, Breyton, Dollomaja, Lemarechal, Petkoski, Sorrentino, Depannemaecker, Hashemi and Jirsa2024, Reference Wang, Dollomaja, Triebkorn, Duma, Williamson, Makhalova, Lemarechal, Bartolomei and Jirsa2025). In Parkinson’s disease, fMRI- and machine-learning-based approaches have been applied to identify clinically relevant DBS parameter settings, and recent adaptive DBS studies demonstrate how implanted sensing combined with symptom-state tracking can guide individualized stimulation strategies (Boutet et al., Reference Boutet, Madhavan, Elias, Joel, Gramer, Ranjan, Paramanandam, Xu, Germann, Loh, Kalia, Hodaie, Li, Prasad, Coblentz, Munhoz, Ashe, Kucharczyk, Fasano and Lozano2021; Oehrn et al., Reference Oehrn, Cernera, Hammer, Shcherbakova, Yao, Hahn, Wang, Ostrem, Little and Starr2024). These neurological applications thus offer valuable test beds for prospectively validating network-based virtual perturbation models in larger implanted cohorts. Within psychiatry, our framework extends this modeling logic to OCD and MDD by simulating DBS-related network responses across candidate targets in silico. Notably, the model prioritized the thalamus and caudate in OCD over the clinically conventional NAc. Rather than indicating clinical superiority of these targets, this divergence may reflect the functional continuum of ventral striatal and CSTC circuitry, where nearby subcortical regions can engage overlapping network gradients (Oldehinkel et al., Reference Oldehinkel, Tiego, Sabaroedin, Chopra, Francey, O’Donoghue, Cropley, Nelson, Graham, Baldwin, Yuen, Allot, Alvarez-Jimenez, Harrigan, Pantelis, Wood, McGorry, Bellgrove and Fornito2023; Webb et al., Reference Webb, Hearne, Tian, Zalesky, Robinson, Hall, Sonkusare, Burgher, Breakspear, Meyer, Horn, Naze, Mosley and Cocchi2025). Empirical fMRI from the same DBS cohort further supports this specificity: baseline functional variations in both regions predict individual therapeutic outcomes following vALIC DBS (Runia et al., Reference Runia, van de Mortel, Smith, Bergfeld, de Kwaasteniet, Luigjes, van Laarhoven, Notten, Beute, van den Munckhof, Schuurman, Denys and van Wingen2026). The thalamus’s prominence further aligns with its role as a high-leverage integrator within the CSTC loop (Hwang, Bertolero, Liu, & D’Esposito, Reference Hwang, Bertolero, Liu and D’Esposito2017). However, model-derived target sensitivity should not be conflated with clinical efficacy. Actual symptom relief depends on additional factors, including structural constraints, stimulation fields, and chronic neuroplasticity, that require future empirical validation.
Our framework employs network deconvolution to transcend the inherent constraints of conventional FC analysis. Network deconvolution estimates a first-order functional network by inverting the adjacency matrix’s transitive closure, thereby separating putative direct functional dependencies from higher-order propagation effects (Chen et al., Reference Chen, Zhang, Li, Huang, Lv, Yu, Chen, Sun and Wang2018; Feizi et al., Reference Feizi, Marbach, Médard and Kellis2013). The resulting sparse, interpretable matrices efficiently simulate DBS-related network responses across biologically relevant circuits (Lin et al., Reference Lin, Pan and Pan2024; Sipes, Arab, Nagarajan, & Raj, Reference Sipes, Arab, Nagarajan and Raj2026). Because this network is estimated from undirected resting-state FC, however, it represents model-based putative direct functional dependencies rather than causal effective connectivity. Consequently, the framework does not determine whether the underlying anatomical pathways are bidirectional or directional. Future studies incorporating perturbational recordings, electrophysiology, structural connectivity, or directed effective-connectivity models will be needed to establish its physiological and causal interpretation. While modern AI architectures (e.g. recurrent neural networks and transformers) achieve state-of-the-art performance in modeling spatiotemporal brain dynamics and estimating latent network interactions (Hollunder et al., Reference Hollunder, Ostrem, Sahin, Rajamani, Oxenford, Butenko, Neudorfer, Reinhardt, Zvarova, Polosan, Akram, Vissani, Zhang, Sun, Navratil, Reich, Volkmann, Yeh, Baldermann and Horn2024; Zhang et al., Reference Zhang, Phi, Li, Gorzek, Zwingenberger, Huang, Zhou, Kingsbury, Raam, Wu, Wei, Kao and Hong2025), their dependence on massive datasets renders them impracticable for rare clinical populations such as DBS cohorts. Our approach addresses this gap: it retains mechanistic interpretability while enabling personalized DBS outcome prediction, offering clear insights into the circuit-level alterations underlying OCD and MDD pathology. By marrying theoretical rigor with clinical utility, this framework delivers a transparent, principled alternative to opaque deep-learning models.
Despite these insights, several factors constrain the generalizability of our current results. First, the relatively small DBS cohorts limited statistical power and increased the uncertainty of the edge-wise analyses. Although the uncorrected group-by-treatment interaction effects showed anatomically structured and directionally similar patterns between the simulated and empirical data, no individual connection survived FDR correction. Accordingly, these edge-level findings should be regarded as exploratory and hypothesis-generating rather than as confirmatory evidence of DBS-induced changes in specific connections. Our primary interpretation is therefore restricted to the correspondence between distributed network-state patterns, rather than inference regarding individual DBS-sensitive connections. Future studies with larger, independent DBS cohorts will be required to determine the reproducibility of these edge-level effects. The modest size of our DBS cohort reduces statistical power for assessing medication effects or subtype-specific responses while also hindering validation of advanced AI approaches like LLMs or deep generative models in neuromodulation (Oliveira et al., Reference Oliveira, Coelho, Guedes, Cattoni, Carvalho and Duarte-Batista2024; Thirunavukarasu et al., Reference Thirunavukarasu, Ting, Elangovan, Gutierrez, Tan and Ting2023) – tools with substantial promise but underexplored risks. Moreover, the limited sample size precluded adequately powered analyses of potential clinical confounds – including medication status, illness duration, baseline symptom severity, and stimulation settings. While individualized FC patterns may partly reflect such clinical heterogeneity, this does not constitute formal adjustment for these factors. These variables may influence both resting-state FC and clinical outcomes, and their independent contributions cannot be estimated in the present study.
Although we applied rigorous preprocessing, fMRI signal quality was constrained by the 1.5T safety protocols; emerging 3T-compatible DBS hardware (Kutscha et al., Reference Kutscha, Mahmutovic, Bhusal, Vu, Chemlali, Hansen, May, Knake, Golestanirad and Keil2025) will enable higher-fidelity data in future studies. Moreover, the model frequency and amplitude parameters represent abstract perturbation terms in the discrete-time dynamical system, and the sinusoidal input should not be interpreted as the physical waveform of clinical DBS. The current parameter search therefore identifies best-fitting perturbation patterns that reproduce empirical DBS-ON FC states, while translation to device-level settings such as stimulation frequency, voltage, current, pulse width, waveform, or electric-field distribution remains unresolved. Bridging this gap requires integrating human fMRI with controlled stimulation protocols and cross-species biophysical calibration (Boutet et al., Reference Boutet, Madhavan, Elias, Joel, Gramer, Ranjan, Paramanandam, Xu, Germann, Loh, Kalia, Hodaie, Li, Prasad, Coblentz, Munhoz, Ashe, Kucharczyk, Fasano and Lozano2021; Hater et al., Reference Hater, Courson, Lu, Diaz-Pier and Manos2026). Finally, while we focused on the NAc, validating our framework across alternative targets (e.g. subthalamic nucleus or internal capsule) is critical to confirm its broad psychiatric utility. Additionally, our MDD simulations used the NAc as a gray-matter proxy for DBS delivered along the vALIC/NAc axis. This approximation was necessary because the present rs-fMRI framework models atlas-defined gray-matter nodes rather than white-matter tracts, stimulation volumes, or electric fields. Although our node-level sensitivity analyses compared NAc perturbation with alternative candidate regions, this analysis should be interpreted as a model-based assessment of gray-matter target sensitivity rather than an explicit reconstruction of vALIC fiber engagement. Future studies integrating electrode localization, electric field or volume-of-tissue-activated modeling, and diffusion tractography will be needed to characterize the distributed white-matter pathways engaged by vALIC DBS and to validate the framework across alternative psychiatric stimulation targets.
Conclusion
This study presents a computational framework for retrospectively modeling DBS-related network-state differences using individualized resting-state FC and abstract virtual perturbations. By linking simulated FC states to empirical DBS-ON/OFF network patterns and exploratory symptom associations, the framework provides a hypothesis-generating tool for studying how baseline functional architecture may shape DBS-related network responses and individualized sensitivity to candidate stimulation targets and model perturbation parameters. Ultimately, this framework provides a biologically grounded, data-driven tool for exploring the personalization of neuromodulation – supporting mechanism-guided therapy and the generation of testable hypotheses regarding quantifiable outcome prediction tailored to each patient’s unique brain connectivity profile.
Supplementary material
The supplementary material for this article can be found at http://doi.org/10.1017/S0033291726106461.
Author contributions
Z.W. and W.T.M. contributed to the conception and design of the study. J.H.L., S.Q.Y., X.Y.C., Q.B.J., H.Y.L., R.H.J., G.W., E.A.F., D.D., W.T.M., and Z.W. contributed to acquisition, post-processing, and analysis of the data. J.H.L., S.Q.Y., and X.Y.C. analyzed all or parts of the data, drafted the text, and prepared the figures with the help of Z.W. Z.W. wrote the manuscript, obtained the funding, and supervised the study. All authors approved the final version of the manuscript.
Funding statement
This study was supported by Brain Science and Brain-like Intelligence Technology – National Science and Technology Major Project (No. 2021ZD0204002), the National Natural Science Foundation of China (No. U25C2016, 82151303), and grants from Peking-Tsinghua Center for Life Sciences.
Competing interests
The authors declare that they have no competing interests.
Open access
