A deep representation learning model to predict response to vagus nerve stimulation

Domain: chemistry
Published: 2026-04-07
Editorial outcome: major_revision

Abstract
Implantable neurotechnologies are increasingly used to reduce seizure burden in epilepsy. Vagus nerve stimulation (VNS), the most common option, is effective for only half of patients, with no means to predict outcome prior to surgery. As a result, many children undergo invasive and costly procedures without benefit. Although T1-weighted magnetic resonance imaging (T1w) is routinely acquired presurgically and may capture structural brain differences relevant to treatment outcome, its high dimensionality relative to sample sizes has limited its utility in predictive modeling. To address this challenge, we present VQ-VNS, a deep representation learning model to predict VNS outcome based on preoperative T1w (n=263). First, we present data from the largest pediatric VNS cohort (n=1046), wherein presurgical clinical data could not predict response (AUC 0.54,p>0.99). Next, VQ-VNS was pretrained on 7433 T1w images to learn compact anatomical representations enabling its classifier to predict VNS response (AUC=0.73,p=0.007). Model predictions localized to serotonin-rich brain regions and inferred large-scale disruptions in network connectivity among non-responders. This biologically-interpretable predictor based on routine structural imaging improves upon current clinical decision-making.

Introduction
Epilepsy is the most common serious neurological disorder in children.1 The medical and psychosocial burden of drug-resistant epilepsy2 (DRE) has driven the increased usage of implantable neurotechnologies to improve seizure control. Vagus nerve stimulation (VNS) is the most commonly recommended device,3 supported by decades of clinical experience in pediatric populations. VNS has been shown to reduce seizure frequency,4 improve quality of life,5 and decrease mortality.2 However, only half of patients realize a benefit after treatment, and clinicians are currently unable to identify responders to VNS prior to surgery,6–8 leading to unnecessary, invasive, and costly procedures. Given that there is little variability in implantation technique,9 heterogenous outcomes are likely related to intrinsic differences among patients – particularly in brain structure and function. There is an unmet need to presurgically identify responders to VNS to personalize treatment options.

T1-weighted (T1w) MRI, a routine component of presurgical evaluation in pediatric epilepsy,10 represents a high-dimensional feature space ( R > 6 million voxels per image) that may contain complex, prognostically-relevant biomarkers. However, it remains underutilized in predictive modelling, as prognostic signatures are likely embedded in subtle, distributed, and nonlinear anatomical patterns11–13 – challenges compounded by the modest sample sizes typical of pediatric imaging datasets. Approaches to identify VNS responders based on coarse imaging features such as the presence of lesions14,15 or structural etiologies of epilepsy15–17 have not been informative and reported inconsistent associations. Region-based morphometry,18 is constrained by atlas-based parcellations that fail to capture diffuse, whole-brain patterns and yields weak associations with outcome. Crucially, there are no studies that have developed predictive models for VNS response on an individual patient level that are suitable for clinical use.8 Consequently, structural imaging has yet to be adopted as a viable predictor of VNS response.

Recent advances in generative modelling have enabled self-supervised learning of compact whole-brain embeddings from high-dimensional imaging data, capturing subtle and distributed anatomical patterns without reliance on predefined features. Architectures such as vector-quantized variational autoencoders19 (VQ-VAEs) and latent diffusion models20 have leveraged these lower dimensional representations to synthesize realistic brain images while modelling interindividual variability.21,22 While primarily developed for image synthesis, these models' principal utility in clinical contexts lies in their ability to learn compressed, semantically rich representations of whole-brain anatomy. Yet, their application to clinical prediction tasks remains underexplored.23 Moreover, most prior work has focused on adult datasets – such as the UK BioBank24 – with limited attention to pediatric populations despite their increased neurodevelopmental variability.25

In addition to technical barriers, clinical deployment of deep learning in neuroimaging is often limited by clinician apprehension about model opacity. Interpretability is increasingly viewed as a prerequisite for clinical adoption.26 Consequently, there is growing interest in addressing this requirement by either developing inherently interpretable models or applying post hoc analyses to explain the outputs of models not explicitly designed for interpretability.27 In the context of implantable neurotechnologies, where treatment decisions often involve surgical risk and uncertainty, predictive models must offer both accuracy and mechanistic transparency, grounded in neurobiological plausibility.

In this work, we present VQ-VNS, a clinically focused deep learning framework for predicting response to VNS directly from presurgical T1w MRI. Our approach combines a VQ-VAE, trained to learn discrete whole-brain anatomical representations, followed by a linear support vector machine28 (SVM) classifier for outcome prediction. To provide interpretability, we integrate gradient-weighted class activation mapping29 (Grad-CAM) and Shapley additive explanations30 (SHAP) to generate voxel-level predictive salience maps. To examine the neurobiological mechanisms underlying model performance, we also related the salience maps to normative neurochemical brain architectures as well as large-scale brain networks disrupted in non-responders using normative connectomic data. Our results demonstrate that interpretable deep representation learning can be leveraged to predict VNS outcomes while offering insight into the neurobiological mechanisms underlying treatment response.

Results
Clinical data do not predict response after VNS
We compiled the largest database of children to undergo VNS (n = 1046) from 11 North American centers through the Connectomics and Vagus Nerve Stimulation Outcome Study (CONNECTiVOS, Fig. 1a).28 Participants were deemed responders if seizure frequency was reduced by ≥80% and non-responders if ≤20% (n = 694, Fig. 1b) to reflect the bimodal distribution of seizure outcomes and mitigate the effects of ambiguous intermediate outcomes.29 The mean [SD; range] age was 10.6 [4.7; 0.8-24.2] years, 54% male and included 51% responders (Fig. 1b-c). This age range was included to capture the full developmental trajectory of brain maturation amongst pediatric and young adult populations that undergo VNS.

We first examined the association of each presurgical clinical variable with response by applying ordinary least squares regression (Fig. 1d; Supplementary Fig. 1). No clinical features were significantly associated with VNS outcome after false discovery rate correction (Supplementary Table 1). Nonlinear embedding of the clinical data using uniform manifold approximation and projection31 (UMAP) also failed to separate responders and non-responders (Fig. 1e; s = –2.7 × 10-4, p = 0.71), which was also corroborated by principal component analysis which revealed low explained variance (<12% for PC1) and no systematic response associated differences (Supplementary Fig. 2, Supplementary Table 2). Labelling the UMAP embeddings by treatment site revealed visually apparent but objectively weak (s = –0.04, p > 0.99) site-specific groupings (Fig. 1f), indicative of subtle inter-site variability in patient selection criteria without strong site-level effects.

Finally, consistent with findings in smaller cohorts,6,8 an SVM model trained on presurgical clinical features – including demographics, epilepsy history, and electroencephalography – was unable to distinguish responders to VNS from non-responders beyond chance-level performance (AUC = 0.54, p > 0.99; Fig. 1g). Furthermore, there was no difference in output current (1.9  0.7 vs. 1.7  0.7; t(101) = -1.4; p = 0.17; n = 103) or duty cycle (21  10% vs. 19  10%; t(101) = -1.3, p = 0.2, n = 103) between responders and non-responders.

Self-supervised representation learning captures anatomical patterns associated with VNS outcome
To study whether structural brain differences could predict outcome, we developed VQ-VNS, a deep representation learning model to identify VNS responders based on standard T1w sequences (Fig. 2). Because each T1w MRI scan contains over six million voxels – far exceeding the sample sizes of most pediatric epilepsy datasets (typically < 100 patients per center)32 – dimensionality reduction is critical to mitigate overfitting and enable effective model training. However, conventional approaches to dimensionality reduction, such as region-based parcellation, are insufficient to model the relationship between structural alterations and response following VNS (Supplementary Tables 3-4). Thus, the VQ-VNS model consists of two stages: 1) a VQ-VAE19 architecture compresses the high-dimensional MRI data into a compact latent representation, and (2) a support vector machine (SVM) classifies patients based on these latent features. The VQ-VAE architecture is well suited to compressing complex anatomical data into compact latent representations.21 We trained this model on large datasets of healthy subjects and persons with epilepsy who had not received VNS (n = 7,433; Fig. 2a–c). These datasets included participants across the lifespan (age 1.9 - 88.9 years) and contained a wide range of anatomical and acquisition-related variability to support the VQ-VAE’s ability to learn robust, generalizable representations. We optimized the VQ-VAE model with a composite loss term – comprising mean absolute error, perceptual loss,33 spectral loss34, commitment loss19 and an adversarial term derived from a discriminator35 model – to generate high-quality reconstructions of the input T1w image. This composite loss function was based on prior experiments21 where extensive ablations were performed to evaluate the contribution of each loss component.

The VQ-VAE’s reconstructions preserved anatomical details and captured biologically meaningful variation (Fig. 3a, Supplementary Fig. 3). We evaluated reconstruction fidelity using established metrics, including multiscale structural similarity36 (MS-SSIM), Fréchet Inception Distance37 (FID), and mean squared error (MSE). The VQ-VAE demonstrated strong performance (Supplementary Table 5), with a mean (SD) MS-SSIM of 0.97 (0.01), FID of 1.5 × 10-3 (1.2 × 10-3), and an MSE of 3.7 × 10-3 (8.6 × 10-4). To assess anatomical plausibility, we segmented both original and reconstructed images using SynthSeg38 and computed structure-wise correlations and Dice-Sørensen coefficients across four tissue classes: grey matter, white matter, subcortical grey matter, and ventricles. The VQ-VAE reconstructions demonstrated high correlation with true volumes across all four classes, with Pearson correlation coefficients ranging from 0.96 to 1.00 (Fig. 3b); furthermore, the reconstructions exhibited strong spatial overlap, yielding mean Dice-Sørensen coefficients ranging from 0.78 to 0.91 (Supplementary Table 6). Reconstruction fidelity was lower in subjects with significant structural abnormalities (Fig. 3a; Supplementary Fig. 4) with a mean (SD) SSIM of 0.95 (0.01), FID of 4.1 × 10-3 (9.8 × 10-4), and MSE of 4.7 × 10-3 (9.2 × 10-3), consistent with the VQ-VAE being trained primarily on anatomically typical brains.

Finally, we assessed whether this learned latent space captured features associated with VNS response without additional training. We encoded each subject’s image from the training subset of children who underwent VNS in the CONNECTiVOS cohort (199 of 263 images, Supplementary Fig. 5), computed the median latent representations separately for responders and non-responders, quantized these representations, and decoded them to generate synthetic archetypal brains (Fig. 3c). Compared to the archetypal responder, the resulting archetypal non-responder showed reduced grey matter, white matter and subcortical volumes, together with ventricular enlargement (Fig. 3d). These differences suggest that the latent representations captured structural features relevant to treatment outcome without direct optimization for VNS outcome prediction and that VNS non-responders exhibit lower brain volumes (i.e. atrophy) compared to responders.

Latent representations support accurate prediction and spatial interpretation of VNS outcome
To evaluate the predictive utility of the latent space for response following VNS, we trained a linear SVM classifier (Fig. 2e) using quantized latent vectors from the frozen VQ-VAE encoder. To reduce dimensionality, each four-dimensional latent tensor (x, y, and z spatial dimensions × feature channels) was spatially averaged across the x, y, and z axes, yielding a single 32-dimensional21 feature vector per subject. Each element of this vector therefore captures the global magnitude of the respective latent feature learned from the underlying T1w image. We then applied two-sided t-test-based feature selection on the training data (n = 199) to retain 27 of the 32 latent features (Supplementary Table 7), which were used to train the SVM, which demonstrated the strongest predictive performance among multiple classifiers (Supplementary Fig. 6). To assess the stability and significance of the learned features, we first trained 1000 SVM classifiers on randomly partitioned subsets of the VNS cohort. The mean AUC across models was significantly greater than that expected under the null distribution (mean AUC = 0.61; t(1998) = 28.6, p = 1e-150; Fig. 3e), suggesting that the latent representations carried measurable predictive signal.

We then trained a final SVM classifier using the training (n = 199) and validation (n = 25) subsets of the VNS cohort and tested on an independent held-out set (n = 39). The resulting model achieved an AUC of 0.73 in predicting response to VNS (permutation p = 0.007; Fig. 3f; Supplementary Fig. 7; Supplementary Table 8). The model maintained its predictive capacity across diverse age groups and treatment centers (Supplementary Fig. 7C). Model performance remained robust with the inclusion of subjects with ambiguous outcomes binarized at the 50% seizure reduction threshold (n = 53; AUC = 0.67, permutation p = 0.02; Supplementary Fig. 8). To evaluate the classifier’s ability to generalize across acquisition sites, we performed site hold-out validation, in which all images from one site were excluded during training and used exclusively for testing (Supplementary Table 9). Sites contributing small sample sizes were combined to create a single evaluation set. The classifiers achieved above chance accuracy (0.61 to 0.67), suggesting site-invariant structural features predictive of VNS response.

To evaluate the specificity of the SVM to predicting VNS outcomes, we tested its performance on an independent cohort of children who underwent temporal lobectomy for temporal lobe epilepsy39 (n = 102; mean age [SD] = 12.8 [3.9] years; 69.6% seizure freedom). The SVM model did not predict outcomes above chance level after temporal lobectomy (AUC = 0.57, permutation p = 0.12; Fig. 3f), suggesting that its decision boundary did not generalize to this population.

To visualize brain regions contributing most to model predictions, we implemented a hybrid salience attribution method integrating SHAP30 values from the SVM with Grad-CAM29 from the VQ-VAE encoder. Subject-specific SHAP values corresponding to the 27 selected latent features were expanded to a full 32-dimensional vector by assigning zeros to the unselected features. We computed the dot product between this SHAP vector and each subject’s latent representation, and backpropagated the resulting scalar through the VQ-VAE encoder. Gradients and activations were then extracted from the final convolutional layer to produce voxel-wise salience maps. This approach enabled SVM-derived feature importance to guide spatial salience attribution in the original image space. We identified several regions of high predictive value including the thalami, brainstem, and peri-Rolandic regions (Fig. 4a).

Model salience localizes to serotonergic networks disrupted in VNS non-responders
To better understand how model predictions relate to the chemoarchitecture of the brain, we correlated the VQ-VAE model’s salience map to the topology of brain receptor and transporter densities across nine different neurotransmitter systems from over 1200 healthy individuals.40 Neurotransmitter receptors are heterogeneously distributed across the neocortex and may be modulated by implanted technologies to modify neural states and reshape network-level communication.41 Each correlation was assessed against a null distribution generated via spin permutation testing, in which a single rotational affine was applied across all subjects’ latent salience maps for each permutation, prior to upscaling and averaging (Fig. 4b). We generated 1,000 such null maps per comparison to compute empirical p-values and account for spatial autocorrelation. The VQ-VAE’s attention was concentrated within brain regions enriched in serotonin transporter (5-HTT) expression (r = 0.26, p = 0.04; Fig. 4c).

The anatomically distributed pattern of salient regions also suggested involvement of broader brain networks in VNS response. To further interpret model predictions, we leveraged lesion network mapping, an approach which links focal structural abnormalities to disruptions in large-scale brain circuitry using normative connectomic data.42 This allowed us to test whether brain atrophy in non-responders, relative to 2,382 typically developing children (age 1.9 to 25.0 years, 63% male), converged upon a common network aligned with the VQ-VNS model’s saliency profile. We first demonstrated that cortical thickness measurements did not contain site-specific batch effects (Supplementary Fig. 9) and that group-level spatial patterns of atrophy alone did not clearly distinguish responders from non-responders (Supplementary Fig. 10). Regional brain atrophy maps converged upon large-scale structural (Fig. 4d, Supplementary Fig. 11) and functional (Fig. 4e) networks that were preferentially disrupted in non-responders to VNS. The VQ-VNS model’s saliency map showed significant spatial concordance (Supplementary Fig. 12a) with these structural (p<0.001, Supplementary Fig. 12b) and functional (p=0.04, Supplementary Fig. 12c) networks following spin permutation tests. Network findings were also concordant with brain regions to which vagus afferent fibres are known to project through the nucleus tractus solitarius,32 including areas that are activated by VNS.43

Methods
VNS Data Collection
Data were collected through a multi-institutional registry as part of the CONNECTiVOS study (Connectomic profiling and Vagus nerve stimulation Outcomes Study).59 This study involved children with epilepsy and was conducted across multiple North American sites, with local investigators engaged in design, data collection, and authorship. This study was approved by the institutional review boards of all participating sites (REB Approval Number 1000060282PID) and carried out in accordance with the Declaration of Helsinki.60 Data were collected from 1046 patients who underwent VNS implantation at centers in Canada and the United States. Detailed inclusion and exclusion criteria are described in the CONNECTiVOS protocol. All data were collected under local ethics approval with informed consent and assent where appropriate, and were de-identified prior to analysis. No material transfer, resource-poor settings, or elevated researcher risk were involved. We did not prespecify indications for VNS implantation, and all programming parameters were at the discretion of treating physicians.

For each participant, we collected continuous data on age at VNS implantation, duration of epilepsy, and percent seizure reduction at last follow-up based on parent-reported seizure frequency. Binary variables included: sex, presence of a genetic diagnosis, developmental delay, family history of epilepsy, history of infantile spasms, congenital or acquired epilepsy etiology, and prior epilepsy surgery. EEG data included variables indicating focal versus generalized seizure onset and lateralization of interictal epileptiform discharges (IEDs). Preoperative imaging data including High-resolution 3D T1-weighted (T1w) were collected both prospectively and retrospectively.

T1w MRI for the VNS cohort were acquired at either 1.5T or 3T. Scanner manufacturers include Philips, Siemens, General Electric, and Toshiba. Scanner models included Achieva, Skyra, Trio TIM, Ingenuity, Signa HDxt, Discovery MR750, Ingenia, Intera, Symphony, Prisma, Avanto, Signa Excite, Symphony TIM, Trio, Titan, Discovery MR750w, and Genesis Signa. Key parameters for T1-weighted imaging were: flip angle, 12.3° (8.0° – 40.0°); echo time, 3.7 ms (1.7 – 9.2 ms); repetition time, 7.9 ms (1.1 – 25.0 ms); and mean voxel size, 0.7x0.7x0.7mm.

Pretraining Data Compilation
We compiled data from a cohort of healthy subjects and persons with epilepsy from the Autism Brain Imaging Data Exchange (ABIDE, n = 412),61 the Healthy Brain Network (HBN, n = 2081),62 the Calgary Preschool MRI (n = 265),63 the Cambridge Centre for Ageing and Neuroscience (Cam-CAN, n = 649),64 and the Amsterdam Open MRI Collection (AOMIC, n = 923),65 SickKids internal normative dataset (n = 1260), The Imaging Database for Epilepsy And Surgery (IDEAS, n = 542),66 Focal Cortical Dysplasia Type II (FCD-2, n = 170),67 Temporal Lobe Epilepsy – UNAM (UNAM-TLE, n = 66), Human Connectome Project Young Adult (HCP, n = 1065)68 datasets for a combined dataset of 7433 images. To balance retaining a sufficient training set size with sufficient subjects for validation and robust assessment on a held-out test cohort, these data were split into training, validation and test sets in ratios of 75%, 10%, and 15% respectively, which is commonly used in the context of clinical datasets

The SickKids internal normative T1w images were acquired at either 1.5T or 3T across various scanners (Signa HDxt, GE Healthcare; Achieva, Philips Healthcare; Magnetom Skyra, Siemens Healthineers) with a dedicated head-coil. Images were acquired either in the sagittal or the axial plane: TR/TE, 6-9/2-5 ms; flip angle: 10-15 degrees, mean voxel size, 1.0x1.0x1.5 mm.

All T1w images were brain extracted,38,69 70resampled to 1mm3 isotropic resolution,71, bias field corrected, intensity normalized with Fuzzy C-means,72 and rigidly registered73 to a 1mm MNI brain extracted template (Fig. 2a).74

Analysis of clinical data
Participants were classified as non-responders if they had ≤20% seizure frequency reduction, and as responders for ≥80% seizure frequency reduction based on the bimodal distribution of treatment outcomes (Fig. 1b) and to mitigate misclassification risk inherent to seizure reporting from participants with ambiguous outcomes.75 Based on the outcome threshold, 694 of the 1046 subjects were included. We used ordinary least squares regression to model the univariable relationship between each clinical variable and VNS response. Counts of missing values for each clinical variable are presented in Supplementary Table 10. False discovery rate (FDR) correction was applied to adjust for multiple comparisons.

To model multivariable and non-linear relationships, we applied Uniform Manifold Approximation and Projection (UMAP).31 Silhouette scores76 quantified the clustering quality of UMAP embeddings by outcome and site, with significance assessed via one-sided one-sample t-tests. To further assess whether the observed manifold structure reflected biologically or methodologically meaningful variance, we additionally performed principal component analysis (PCA) on the same dataset. Loadings of the top principal components were examined to identify the metrics contributing most to overall variance within the clinical data. Scores for the leading components were then compared between responders and non-responders to evaluate whether these reduced dimensions captured differences associated with VNS outcome.

To assess the predictive value of clinical data, a support vector machine (SVM)36 was trained on 80% of the dataset (n = 550) to predict VNS outcome. The remaining 20% (n = 144) was held out as an independent test set to evaluate model performance after training. To ensure that both training and test sets contained only clearly defined responders and non-responders, subjects with ambiguous (intermediate) outcomes were excluded prior to this split. Recursive feature elimination with cross-validation (RFECV)77 was employed for feature selection. The clinical SVM model’s performance was quantified using area under the curve (AUC) and confusion matrices. We additionally evaluated model performance with the inclusion of subjects with ambiguous outcomes (n = 352) and observed no improvement.

Cortical Thickness Estimation and Regional-based Modelling of Response
Cortical thickness was estimated from T1-weighted images using FreeSurfer71 to derive individualized thickness maps used in all downstream atrophy and network-based analyses. Cortical surface reconstruction in FreeSurfer v7.4.2 began with an affine transformation aligning images to the MNI305 atlas, followed by B1 bias field estimation from white matter intensity variations and subsequent intensity inhomogeneity correction. Intensity normalization was then performed to correct for cross-site and cross-scanner fluctuations that would otherwise disrupt subsequent intensity-based segmentation steps. The skull was stripped, and white matter boundaries were automatically identified based on intensity and neighbour constraints. The white matter surface then underwent smoothing through tessellation and was expanded outwards to identify the boundaries of the pial surface. Cortical thickness was calculated as the average distance from a point on the smoothed white matter surface to its corresponding point on the pial surface and vice versa. The surfaces and parametric maps were spherically inflated and registered to a standard spherical atlas.78 This was followed by region-wise regression to VNS outcome with a generalized linear model accounting for age and estimated total intracranial volume (eTIV). The parcellations were derived from the Desikan-Killiany Atlas. Although this pipeline was established in adults, it has been validated for use in children,79 and all surface reconstructions were visually examined for accuracy.

Vector-Quantized Variational Autoencoder
T1-weighted MRI scans are inherently high-dimensional, making it challenging to train predictive models directly on raw voxel intensities – particularly in the context of relatively small clinical imaging datasets. This necessitates an approach to represent these data in a low dimensional space without loss of detail necessary for downstream classification tasks. We implemented a vector-quantized variational autoencoder (VQ-VAE)19 from the MONAI library to generate these lower dimensional latent representations. Compared to classical convolutional classifiers that require pretraining on surrogate labels such as age or diagnosis, using an autoencoder allowed the model to learn representations based solely on structural anatomy, without bias toward any specific pretraining objective.

During training, we employed a composite loss function incorporating voxel-wise accuracy, spectral alignment, perceptual similarity, encoder-quantizer alignment (commitment), and adversarial objectives (Summarized in Supplementary Table 11).21 Voxel-wise error (ℒ ) was quantified with mean absolute error (MAE) between the original 𝑋 and reconstructed 𝑋̂ images
ℒ = MAE(𝑋,𝑋̂) (1)
𝑣𝑜𝑥
Spectral alignment was enforced using a frequency loss function (ℒ ) computed as the mean squared error (MSE) between the fast Fourier transforms34 of the original and reconstructed images:
ℒ = MSE(FFT(𝑋),FFT(𝑋̂)) (2)
𝑓𝑟𝑒𝑞
Perceptual similarity was assessed using MONAI's perceptual loss function (ℒ ), which leveraged the Learned Perceptual Image Patch Similarity (LPIPS)33 metric with AlexNet80 as the feature extractor. Given that the 2D AlexNet model cannot directly process 3D volumes, we computed the LPIPS metric on individual slices extracted from the axial (ax), coronal (cor), and sagittal (sag) views, subsequently summing these per-slice scores as in:
ℒ = LPIPS(𝑋 ,𝑋̂ )+ LPIPS(𝑋 ,𝑋̂ )+ LPIPS(𝑋 ,𝑋̂ ) (3)
𝑝𝑒𝑟𝑐𝑒𝑝 𝑎𝑥 𝑎𝑥 𝑐𝑜𝑟 𝑐𝑜𝑟 𝑠𝑎𝑔 𝑠𝑎𝑔
Commitment loss (ℒ ) in the VQ-VAE encouraged the encoder outputs to commit to discrete codebook entries. Since the codebook was updated using an exponential moving average (EMA) approach,16 gradients should not flow directly into the codebook vectors. Accordingly, we applied a stop-gradient operator to the codebook entry:
ℒ = 𝛽‖𝑧 (𝑋)−𝑠𝑔(𝑒)‖ (4)
𝑐𝑜𝑚𝑚𝑖𝑡 𝑒 2
Where 𝛽 is a scaling constant, 𝑧 (𝑋) is the encoder’s output, 𝑒 is the closest codebook vector, and 𝑠𝑔(∙) denotes the stop gradient operator.

Adversarial loss (ℒ ) guided the generator to produce outputs that are indistinguishable from real images by a patch-based discriminator 𝐷,35 thereby retaining local consistency and fine structural detail. The loss employed a least-squares criterion with mean reduction, targeting discriminator outputs of 1 for generated samples. No activation was applied to the discriminator outputs prior to loss computation:
ℒ = (𝐷(𝑋̂)−1) 2 (5)
𝑎𝑑𝑣
Together the expression representing the total loss term was,
ℒ = ℒ +ℒ +ℒ +ℒ +ℒ (6)
𝑡𝑜𝑡𝑎𝑙 𝑣𝑜𝑥 𝑓𝑟𝑒𝑞 𝑝𝑒𝑟𝑐𝑒𝑝 𝑐𝑜𝑚𝑚𝑖𝑡 𝑎𝑑𝑣
Training stability across loss terms with differing magnitudes was maintained through the application of gradient normalization, a method inspired by GradNorm.81 This involved rescaling gradients according to their norms, as detailed below:
𝑔̃ = (7)
||𝑔||+ 𝜀
where 𝑔̃ represents the normalized gradient, 𝑔 is the original gradient, 𝑤 is a scaling factor, and 𝜀 is a small constant to prevent numerical instability. This approach rescaled all gradients based on their norms, effectively balancing their contribution to the training update and preventing dominance by large gradients or attenuation of smaller ones. The scaling factors were empirically selected.

The discriminator was trained using mean squared error (MSE) loss between its predicted scores for real and generated images, denoted as ℒ :
𝑑𝑖𝑠
1
ℒ = [MSE(𝐷(𝑋))+MSE(𝐷(𝑋̂))] (8)
𝑑𝑖𝑠 2
In addition to gradient normalization, we applied LeCam regularization82 to further stabilize adversarial training. This method limits how much the discriminator’s outputs can fluctuate, helping prevent mode collapse and improving generalization.

Let 𝐷(𝑋) and 𝐷(𝑋̂) denote the discriminator logits for real and generated images, respectively. We defined LeCam anchors 𝐴 and 𝐴 as exponential moving averages of real and fake logits:
𝑟 𝑓
𝐴 ← 𝛽𝐴 +(1−𝛽)𝐸[𝐷(𝑋)], 𝐴 ← 𝛽𝐴 +(1−𝛽)𝐸[𝐷(𝑋̂)] (9)
𝑟 𝑟 𝑓 𝑓
where 𝛽 is a momentum parameter and E represents the average of 𝐷(𝑋̂). The LeCam loss (ℒ ) penalized deviations of real logits from the fake anchor and fake logits from the real anchor:
ℒ = 𝐸[( 𝐷(𝑋)−𝐴 ) 2 ]+𝐸[( 𝐷(𝑋̂)−𝐴 ) 2 ] (10)
𝐿𝑒𝐶𝑎𝑚 𝑓 𝑟
This regularization discourages abrupt changes in discriminator outputs over training and prevents overconfident separation of real and generated inputs. The final discriminator loss combined the adversarial term with the LeCam penalty, scaled by a tunable coefficient:
ℒ = ℒ +𝜆 ∙ℒ (11)
𝑎𝑑𝑣 𝑑𝑖𝑠 𝐿𝑒𝐶𝑎𝑚 𝐿𝑒𝐶𝑎𝑚
where ℒ represents the discriminator loss at a particular instance and 𝜆 represents the LeCam anchor weight.

Vector-Quantized Variational Autoencoder Implementation
The input to the VQ-VAE model was a 3D volume of size 176 × 208 × 176. The encoder comprised four convolutional downsampling blocks with 128, 128, 256, and 256 feature maps, respectively, yielding a latent representation of size 11 × 13 × 11 with 32 channels. Quantization was performed using a codebook with 2,048 entries and embedding dimension 32. The resulting quantized latent had the same shape as the encoded latent. The decoder consisted of four transposed convolutional layers to reconstruct the final image.

The discriminator was a PatchDiscriminator from MONAI, consisting of four convolutional layers with 64 channels each, followed by a classifier head with a single output channel. Training was performed for 140 epochs using the Adam optimizer83 with a learning rate of 1e-5. The generator used 𝛽 = (0.9, 0.999), and the discriminator used 𝛽 = (0.5, 0.99). Gradient loss terms were weighted by empirically selected scaling factors: 0.1, 0.5, 1.0, 1.0, and 1.0 for ℒ ,ℒ ,ℒ ,ℒ and ℒ respectively. Training was conducted on a single node in the SickKids High Performance Computing cluster with an Nvidia H100 NVL GPU (96 GB VRAM), 18 CPU cores and 48 GB of RAM.

Reconstruction Quality Evaluation
To evaluate the fidelity of VQ-VAE reconstructions, we computed three complementary image-level metrics – mean squared error (MSE), multiscale structural similarity (MS-SSIM),36 and Fréchet Inception Distance (FID)37 between the original and the reconstructed 3D T1-weighted images. MSE and MS-SSIM were computed using MONAI (MSEMetric, MultiScaleSSIMMetric), and FID was calculated using features extracted from a pretrained ResNet-50 model,84 following global average pooling across spatial dimensions. To assess the VQ-VAE model’s robustness to anatomical abnormalities, we performed a secondary evaluation on 106 subjects with abnormal imaging findings. Reconstruction fidelity was evaluated using the same three metrics – MSE, MS-SSIM, and FID.

To assess anatomical plausibility, all reconstructed and original images were segmented using SynthSeg.38 Structure-wise agreement was evaluated across four tissue classes – grey matter, white matter, subcortical grey matter, and ventricles – by computing the Pearson correlation between reconstructed and original volumes across subjects.

Training the SVM Classifier to Predict VNS Response
Response classification used the same thresholds as applied to the clinical data: ≥80% seizure reduction for responders and ≤20% for non-responders. Subjects with significant anomalies, prior neurosurgical intervention (Supplementary Fig. 4a-d), or conflicting outcome assessments, were excluded. After filtering, the final dataset comprised 263 subjects (Supplementary Fig. 5). The 263 patients included in model training did not differ significantly from the excluded cohort in terms of demographic or clinical characteristics (Supplementary Table 12).

An SVM classifier was trained on vector-quantized variational autoencoder (VQ-VAE) latent embeddings to predict binary treatment outcomes. The VQ-VAE first encoded T1-weighted MRI scans into a compressed latent representation, which was averaged across its spatial dimensions and flattened to yield a one-dimensional 32 feature vector per subject.

Latent features were evaluated via two-sided t-tests, and only those showing significant group-level differences (p < 0.05) were retained (27 of 32, Supplementary Table 7). To assess the stability and significance of the learned features, we first trained 1000 SVM classifiers with class balancing and a high regularization parameter (C=5000) on randomly partitioned subsets of the VNS cohort. To assess generalizability, a final SVM was trained with a training, validation set, and held-out test set. To minimize data imbalance, participants were split into training (n = 199, responders = 100, non-responders = 99), validation (n = 25, responders = 13, non-responders = 12), and held-out test sets (n = 39, responders = 19, non-responders = 20) using stratified sampling to preserve approximately equal proportions of responders and non-responders in each subset. Model performance was evaluated using area under the curve (AUC) and confusion matrices across these sets. Robustness was assessed by permutation testing with 1,000 label shuffles to generate a null distribution of AUC scores, enabling empirical significance estimation. This procedure tests whether the VQ-VNS model’s observed performance exceeds chance expectations, thereby confirming that predictions reflect a learned mapping between input features and outcome.

To evaluate the SVM classifier’s ability to generalize across acquisition sites, we performed site-hold-out validation, in which all images from one site were excluded during training and used exclusively for test. To enable evaluation of smaller test sizes, sites were grouped by geographic region. To further assess robustness, we tested the VQ-VNS model’s performance in subjects with intermediate outcomes - seizure reductions between 20% and 80% - without further model training. For this analysis response was binarized at the conventional >50% threshold. Similarly, model performance was evaluated using AUC and confusion matrices. The VQ-VAE encoder was frozen during SVM training. Latent embeddings were generated using an Nvidia RTX A6000 (Ampere, 48 GB VRAM) GPU, after which SVM training and inference were performed on CPU. In inference mode – including VQ-VAE latent generation and SVM prediction – the average runtime per subject was 0.19 seconds.

To evaluate the specificity of the VQ-VNS model’s performance, we tested its ability to predict seizure freedom in a separate cohort of patients with temporal lobe epilepsy39 (TLE) who underwent temporal lobectomy (n = 102; mean age [SD] = 12.8 [3.9] years; 69.6% seizure freedom). This retrospective TLE cohort was selected because high-quality preoperative T1-weighted imaging and detailed postoperative outcome data were available, enabling an assessment of whether the VQ-VNS model, training exclusively on VNS recipients, generalizes to an independent surgical population. In this cohort, approximately 70% of patients achieved seizure freedom following surgery, which enabled the evaluation of model specificity. T1-weighted images were acquired using Philips and GE scanners (Achieva, Genesis Signa, Signa HDxt, Signa Excite) with field strengths of 1.5T or 3T. Imaging parameters were: flip angle = 13.6° (range: 8.0°–30.0°), echo time = 4.3 ms (2.3–4.6 ms), repetition time = 8.7 ms (5.0–25.0 ms), and mean voxel size = 0.9 × 0.9 × 1.4 mm. Images with gross abnormalities were excluded as was performed in the VNS cohort.

Model Interpretation Using Gradient-based Salience Mapping
To identify brain regions relevant for outcome prediction, we generated voxel-wise heatmaps using a hybrid approach inspired by both Shapley additive explanations85 (SHAP) and gradient-weighted class activation mapping29 (Grad-CAM). The VQVAE encoded each subject’s image into its quantized latent representation. The feature vector for each subject was passed to the trained SVM to predict response. SHAP values were then computed using a model-agnostic explainer with the median training feature vector as the reference. For each subject, the resulting SHAP vector was used in a dot product with their input feature vector. This scalar output was backpropagated through the encoder, and gradients were extracted from the final convolutional layer to generate latent salience maps based on gradient-weighted activations.

The latent salience maps were upsampled to the original image resolution and then spatially normalized to MNI space by applying the subject-specific nonlinear transformation from model T1 space to MNI space, previously estimated using ANTs.86 Individual maps were averaged across all subjects, and the resulting group-level salience map was converted to absolute values to reflect feature importance. The signed average salience map was retained for receptor and transporter expression analyses, which are sensitive to the direction of feature contributions. Visualizations of the salience map are log transformed to aid in identification of the most salient regions.

Mapping neurotransmitter receptor expression within salience map
To better understand the biological relevance of the regions in the VQ-VAE model’s salience map, we studied the associations to the normative chemoarchitecture of the brain using PET-derived neurotransmitter receptor and transporter expression profiles across nine different neurotransmitter systems, as compiled from previously published datasets.87 These datasets were acquired from over 1200 healthy individuals. Spearman correlations were computed between the signed average salience map and each receptor expression map.

Each correlation was permutation-corrected using a latent-space spin procedure analogous to cortical surface map rotations.88 Subject-level latent salience maps were rotated with random affine transformations prior to upsampling and spatial normalization. These upsampled maps were averaged to generate null group-level salience maps, which were correlated with each receptor expression map. Empirical p-values were defined as the proportion of null correlations exceeding the observed correlation. To account for multiple comparisons across all expression profiles, p-values were Bonferroni-corrected. Correlations with corrected p < 0.05 were reported.

Cortical Thickness based Lesion Network Mapping
To study whether model predictions align with large-scale brain network disruptions related to brain atrophy in non-responders to VNS, we employed a lesion network mapping (LNM) approach. Areas of cortical atrophy (relative to healthy age- and sex-matched subjects) were considered the “lesions”.

To evaluate potential site-related batch effects in the cortical thickness data, we applied t-distributed stochastic neighbor embedding89 (t-SNE) to all cortical vertex values and estimated total intracranial volumes (eTIV), projecting them into two dimensions. To preserve local neighbourhood structure in high-dimensional data, t-SNE was selected. Although t-SNE does not preserve global distances well, it offers an established and effective visualization tool for our goal of exploring relative clustering and assessing potential batch effects. Embeddings were visualized by site (Supplementary Fig. 9a), outcome, and age (Supplementary Fig. 9b–c). Age showed the strongest correlation with t-SNE component 1 (Supplementary Fig. 9d). As no consistent clustering by site was observed, no further harmonization was applied.

Mapping regions of brain atrophy to a common network with normative structural connectivity
We sought to delineate whether spatially distributed atrophic regions in individuals map onto a common brain network associated with VNS response. To this end, we employed LNM.90 We generated parametric maps to quantify each subject’s vertex-wise cortical thickness deviations from a normative population of 2382 children (Hospital for Sick Children, n=393, ABIDE,62 n=338, HBN,62 n=1371, Calgary Preschool,63 n=251, Cam-CAN,64 n=29, age 1.9 to 25.0 years, 63% male). We fit first-order linear models for age and eTIV for each sex at each vertex in the normative dataset and applied them to the VNS cohort. Cohort residuals were then normalized to the normative residual distribution to derive age-, eTIV-, and sex-adjusted z-scores, which we termed 𝑣-scores. Each subject’s 𝑣-map (full width at half maximum (FWHM) = 3mm) was thresholded at 𝑣 ≤ −1.64 and binarized to generate the “lesion” map for LNM. Subjects were excluded if the atrophic vertex count exceeded 5% (Supplementary Fig. 10a) of the total brain surface (n = 5) to minimize the confounding effects of global atrophy. After applying this threshold and FreeSurfer segmentation quality-based exclusion criterion, 115 of 263 subjects were retained for lesion network mapping analyses. The remaining subjects were excluded either for failed FreeSurfer reconstruction or incomplete sex information.

To perform lesion network mapping (LNM), we used a subset of the Human Connectome Project (HCP) diffusion MRI dataset68 (n = 218, age < 25). Diffusion data were preprocessed using TractoFlow, 91 which performs denoising, motion and eddy current correction, B0 normalization, and fibre response function (FRF) estimation. Tractography was performed using Particle Filtering Tractography (PFT)92 with interface-based seeding (60 seeds/voxel), a minimum streamline length of 5 mm, maximum of 200 mm, and spherical harmonics order 8. Only the b=1000 s/mm² shell was used. Mean FRFs were disabled to improve multi-site harmonization.

Cortical surfaces were generated using FreeSurfer v6.0 as part of the freesurfer_flow pipeline. Surface-Enhanced Tractography93 (SET) was then applied using probabilistic fODFs from TractoFlow, with an anatomically informed stopping criterion at the grey–white junction to minimize gyral bias. Streamlines were further filtered using convex optimization modeling for microstructure informed tractography94 (COMMIT) to retain only those contributing meaningfully to the diffusion signal.

Processing was implemented via Nextflow95 workflows and Singularity96 containers, using publicly available images from the Sherbrooke Connectivity Imaging Lab (https://scil.usherbrooke.ca/).

Each VNS subject’s binarized atrophy map was used to filter HCP tractograms, retaining only streamlines intersecting the atrophic cortex. The endpoints of these streamlines were projected to MNI space (1 mm³ resolution).97 Atrophy networks were defined by combining endpoint maps across all HCP subjects and retaining voxels present in ≥20% of subjects. This thresholding approach was conceptually adapted from the disconnectome framework,98 but tailored to the endpoint-based tractography analyses used in our study. Voxel-wise t-tests were performed comparing responders and non-responders across all non-zero voxels, with statistical significance determined via 3DClustSim99 (voxel-wise p < 0.05, cluster-level α = 0.05).

Mapping regions of brain atrophy to a common network with normative functional connectivity
To assess functional connectivity patterns associated with regional atrophy, we used the same atrophy maps described above as seed regions in a normative resting-state fMRI connectome. Preprocessed fMRI data were obtained from the Brain Genomics Superstruct Project (n = 1000).100 Preprocessing was conducted by the data providers and included removal of the first 4 volumes, slice-time and motion correction (MCFLIRT), exclusion of frames with framewise displacement > 0.2 mm or DVARS > 50, co-registration to anatomical images, regression of nuisance signals (motion, global signal, CSF, and white matter), bandpass filtering (0.0009–0.08 Hz), and projection to MNI-152 2 mm space.

For each VNS subject, the average BOLD timeseries from the subject’s atrophy map was extracted and used as a seed in voxel-wise linear regression across the normative dataset using AFNI’s 3dDeconvolve.99 Correlation coefficients were Fisher-transformed to Z-scores, and voxelwise t-statistics were computed, yielding one statistical map per subject (n = 115). These maps were compared between responders and non-responders using a voxel-wise generalized linear model implemented in FSL’s randomise,101 controlling for total atrophy volume. The resulting group-level statistical map was used in subsequent analyses assessing spatial concordance with the VQ-VAE model salience map.

Assessment of Network Concordance with Salience Map
To evaluate the spatial concordance between predictive salience maps and normative brain networks, we performed a voxel-wise overlap analysis using response-associated functional and structural connectivity networks derived from resting-state fMRI and diffusion tractography, respectively. All network maps were aligned to MNI space and resampled to 1 mm isotropic resolution.

We computed the mean of the element-wise product between the salience map and the structural and functional network maps to quantify the empirical spatial concordance. To test whether the spatial concordance was statistically significant, we leveraged the aforementioned latent-space spin permutation procedure, adapted from methods for rotating cortical surface maps to maintain spatial autocorrelation76 (Supplementary Fig. 12a). Subject-level latent salience maps were randomly rotated using affine transformations (consistent affine across all subjects per permutation), upsampled to image space via trilinear interpolation, and registered to MNI space using each subject’s pre-computed nonlinear warp. The resulting rotated maps were averaged across subjects to generate null group-level salience maps (n = 1000). For each, alignment was quantified using the same structural and functional procedures applied to the observed data, yielding a null distribution of concordance values. Statistical significance was defined as p < 0.05 based on the proportion of null values exceeding the observed concordance, with the 95th percentile visualized for context.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Aims:\nIn this manuscript, the authors seek to produce a predictive model of VNS response using T1w MRI. Specifically, this model\nprovides reasonable predictive power on paediatric cases binarising response between <=20 and >=80 seizure reduction,\nalthough this effect is reduced when using more traditional outcomes (<50 and >=50 seizure reduction, respectively). They\nidentify areas key to their model’s prediction of response (thalamic, brainstem, insula, and peri-Rolandic) notably these\ncorrelate well with regions known to be affected by VNS (VagAN). This work paves the way for the use of T1w MRI in\npresurgical evaluation for VNS in paediatric cohorts.\nThis study has many positives, and I congratulate the authors for this impressive work. The study is novel in both its\nmethods, data, and hypotheses. The sample size is huge compared to previous work, and the growing use of VNS makes it\ntimely. I consider the field of VNS research in epilepsy to be extremely worthwhile, and hope my comments can help\nimprove the authors improve this important study.\nMajor Points:\n1. The methodology in this paper is interesting as an application of modern machine learning tools. However, it would also\nbe good to know if simpler statistical analyses of basic measures are similarly informative (e.g. volume/surface\narea/thickness measures from freesurfer). Looking at figure 3c, perhaps thalamic volume alone is already somewhat\ninformative?\n2. If I understood correctly, the 39 patients in the withheld test data was chosen randomly. This withheld data of 39 patients\nrepresents only one of the possible sets of 39 patients that could be held-out. It would give more confidence in the approach\nif this could be permuted many times and error bars plotted on figure 3b, along with confidence intervals for AUCs.\n3a Following on from point 2 regarding test data. I think it would be worthwhile to test generalisability across data from\ndifferent hospitals. I.e. hold out data from one hospital and train on the others, then report test accuracy on the withheld\nhospital. This would be useful to know since it will tell us if the model can be used/transferred to other hospitals, which\nwould be clinically important.\n3b. Continuing from point 3a, are there any hospitals with data to which the model/approach does not generalise to very\nwell? If so why? (perhaps they have older/younger demographics? Poorer scan quality?)\n4. The authors focus their evaluation mainly on AUC. This is important, and I have no problem with this. I think it might also\nbe worthwhile reporting if the specificity at 100% sensitivity and vice versa. I think this would be worthwhile because this\nmay be useful for either ruling out definite non-responders, or perhaps as a screening tool for identifying/prioritising likely\nresponders.\n5. The regions identified as having a “Higher” feature value for prediction using the gradient-based saliency mapping are not\nprovided with a quantitative backing. Providing a numerical underpinning to justify the chosen regions would provide greater\nconfidence and understanding of whether these regions are uniquely involved or if the whole brain is still required. N.b. this\ncould potentially be partly addressed through point 1 above.\nMinor Points:\nReference 8 and reference 32 are the same (I have not checked, but it is possible this occurs more than once).\nIt is surprising given prior literature that you find so many individuals with >=80% and <=20% seizure reduction. Historically\nthese rates appeared far lower. Do you have any knowledge as to why this combined North American cohort appears to\nhave significantly more individuals at both extremities of response?\nIt is not identified how seizure reduction was calculated. Assumedly, this is through seizure diaries or discussed at follow-up\nappointments which should be clarified.\nWithin the methods you refer to training cohorts both as training and discovery cohorts consistent naming would be clearer.\n“Subjects with significant anomalies, prior neurosurgical intervention (Supplementary Fig. 3a-d), or conflicting outcome\nassessments, were excluded. After filtering, the final dataset comprised 263 subjects”. The exact number dropped for each of\nthese reasons would be interesting. A PRISMA flow diagram would be helpful as a supplementary figure.\nfigure 4d; top row, middle panel, the Insula arrow terminates in the putamen at the moment.\n“ Each four-dimensional latent tensor (spatial dimensions × feature channels) was spatially averaged across the x, y, and z\naxes, yielding a single 32-dimensional feature vector per subject. “ I struggled to follow/understand the rationale for the\naveraging across the spatial axes xyz. Please can you explain this in more detail? Perhaps in an expanded supplementary\nmethods.\nfigure 2g -I did not follow what is relative volume. relative to what? The responder blue bars on the chart are all at 100%, so I\nassume it’s relative to the responder. Perhaps it would be helpful to plot relative to healthy control and show a synthetic\ncontrol in figure 2f for reference\nI thought the validation in the TLE cohort was a brilliant idea. I was wondering what the interpretation would be if the model\nhad generalised to the TLE cohort. Do you think that would indicate shared mechanisms of response to VNS and also\nresponse to TLE surgery? As such, should we conclude that given the model did not generalise to TLE the mechanisms\ndiffer?", "reviewer_2": "Recommendation:\nMajor Revisions\nTitle:\nA Deep Representation Learning Model to Predict Response to Vagus Nerve Stimulation\nOverview and general recommendation:\nThis work investigated the use of T1-weighted MRI data for predicting clinical response to implantable vagus nerve\nstimulation (VNS) therapy in individuals with drug-resistant epilepsy. The authors leveraged their access to a network of\nclinical neuroimaging databases to aggregate a large cohort of preoperative and postoperative (resections etc.) T1w-MRI\nscans, which were then used for modeling development. A vector quantized variational autoencoder was used to model the\nlatent space of input imaging data. Subsequently, an SVM classifier was trained on the reduced latent features to predict\ntreatment outcome, as binary classes: responder(>80% seizure reduction) or non-responder (<20% seizure reduction).\nModeling results indicate that with this dataflow, the SVM can predict VNS response above chance levels (AUC = 0.73).\nGiven the lack of robust, reliable presurgical bio/markers of response to implantable VNS, this work represents a valuable\nstep toward the development of preoperative, neuroimaging-informed clinical decision making in treating pediatric epilepsy\nwith VNS. In addition to aggregating a large cohort, the authors further attempt to localize relevant anatomical and functional\ncorrelates of response by examining the model’s learned representations, informing further clinical explorations and\nhighlighting the utility of presurgical T1w MRI scans.\nThe manuscript is clinically valuable and generally well-written with most methods detailed. However, given the difficulty in\nreproducing this work (closed data due to patient privacy), there are several critical aspects of the methodology, evaluation,\nand presentation that require significant clarification and refinement to improve transparency, and reproducibility. These\nfactors along with some unsupported discussion claims detract from the main message of the manuscript and require\ntitration.\nMajor Concerns\nFeedback\n• The organizations of the manuscript and its comparisons are disjointed, and their presentation is confusing. The\nmanuscript would benefit from improved organization. For example, combine all modeling steps into one figure with subparts\nand clearly label inputs, sizes, transformations etc. Outcomes can be ported to another figure.\n• The work centers on presenting a framework for predicting VNS response called VQ-VNS, which presumably includes the\nencoding/decoding model (VQ-VAE) and then the SVM. The use of “VQ” in naming the framework leaves room for confusion\nas it’s suggestive of only the encoding/decoding aspect of modeling (Is this framework specific to only using an SVM post\nencoding?). I suggest renaming or explaining explicitly what the framework is and isn’t for clarity.\n• There are several models fit throughout the manuscript (OLS, Linear-SVM, etc.). The manuscript would benefit from a\nstructured comparison between a progression of models fit on clinical variables then latent features. For example, fitting a\nsimple logistic regression model on the clinical variables, and SVM, and SVM with anatomical information, then SVM as\nwell as other ensemble/ explainable models on latent space features. The AUC and other metrics for this should be\ncompared within one figure.\n• Throughout the manuscript avoid nonspecific references like “the model,” be specific with model type or model variation as\nit becomes confusing.\n• The lack of cross-validation or presentation of variation throughout the manuscript is concerning. Please apply these\ntechniques where possible and report outcomes with quantified variation.\n• Figure 1c includes subjects over the age of 20, can you explain their inclusion? The manuscript could benefit from a\ndeeper/more explicit discussion of the heterogenous age distribution of subjects and how this variation could affect\nindications of atrophy and findings related to serotonin networks.\n• Please revise the claims on lines 206-208. The model’s poor performance on out-of-distribution data generally does not\nsuggest any specificity to VNS, especially given that the model was pre trained on healthy non-VNS patients. This likely\nreflects domain shift or lack of representativeness in training data, and less likely biological specificity.\n• Throughout the manuscript please be explicit about dimensions and sizes where possible. For example, lines 194-195:\n“spatial dimension” could be made clearer with subscripts x,y,z.\n• Explain the choice and dimensionality of the 32-dimension feature space on line 196.\nFigures\n• The manuscript would greatly benefit from a figure illustrating the full model development process: from pretraining, to\nfeature extraction, to SVM training and prediction. Each step should also include dimensionality of data flow, especially\ninput, output, codebook, and latent volumes.\n• Figure 1g and Figure 3b both have labels for “VNS” but demonstrate modeling results for two different models. Please use\nbetter labels in figure legends to better distinguish between different types of models used.\n• Figure 2g borders on visual exaggeration. Please plot x-axis from zero and have a zoomed-in subpanel or use a different\nplot type. Please also clarify how data from this figure were acquired. For example if this is an aggregate measure across\nsubjects, please show variations across generations for each region.\nModeling\n• Explain the choice of using an SVM as the classifier for the embeddings rather than an ensemble of models or other\nalternatives (random forest etc.). I strongly suggest comparing different modeling strategies, where applicable.\n• Please include a full suite of model assessment metrics where applicable (Accuracy, F1, Precision, Recall, etc.)\n• For each subset of data (test, validation, test) please explicitly provide the distribution of responders and non-responders\nand be explicit about data imbalances and highlight mitigation strategies.\n• Add reconstruction scores to the images shown\n• In supplemental figure 3, please double check that the cases shown in a-d do in fact correspond to the inputs in panel e\n• Ensure all mathematical equations are fully defined and explained with all parameters clearly labeled.\n• In the UMAP embeddings the manifold shows some separation of data. Were those data points inspected further to\ndetermine which factors could have contributed to their separation? If so, please explain or suggest possible separation\nfactors.\n• Methodological details tend to be lacking or missing for secondary comparisons. For example, in examining subjects with\nbetween 20%-80% seizure reduction, evaluating reconstruction with abnormal imaging etc. Please provide additional details\nand organize carefully.\n• For lines 527-530, please reword for improved clarity.\n• Please clarify/detail assessment steps used for subjects between 20%-80% seizure reduction.\n• In evaluating the cortical thickness estimation, explain the choice behind using t-SNE vs other techniques. What are its\nbenefits and shortcomings?\n• Explain the choice between permuting model outcomes versus other model comparison techniques.\n• In order to improve potential clinical applicability, consider framing predictions as probability scores of response (e.g.\nP(Response|Input)) rather than binary outputs.\nData and Code\n• I strongly suggest porting the model weights over to a more stable, permanent repository rather than its current form, which\nis inaccessible. Some domain appropriate archives include OpenNeuro, Zenodo, etc.\n• This work would greatly benefit from a demonstrative Jupyter notebook. This could potentially be shown with anonymized\nor synthetically generated responder/non-responder T1w MRI scans, which are then used to demonstrate feature extraction,\nvisualization, and then prediction with the trained classifiers.\nMinor Concerns\nFigures\n• Please unionize all the font types across figures and sub panels\n• Double check that figures conform to the journal's guidelines (i.e. font sizes etc.)\n• For figure 1d, discretize colors, space them out more and add the rounded percentages to the corresponding block for\nclarity. For example, if “Development Delay” accounts for 90% of subjects please label the bar with the corresponding\npercentage. Under each label add the type of outcome (i.e. binary, continuous, etc.)\n• Figure 1d please make axes square for clearer visual comparison\n• Figure 1c, please add the sample size for each subgroup (responders/non-responders). Can be added inside the plot\nunder each group or centered under the corresponding label.\nResults\n• Add min and max age ranges where applicable", "reviewer_3": "Clinical\nResearchers developed VQ-VNS, a deep learning model that predicts the effectiveness of vagus nerve stimulation (VNS) for\nchildren with epilepsy using routine preoperative MRI scans. Traditional clinical data could not predict treatment response,\nbut VQ-VNS achieved significant predictive accuracy and highlighted serotonin-rich brain regions and connectivity\ndisruptions in non-responders. These results are clinically relevant and highly translational, offering a biologically\ninterpretable tool to guide decisions and reduce unnecessary invasive procedures.\nMethods\nPatient grouping: Participants were divided into >80% seizure reduction vs. <20%. What about patients between these cut-\noffs? Excluding ~350 patients (694/1046 included) seems a large data loss. Since you collected seizure frequency reduction\nas a continuous variable, would correlating with this continuous measure be more informative?\nClinical outcomes: In the analysis section, you note that including patients with ambiguous outcomes did not improve SVM\nperformance. This is an important observation — why not include them in the main analysis if results are unchanged?\nPlease clarify the added value of excluding them.\nClinical variables: Have magnet effects and quality of life (QOL) been addressed? Both are well-established as clinically\nrelevant (Danthine V., Boon P.), and could enrich the dataset. If available, these should be included.\nSVM classifier: Only 263 patients were included after excluding those with prior surgery or significant abnormalities. Please\ndescribe their clinical characteristics. Were they representative of the larger cohort (e.g., intellectual disability, DEEs,\nepilepsy duration)? The generalizability of the model is limited if large patient groups are excluded.\nTLE cohort: It is unclear why specificity analyses were limited to a TLE subgroup. Were these patients outside the\nCONNECTIVOS trial? Were they seizure free? Why not include other surgical subgroups (e.g., ETLE)?\nMethodological consistency: For the VNS analysis, patients with brain lesions were excluded. Is this a methodological\nrequirement, and if so, why would it not apply to the surgical cohort? Would comparison to a DBS cohort (ANT, CM) be more\nappropriate than only TLE?\nDefinitions: How was developmental delay defined, and does this overlap with DEE cases?\nParameters: The discussion states VNS was “optimally programmed.” Please define this — were parameters such as rapid\ncycling or saturation cycling considered (e.g., Norpulse study)?\nData splitting: The dataset was split 80/10/10 for training/validation/test. What is the rationale for this choice? Did results vary\nwith different splits (e.g., 50/25/25)? Please provide references or justification.\nPET-derived neurotransmitter maps: The description of PET-based receptor/transporter profiles across nine systems is\nintriguing. Were these data specifically acquired for this study, or part of a pre-existing dataset? This is not fully described in\nthe “pretraining data compilation” section.\nResults / Discussion\nThe finding of serotonergic large-scale network disruptions is novel and important. However, the lack of results for the locus\ncoeruleus–noradrenergic system is surprising, given its strong support in both human (P300, pupillometry) and animal\nstudies (Raedt, Krahl). Could the authors comment on this apparent discrepancy?\nThe statement that “beyond optimal programming, there are no associations between increased charge density and\noutcome” (ref 59) may need qualification. Recent work, including Fahoum et al. (Brain Stimulation, 2022), suggests\noptimal intensity correlates with outcomes, especially in pediatric cohorts. Please reconcile with this evidence.\nDetails\nThe section on region-based morphometry would be clearer if you specify whether this refers to VBM, SBM, or other\napproaches.\nThe statement “region-based morphometry is constrained by atlas-based parcellations…” requires a proper reference. If\nreference 15 already supports this, please move it to the end of the sentence; if not, an additional reference is needed.\nConclusion:\nIn its current form, the manuscript is technically strong and demonstrates the use of highly specialized methodologies.\nHowever, the presentation may be difficult for clinicians and readers outside the computational field to fully grasp. If the\nauthors aim for a high-impact journal and a broad audience, we encourage them to simplify and clarify the description of the\nmethods, and to place greater emphasis on the biological and clinical interpretation of their findings — particularly the role of\nthe implicated networks and neurotransmitter systems. This would considerably enhance the accessibility, readability, and\noverall impact of the work.\nSummary Methodological/Engeneering\nThe manuscript introduces the VQ-VNS framework, a methodology designed to predict and classify responses to vagus\nnerve stimulation (VNS) in pediatric populations. The framework builds on the VQ-VAE architecture and incorporates\ninterpretability mechanisms by adapting SHAP and Grad-CAM. Identifying likely responders has the potential to minimize\nthe need for invasive procedures.\nA key limitation of the study is that the proposed methodology is presented with substantial technical complexity, but with\nlimited benchmarking against alternative approaches.\nMajor Comments\n- The authors reference several alternative approaches for responder identification [11–15]. A comparative table would help\nquantify the accuracy improvements of the proposed framework, particularly in relation to other T1-weighted MRI–based\npredictive models.\n- The manuscript states: “The resulting archetypal non-responder showed reduced grey matter, white matter, and subcortical\nvolumes, together with ventricular enlargement.” If such morphological alterations are discriminative, why would automated\nsegmentation of T1-weighted images not provide an equally effective classification of responders? A comparison with\nstandard volumetric and morphometric approaches would strengthen the contribution of the proposed framework.\n- The original spin permutation procedure rotated cortical surface maps while preserving contralateral symmetry. Was this\nconstraint maintained in the present experiments? Furthermore, the original test statistic excluded non-cortical regions. Were\nsuch regions incorporated in the current analysis? These methodological details require clarification.\nMinor Comments\n- The statement that prognostic signatures are “likely embedded in subtle, distributed, and nonlinear anatomical patterns”\nwould benefit from supporting references.\n- The rationale for selecting an SVM classifier should be expanded, in comparison with other machine learning methods\n(e.g., random forests, k-nearest neighbors).\n- The contribution of individual components within the composite loss function could be clearer. Ablation experiments\ncomparing performance with and without each term could justify their inclusion.\n- The claim that “the model’s reconstructions preserved fine anatomical details” appears somewhat contradictory, since VQ-\nVAE models are known to lose fine-scale structural information. A more precise description of which features are preserved\nor lost would enhance interpretability.\n- The observation that “reconstruction fidelity was lower in subjects with significant structural abnormalities” raises the\nquestion of whether this reduction affects classification performance. Are response predictions less reliable for patients with\npronounced abnormalities?\n- Additionally, the reconstruction fidelity of the CONNECTiVOS pediatric epilepsy cohort should be reported.\n- Gradient-based saliency mapping seems to implicate more regions than those emphasized in the text (thalami, brainstem,\ninsular cortex, peri-Rolandic regions). How were the reported regions identified and selected (e.g., thresholding criteria, local\nmaxima)?\n- Differences between synthetic responders and non-responders in Figure 2F are visually subtle. Enhancing or highlighting\nthese differences would improve interpretability.\n- “The model’s attention was concentrated within brain regions enriched in serotonin transporter (5-HTT) expression”. Can\nyou quantify this result?\n- “To assess the predictive value of clinical data, a support vector machine (SVM) was trained on 80% of the dataset\n(discovery cohort, n = 550) to predict VNS outcome, with the remaining 20% held out for validation (n = 144)”. Was there an\nindependent test set? Clarification is needed.\n- “Structure-wise agreement was evaluated across four tissue classes – grey matter, white matter, subcortical grey matter,\nand ventricles – by computing the Pearson correlation between reconstructed and original volumes across subjects”. This\nonly accounts for the volume size of each tissue type and not their location. Adding the DICE scores or an equivalent would\nprovide a more comprehensive assessment.\n- “Atrophy networks were defined by combining endpoint maps across all HCP subjects”. Why were endpoints chosen\ninstead of complete streamlines?\n- A visualization of the cortical atrophy map should be included for completeness.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We thank the reviewer for their enthusiasm for this work and for recognizing its timeliness and\nclinical utility. We appreciate the valuable constructive feedback they provided, which has\nsignificantly strengthened our report.\nMajor Points:\n1. The methodology in this paper is interesting as an application of modern machine learning\ntools. However, it would also be good to know if simpler statistical analyses of basic measures\nare similarly informative (e.g. volume/surface area/thickness measures from freesurfer). Looking\nat figure 3c, perhaps thalamic volume alone is already somewhat informative?\nThe reviewer’s comment regarding the potential informativeness of simpler statistical analyses\nand more basic imaging metrics is well-received. Indeed, it is the inadequacy of these models to\nconsider multi-dimensional representations and interactions that inspired the current\nrepresentation learning approach to VNS outcome prediction. In the revised manuscript, we have\nincluded our analysis of basic imaging-related metrics using conventional modelling techniques,\nwhich further highlight the value of novel approaches like VQ-VNS.\nMore specifically, we conducted additional analyses using standard cortical thickness and\nsubcortical volume features derived from FreeSurfer. Each measure was modelled using a\ngeneralized linear model to test its association with VNS outcome, while controlling for age and\nintracranial volume. Interestingly, the thickness of the left posterior cingulate cortex (also\nfeaturing prominently in the VQ-VNS salience map) demonstrated significant associations with\nresponse to VNS prior to FDR correction. No cortical or subcortical morphometric features\nsurvived FDR correction (Supplementary Table 3-4). Specifically, the thalamus – measured as a\nsingular structure – did not differ in volume among responders and non-responders, despite being\nthe most prominent feature in the VQ-VNS salience map. These findings confirm that\nconventional morphometric measures carry insufficient discriminative information to predict\nVNS outcome, reinforcing the need to employ deep representation learning to capture more\ncomplex, distributed, and interacting MRI features.\nFor completeness, we also applied t-distributed stochastic neighbour embedding (t-SNE) to all\ncortical vertex values and visualized embeddings by site, outcome, and age (Supplementary Fig.\n9). Comparison of atrophy size and spatial distribution between responders and non-responders\nalso demonstrated no significant group differences (Supplementary Figure 10).\nThese findings are concordant with the central thesis of the paper, which is that latent space\nrepresentations based on generative models designed to learn compact representations of brain\npatterns are useful and necessary in understanding complex associations with VNS outcomes.\nWe’ve mentioned as such in the manuscript as well.\n“To study whether structural brain differences could predict outcome, we developed VQ-VNS, a\ndeep representation learning model to identify VNS responders based on standard T1w\nsequences (Fig. 2). Because each T1w MRI scan contains over six million voxels – far exceeding\nthe sample sizes of most pediatric epilepsy datasets (typically < 100 patients per center)1 –\ndimensionality reduction is critical to mitigate overfitting and enable effective model training.\nHowever, conventional approaches to dimensionality reduction, such as region-based\nparcellation, are insufficient to model the relationship between structural alterations and\nresponse following VNS (Supplementary Tables 3-4). Thus, the VQ-VNS model consists of two\nstages: 1) a VQ-VAE2 architecture compresses the high-dimensional MRI data into a compact\nlatent representation, and (2) a support vector machine (SVM) classifies patients based on these\nlatent features.”\nSupplementary Table 3. Generalized linear model of Cortical Thickness to VNS Response\nRegion Coef. P P\nfdr\nL R L R L R\nBanks of\nsuperior sulcus 0.45 0.16 0.39 0.75 0.92 0.93\nCaudal anterior cingulate 0.21 0.57 0.65 0.23 0.93 0.92\nCaudal middle frontal 0.38 1.24 0.48 0.04 0.92 0.80\nCuneus -0.11 0.21 0.83 0.69 0.93 0.93\nEntorhinal 0.47 -0.12 0.17 0.72 0.92 0.93\nFrontal pole 0.26 0.15 0.35 0.58 0.92 0.92\nFusiform 0.67 0.35 0.28 0.56 0.92 0.92\nInferior parietal 0.58 0.98 0.30 0.10 0.92 0.92\nInferior temporal 0.44 0.34 0.40 0.52 0.92 0.92\nInsula 0.41 0.13 0.50 0.82 0.92 0.93\nIsthmus cingulate 0.79 0.15 0.13 0.78 0.92 0.93\nLateral occipital 0.48 -0.03 0.41 0.96 0.92 0.99\nLateral orbitofrontal 0.68 -0.09 0.22 0.87 0.92 0.94\nLingual -0.08 0.02 0.89 0.97 0.94 0.99\nMedial orbitofrontal 0.1 0.17 0.84 0.74 0.93 0.93\nMiddle temporal 0.36 1.21 0.51 0.05 0.92 0.80\nParacentral 0.35 0.13 0.51 0.81 0.92 0.93\nParahippocampal 0.21 -0.47 0.58 0.36 0.92 0.92\nParsopercularis 0.34 0.4 0.61 0.49 0.93 0.92\nParsorbitalis -0.09 0.51 0.84 0.26 0.93 0.92\nParstriangularis 1.15 0.10 0.05 0.85 0.8 0.93\nPericalcarine -0.3 -0.37 0.55 0.49 0.92 0.92\nPostcentral 0.19 0.00 0.76 1.00 0.93 1.00\nPosterior cingulate 1.31 0.79 0.04 0.21 0.8 0.92\nPrecentral 0.18 0.46 0.77 0.41 0.93 0.92\nPrecuneus 0.20 0.36 0.75 0.56 0.93 0.92\nRostral anterior cingulate 0.87 0.41 0.09 0.37 0.92 0.92\nRostral middle frontal 0.67 0.54 0.29 0.36 0.92 0.92\nSuperior frontal 0.11 0.73 0.83 0.19 0.93 0.92\nSuperior parietal -0.15 -0.06 0.80 0.92 0.93 0.96\nSuperior temporal 0.56 0.41 0.34 0.44 0.92 0.92\nSupramarginal 0.38 0.61 0.50 0.28 0.92 0.92\nTemporal pole 0.26 0.38 0.38 0.19 0.92 0.92\nTransverse temporal -0.49 0.19 0.26 0.64 0.92 0.93\nSupplementary Table 4. Generalized linear model of Subcortical Volumes to VNS Response\nRegion Coefficient P-value uncorrected P-value FDR corrected\nLeft Right Left Right Left Right\nAccumbens area 0.13 0.05 0.91 0.96 0.99 0.99\nAmygdala -0.34 -0.13 0.51 0.82 0.99 0.99\nCaudate -0.26 -0.3 0.33 0.26 0.99 0.99\nCerebellum Cortex 0.01 0.00 0.67 0.86 0.99 0.99\nHippocampus 0.28 0.30 0.27 0.20 0.99 0.99\nPallidum 0.17 0.03 0.72 0.94 0.99 0.99\nPutamen 0.00 -0.16 0.99 0.40 0.99 0.99\nThalamus -0.04 -0.13 0.80 0.41 0.99 0.99\nVentralDC -0.04 0.06 0.91 0.86 0.99 0.99\nBrain-Stem -0.11 0.09 0.99\nSupplementary Figure 9: Site Effect Analysis. (a) t-SNE embedding of all cortical thickness vertices and eTIV\ndata projected in 2-dimensions labelled by site. There is no visible clustering of the data by site, suggesting the\nabsence of batch effects related to image acquisition. (b) t-SNE embeddings visualized and labelled by VNS\nresponse on the left show no clear separation between the groups. (c) t-SNE embeddings labelled by age are shown\non the right and indicate an anti-correlation with the first t-SNE component. (d) The relationship between age and\nthe first t-SNE component is visualized, demonstrating a significant anti-correlation (r = -0.45, p = 6.7e-14),\nsuggesting that age is a major contributor to variability in cortical thickness data.\nSupplementary Figure 10: Size and Spatial Distribution of Atrophy. (a) Distribution of atrophy size across all\nsubjects, relative to the 95th percentile and the 5% exclusion threshold. (b) Comparison of atrophy size between\nresponders and non-responders shows no significant group difference. (c) Spatial distribution of atrophy across the\nfull cohort. (d) A small focal cluster of atrophy is preferentially observed in non-responders, with no reciprocal\npattern in responders.\n2. If I understood correctly, the 39 patients in the withheld test data was chosen randomly. This\nwithheld data of 39 patients represents only one of the possible sets of 39 patients that could be\nheld-out. It would give more confidence in the approach if this could be permuted many times\nand error bars plotted on figure 3b, along with confidence intervals for AUCs.\nWe thank the reviewer for this thoughtful suggestion. As requested, we performed repeated-\nholdout cross-validation by permuting the train/test split 1000 times and computing AUCs across\nthese permutations. The resulting distribution of AUCs confirmed that latent features predict\nVNS response at AUCs significantly above chance (p < 0.001; mean = 0.613).\nImportantly, our original withheld dataset did not represent an outlier within this randomly\npermuted distribution (p ≈ 0.2, based on 1,000 random splits). These results support the stability\nof our model’s ability to predict VNS outcomes beyond what would be expected by chance\nalone.\nWe certainly appreciate that both single holdout evaluation and repeated-holdout cross-validation\nare valid approaches to assessing generalization, each with distinct advantages and drawbacks.\nThe single held-out test set provides a stringent measure of external generalization to unseen\ndata, while repeated-holdout validation offers a measure of stability across multiple random\npartitions. We have now presented both approaches in the revised manuscript.\n“To evaluate the predictive utility of the latent space for response following VNS, we\ntrained a linear SVM classifier (Fig. 2e) using quantized latent vectors from the frozen VQ-VAE\nencoder. To reduce dimensionality, each four-dimensional latent tensor (x, y, and z spatial\ndimensions × feature channels) was spatially averaged across the x, y, and z axes, yielding a\nsingle 32-dimensional3 feature vector per subject. Each element of this vector therefore captures\nthe global magnitude of the respective latent feature learned from the underlying T1w image. We\nthen applied two-sided t-test-based feature selection on the training data (n = 199) to retain 27\nof the 32 latent features (Supplementary Table 7), which were used to train the SVM. We used an\nSVM for this purpose as it demonstrated the strongest predictive performance among multiple\nclassifiers (Supplementary Fig. 6). To assess the stability and significance of the learned\nfeatures, we first trained 1000 SVM classifiers on randomly partitioned subsets of the VNS\ncohort. The mean AUC across models was significantly greater than that expected under the\nnull distribution (mean AUC = 0.61; t = 28.6, p = 0.001; Fig. 3f), suggesting that the latent\nrepresentations carried reproducible predictive signal. We then trained a final SVM classifier\nusing the training (n = 199) and validation (n = 25) subsets of the VNS cohort and tested on an\nindependent held-out set (n = 39). The resulting model achieved an AUC of 0.73 in predicting\nresponse to VNS (permutation p = 0.007; Fig. 3g; Supplementary Fig. 7; Supplementary Table\n8). The model maintained its predictive capacity across diverse age groups and treatment\ncenters (Supplementary Fig. 7C). Model performance remained robust with the inclusion of\nsubjects with ambiguous outcomes binarized at the 50% seizure reduction threshold (n = 53;\nAUC = 0.67, p = 0.02; Supplementary Fig. 8).”\nFigure 3. Performance of the VQ-VNS Quantized Autoencoder Model for T1w Reconstruction and\nClassifier for Predicting VNS Outcomes. (f) 1000 models trained using repeated hold-out cross\nvalidation perform significantly better than chance. (g) Receiver operating characteristic (ROC) curve\nshowing the predictive performance of the SVM classifier on the held-out test set (AUC = 0.73, p =\n0.007), which was specific to VNS and did not generalize to outcomes following temporal lobectomy\n(AUC = 0.57, p = 0.12).\n3a Following on from point 2 regarding test data. I think it would be worthwhile to test\ngeneralisability across data from different hospitals. I.e. hold out data from one hospital and train\non the others, then report test accuracy on the withheld hospital. This would be useful to know\nsince it will tell us if the model can be used/transferred to other hospitals, which would be\nclinically important.\nTo examine the generalizability of this model across hospitals, we have trained additional\nclassifiers with one hospital held out as the test set. The resulting performance metrics are\nsummarized below and have been included in the revised manuscript.\nWhile this analysis provides an additional perspective on model generalizability, it is constrained\nby multiple limitations inherent to our dataset. Most notably, the unequal contribution of data\nacross sites introduces variability in both training and testing conditions: omitting large sites\nsuch as SickKids substantially reduces the training sample and impairs model learning, whereas\nsmaller sites such as UCLA yield limited test samples where performance estimates may be\nunreliable or overly influenced by individual cases. Furthermore, excluding certain hospitals\nleads to class imbalance that affects both model training and the stability of evaluation metrics.\nTo partially mitigate these issues, we combined smaller sites (Indiana, Miami, UBC, Montreal,\nSt. Louis, Texas, and Wisconsin) into pooled test sets to provide a more stable evaluation.\n“To evaluate the classifier’s ability to generalize across acquisition sites, we performed\nsite hold-out validation, in which all images from one site were excluded during training and\nused exclusively for testing (Supplementary Table 9). To mitigate the effect of evaluating the\nmodel on sites with small sample sizes, these sites were combined to create a single evaluation\nset. The classifiers achieved above chance accuracy (0.61 to 0.67), suggesting site-invariant\nstructural features predictive of VNS response.”\nSupplementary Table 9. Site Hold-Out Model Performance\nTest site Test Size Accuracy Precision Recall F1 Specificity\nPittsburgh 33 0.67 0.67 0.71 0.69 0.62\nUBC and 23 0.65 0.50 0.38 0.43 0.80\nMontreal\nWisconsin 54 0.61 0.63 0.66 0.64 0.56\nIndiana, St.\nLouis\nTexas and 55 0.62 0.64 0.78 0.70 0.39\nMiami\nUCLA 13 0.62 0.62 0.71 0.67 0.50\n3b. Continuing from point 3a, are there any hospitals with data to which the model/approach\ndoes not generalise to very well? If so why? (perhaps they have older/younger demographics?\nPoorer scan quality?)\nTo examine potential demographic or site-specific effects on model performance, we compared\nage and site distribution of correct and incorrect predictions (Supplementary Figure 6C). Model\naccuracy was consistently above chance across participating sites with four or more patients.\nFurther, the age distribution did not differ between correct and incorrect predictions.\nImportantly, the deep learning model was trained on a substantially larger multi-institutional\ndataset of MRIs from across the lifespan, which provided representation from diverse scanners\nand image qualities. This training helps mitigate the risk that performance differences are driven\nby site-specific technical factors or demographics.\nWe have added this information to the manuscript as below,\n“We then trained a final SVM classifier using the training (n = 199) and validation (n = 25)\nsubsets of the VNS cohort and tested on an independent held-out set (n = 39). The resulting\nmodel achieved an AUC of 0.73 in predicting response to VNS (permutation p = 0.007; Fig. 3g;\nSupplementary Fig. 7; Supplementary Table 8). The model maintained its predictive capacity\nacross diverse age groups and treatment centers (Supplementary Fig. 7C). Model performance\nremained robust with the inclusion of subjects with ambiguous outcomes binarized at the 50%\nseizure reduction threshold (n = 53; AUC = 0.67, p = 0.02; Supplementary Fig. 8).”\nSupplementary Figure 7C: VQ-VNS model performance did not differ across age or sites.\nBoxplots (left) show the distribution of age at scan for correctly and incorrectly predicted\nsubjects, demonstrating no systematic age bias. Bar plots (right) depict distribution of correct\nand incorrect predictions across participating sites. Performance was consistently above chance\nacross all acquisition sites.\n4. The authors focus their evaluation mainly on AUC. This is important, and I have no problem\nwith this. I think it might also be worthwhile reporting if the specificity at 100% sensitivity and\nvice versa. I think this would be worthwhile because this may be useful for either ruling out\ndefinite non-responders, or perhaps as a screening tool for identifying/prioritising likely\nresponders.\nIn the revised manuscript, we have reported the sensitivity and specificity of the VQ-VNS\nmodel.\nSupplementary Table 8. Model Assessment Metrics\nModel Accuracy Precision Recall F1 Specificity\nClinical Data 0.53 0.53 0.51 0.52 0.54\nTest set 0.72 0.70 0.74 0.72 0.70\n50% 0.66 0.67 0.67 0.67 0.65\nBinarization\nTLE 0.56 0.69 0.62 0.65 0.44\nPer the reviewer’s suggestion, we also examined the model performance at extreme operating\npoints. At thresholds achieving ~100% sensitivity, the specificity was 5%, reflecting the large\nnumber of non-responders misclassified when attempting to capture all responders. Conversely,\nat thresholds achieving ~100% specificity, sensitivity fell to 0%, as no responders were\nidentified. These results suggest that the model is not better suited to rule out definite non-\nresponders.\n5. The regions identified as having a “Higher” feature value for prediction using the gradient-\nbased saliency mapping are not provided with a quantitative backing. Providing a numerical\nunderpinning to justify the chosen regions would provide greater confidence and understanding\nof whether these regions are uniquely involved or if the whole brain is still required. N.b. this\ncould potentially be partly addressed through point 1 above.\nIn the revised manuscript, we have added quantitative values (natural log-transformed Grad-\nCAM saliency values) to the figure to provide a numerical underpinning for the relative feature\nimportance across the whole brain. The values represent normalized importance scores and\ntherefore do not have intrinsic or absolute meaning. Instead, they indicate the relative\ncontribution of each region to the model’s prediction. As noted in our response to point 1 above,\nno single region alone was sufficient to predict VNS response after applying FDR correction,\nsuggesting that distributed patterns across multiple regions are required. Further, we have\nincluded below an illustration of the highest-ranking features in Figure R1, thresholded at a\nsaliency value of 0.3, representing the top 10% of values.\nFigure 4: Model Salience Aligns with Large-Scale Serotonergic Brain Networks Disrupted in\nVNS Non-Responders. (a) Gradient-based saliency mapping highlighted the most predictive\nregions, including the thalami, brainstem, and peri-Rolandic cortices.\nFigure R1. Model salience map thresholded at 0.3, representing the top 10% of values.\nMinor Points:\nReference 8 and reference 32 are the same (I have not checked, but it is possible this occurs more\nthan once).\nWe thank the reviewer for noting the duplication of references. We have carefully reviewed the\nfull list of references and corrected the duplicate entries.\nIt is surprising given prior literature that you find so many individuals with >=80% and <=20%\nseizure reduction. Historically these rates appeared far lower. Do you have any knowledge as to\nwhy this combined North American cohort appears to have significantly more individuals at both\nextremities of response?\nVNS outcomes are often reported dichotomously using a 50% threshold, but a finer examination\nof outcomes has often shown a bimodal distribution. For instance, Majkaowska-Zwolinksa et al.\nobserved that among 50% responders at 12 months, 41% achieved > 75% seizure reduction,4\nwhile Orosz et al. found that 40% of children had <25% reduction at 12 months.5 Although this\npattern has been underrecognized due to reporting conventions, our findings are consistent with\nprior pediatric studies that indicate VNS response may cluster at the extremes rather than\ndistribute evenly across categories.\nIt is not identified how seizure reduction was calculated. Assumedly, this is through seizure\ndiaries or discussed at follow-up appointments which should be clarified.\nSeizure reduction was determined from parent-reported seizure frequency. We have updated the\nMethods to reflect this information:\n“For each participant, we collected continuous data on age at VNS implantation, duration of\nepilepsy, and percent seizure reduction at last follow-up based on parent-reported seizure\nfrequency.”\nWithin the methods you refer to training cohorts both as training and discovery cohorts\nconsistent naming would be clearer.\nTo enhance clarity, we have removed the term “discovery cohort” from the revised manuscript.\n“Subjects with significant anomalies, prior neurosurgical intervention (Supplementary Fig. 3a-d),\nor conflicting outcome assessments, were excluded. After filtering, the final dataset comprised\n263 subjects”. The exact number dropped for each of these reasons would be interesting. A\nPRISMA flow diagram would be helpful as a supplementary figure.\nWe thank the reviewer for this suggestion. We have now included the exact number of subjects\nexcluded for each criterion in Supplementary Figure 5.\nSupplementary Figure 5: Participant inclusion and exclusion flow diagram.\nfigure 4d; top row, middle panel, the Insula arrow terminates in the putamen at the moment.\nThank you for noting this. The revised figure attached below.\nFigure 4: Model Salience Aligns with Large-Scale Serotonergic Brain Networks Disrupted in VNS\nNon-Responders. (a) Gradient-based saliency mapping highlighted the most predictive regions,\nincluding the thalami, brainstem, insula, and peri-Rolandic cortices. (b) Methodology for generating null\nsalience maps. Starting from each subject’s latent salience map, a random rotational affine\ntransformation is generated and applied to all subjects. The rotated maps are up-sampled and averaged\nacross subjects to create a null salience map. This process was repeated 1,000 times to construct an\nempirical null distribution. (c) Model salience was concentrated within brain regions enriched in 5-HTT\nexpression (p = 0.04). Lesion network mapping using normative (d) structural and (e) functional\nconnectomic datasets showed that regional atrophy in non-responders disrupted large-scale networks\nthat aligned spatially with the salience map (p < 0.001 and p = 0.04, respectively).\n“ Each four-dimensional latent tensor (spatial dimensions × feature channels) was spatially\naveraged across the x, y, and z axes, yielding a single 32-dimensional feature vector per subject.\n“ I struggled to follow/understand the rationale for the averaging across the spatial axes xyz.\nPlease can you explain this in more detail? Perhaps in an expanded supplementary methods.\nThe averaging across spatial dimensions was applied as a dimensionality reduction step. The\nlatent representation produced by the encoder is a four-dimensional tensor (11 × 13 × 11 × 32),\nwhich remains too large to use directly for subject-level prediction. Each of the 32 channels\nwithin these subject-level latent maps captured distinct tissue-related features (e.g., CSF, gray\nmatter, white matter), albeit in a distributed and abstract manner (Fig. 3e). By spatially averaging\nacross the x, y, and z dimensions, we summarized each channel into a single scalar value,\nyielding a 32-dimensional feature vector per subject. This compact representation reflects the\nmagnitude of each latent feature across the entire brain volume and is suitable for use in the\nlinear SVM classifier.\n“To evaluate the predictive utility of the latent space for response following VNS, we\ntrained a linear SVM classifier (Fig. 2e) using quantized latent vectors from the frozen VQ-VAE\nencoder. To reduce dimensionality, each four-dimensional latent tensor (x, y, and z spatial\ndimensions × feature channels) was spatially averaged across the x, y, and z axes, yielding a\nsingle 32-dimensional3 feature vector per subject. Each element of this vector therefore captures\nthe global magnitude of the respective latent feature learned from the underlying T1w image.”\nFigure 2: Training of the VQ-VNS Model. (a) Preprocessing pipeline of T1-w MRI for pretraining of the\nVQ-VNS Quantized Autoencoder Model for T1w reconstruction. Each raw T1w MRI image undergoes\nhybrid brain extraction (SynthStrip and SynthSeg, reorientation to standard RAS orientation, 1 mm\nisotropic resampling, bias field correction, and rigid registration to the MNI 152 template followed by\ncropping to brain margins. (b)-(c) VQ-VAE architecture for T1w MRI reconstruction, comprising a\nconvolutional encoder, vector quantization bottleneck, and convolutional decoder. (d) Composition of the\ndatasets used for VQ-VAE reconstruction training, which included healthy subjects and persons with\nepilepsy. (e) Schematic of the classification pipeline: the spatial mean of the VQ-VAE quantized latent is\ncomputed, feature selected, and used as input to a linear SVM for VNS response prediction.\nfigure 2g -I did not follow what is relative volume. relative to what? The responder blue bars on\nthe chart are all at 100%, so I assume it’s relative to the responder. Perhaps it would be helpful to\nplot relative to healthy control and show a synthetic control in figure 2f for reference\nIn the original submission, “relative volume” in Figure 2g indeed referred to differences\ncalculated relative to the synthetic, archetypal responder, which is why the responder bars were\nall at 100%. To improve clarity, we have revised the plot to show the differences between the\nsynthetic, archetypal responder and non-responder. Panel g highlights the volumetric differences\nof the synthetic responder and non-responder shown in Figure 2f.\nFigure 3: Performance of the VQ-VNS Quantized Autoencoder Model for T1w Reconstruction and\nClassifier for Predicting VNS Outcomes. (e) The synthetic archetypal non-responder showed reduced\ngrey matter, white matter, and subcortical grey matter volumes, and ventricular enlargement relative to\nthe synthetic archetypal responder.\nI thought the validation in the TLE cohort was a brilliant idea. I was wondering what the\ninterpretation would be if the model had generalised to the TLE cohort. Do you think that would\nindicate shared mechanisms of response to VNS and also response to TLE surgery? As such,\nshould we conclude that given the model did not generalise to TLE the mechanisms differ?\nWe thank the reviewer for this thoughtful comment. We agree that had the model generalized to\nthe TLE surgical cohort, this could have suggested that there are shared predictors or\nmechanisms underlying favorable outcomes to both VNS and resective surgery. The lack of\ngeneralization may suggest that the features predictive of VNS response may be distinct from\nthose driving surgical success. Indeed there are some studies that link subcortical (i.e. thalamic)\nconnectivity to outcomes following TLE.6 The current results are consistent with the possibility\nthat the mechanisms supporting VNS response differ from those related to TLE outcomes.\nReviewer #2 (Remarks to the Author):\nRecommendation:\nMajor Revisions\nTitle:\nA Deep Representation Learning Model to Predict Response to Vagus Nerve Stimulation\nOverview and general recommendation:\nThis work investigated the use of T1-weighted MRI data for predicting clinical response to\nimplantable vagus nerve stimulation (VNS) therapy in individuals with drug-resistant epilepsy.\nThe authors leveraged their access to a network of clinical neuroimaging databases to aggregate\na large cohort of preoperative and postoperative (resections etc.) T1w-MRI scans, which were\nthen used for modeling development. A vector quantized variational autoencoder was used to\nmodel the latent space of input imaging data. Subsequently, an SVM classifier was trained on the\nreduced latent features to predict treatment outcome, as binary classes: responder (>80% seizure\nreduction) or non-responder (<20% seizure reduction).\nModeling results indicate that with this dataflow, the SVM can predict VNS response above chance levels (AUC = 0.73).\nGiven the lack of robust, reliable presurgical bio/markers of response to implantable VNS,\nthis work represents a valuable step toward the development of preoperative, neuroimaging-informed\nclinical decision making in treating pediatric epilepsy with VNS. In addition to aggregating a\nlarge cohort, the authors further attempt to localize relevant anatomical and functional correlates\nof response by examining the model’s learned representations, informing further clinical\nexplorations and highlighting the utility of presurgical T1w MRI scans.\nThe manuscript is clinically valuable and generally well-written with most methods detailed.\nHowever, given the difficulty in reproducing this work (closed data due to patient privacy), there\nare several critical aspects of the methodology, evaluation, and presentation that require\nsignificant clarification and refinement to improve transparency, and reproducibility. These\nfactors along with some unsupported discussion claims detract from the main message of the\nmanuscript and require titration.\nWe thank the reviewer for their thorough and thoughtful evaluation of our manuscript. We\nappreciate the recognition of clinical value and novelty of our approach, as well as the\nconstructive feedback aimed at improving the clarity, transparency, and reproducibility of our\nwork. We have carefully considered each comment and have revised the manuscript accordingly\nto enhance methodological detail, strengthen the justification of our claims, and improve overall\npresentation.\nMajor Concerns\nFeedback\n• The organizations of the manuscript and its comparisons are disjointed, and their presentation is\nconfusing. The manuscript would benefit from improved organization. For example, combine all\nmodeling steps into one figure with subparts and clearly label inputs, sizes, transformations etc.\nOutcomes can be ported to another figure.\nThank you for this suggestion. We have reorganized the manuscript to improve clarity and flow.\nSpecifically, we now introduce an overview of the full modeling pipeline to provide readers with\na clear, high-level understanding of the approach (Figure 2). Subsequently, we walk the reader\nthrough the two key components of the model, the VQ-VAE and the SVM, in a logical sequence.\nTo ensure that the VQ-VAE latent representations captured the necessary anatomical\ninformation, we evaluated the full reconstruction pipeline by comparing input and reconstructed\nimages. We then assessed whether these latent representations encoded meaningful group\ndifferences by training an SVM to distinguish responders from non-responders. With this\nintegrated structure, the manuscript now follows a clear narrative from motivation to model\ndevelopment, validation, and finally to model interpretation.\nFigure 2: Training of the VQ-VNS Model. (a) Preprocessing pipeline of T1-w MRI for pretraining of the\nVQ-VNS Quantized Autoencoder Model for T1w reconstruction. Each raw T1w MRI image undergoes\nhybrid brain extraction (SynthStrip and SynthSeg, reorientation to standard RAS orientation, 1 mm\nisotropic resampling, bias field correction, and rigid registration to the MNI 152 template followed by\ncropping to brain margins. (b)-(c) VQ-VAE architecture for T1w MRI reconstruction, comprising a\nconvolutional encoder, vector quantization bottleneck, and convolutional decoder. (d) Composition of the\ndatasets used for VQ-VAE reconstruction training, which included healthy subjects and persons with\nepilepsy. (e) Schematic of the classification pipeline: the spatial mean of the VQ-VAE quantized latent is\ncomputed, feature selected, and used as input to a linear SVM for VNS response prediction.\n• The work centers on presenting a framework for predicting VNS response called VQ-VNS,\nwhich presumably includes the encoding/decoding model (VQ-VAE) and then the SVM. The use\nof “VQ” in naming the framework leaves room for confusion as it’s suggestive of only the\nencoding/decoding aspect of modeling (Is this framework specific to only using an SVM post\nencoding?). I suggest renaming or explaining explicitly what the framework is and isn’t for\nclarity.\nIn response to the previous comment, we have explicitly defined VQ-VNS as a two-stage\nframework that integrates (1) a Vector Quantized Variational Autoencoder (VQ-VAE) to learn\ncompact, anatomically meaningful latent representations from T1-weighted MRI, and (2) a\nSupport Vector Machine classifier trained on these latent representations to predict clinical\nresponse to VNS.\n“To study whether structural brain differences could predict outcome, we developed VQ-VNS, a\ndeep representation learning model to identify VNS responders based on standard T1w\nsequences (Fig. 2). Because each T1w MRI scan contains over six million voxels – far exceeding\nthe sample sizes of most pediatric epilepsy datasets (typically < 100 patients per center)1 –\ndimensionality reduction is critical to mitigate overfitting and enable effective model training.\nHowever, conventional approaches to dimensionality reduction, such as region-based\nparcellation, are insufficient to model the relationship between structural alterations and\nresponse following VNS (Supplementary Tables 3-4). Thus, the VQ-VNS model consists of two\nstages: 1) a VQ-VAE2 architecture compresses the high-dimensional MRI data into a compact\nlatent representation, and (2) a support vector machine (SVM) classifies patients based on\nthese latent features.”\n• There are several models fit throughout the manuscript (OLS, Linear-SVM, etc.). The\nmanuscript would benefit from a structured comparison between a progression of models fit on\nclinical variables then latent features. For example, fitting a simple logistic regression model on\nthe clinical variables, and SVM, and SVM with anatomical information, then SVM as well as\nother ensemble/ explainable models on latent space features. The AUC and other metrics for this\nshould be compared within one figure.”\nWe thank the reviewer for this thoughtful suggestion. While we agree that a structured\nprogression of models could be informative in some contexts, this was not the case in our study.\nWe did fit a model using only clinical variables, but its performance was not predictive of\noutcome and therefore was not pursued further. For this reason, we focused our analyses on\nmodels based on latent features extracted from T1w images, which provided meaningful\npredictive signal.\nIn Supplementary Table 8, we now include assessment metrics for the clinical dataset, the held-\nout test set, the additional patients with intermediate outcomes, and the TLE cohort for\ncomparison.\nSupplementary Table 8. Model Assessment Metrics\nModel Accuracy Precision Recall F1 Specificity\nClinical Data 0.53 0.53 0.51 0.52 0.54\nTest set 0.72 0.70 0.74 0.72 0.70\n50% 0.66 0.67 0.67 0.67 0.65\nBinarization\nTLE 0.56 0.69 0.62 0.65 0.44\n• Throughout the manuscript avoid nonspecific references like “the model,” be specific with\nmodel type or model variation as it becomes confusing.\nWe appreciate the reviewer’s suggestion to avoid nonspecific terminology to improve clarity.\nWe have revised the manuscript throughout to consistently specify the exact model type or\nvariation being referenced (e.g. clinical SVM, VQ-VAE, or VQ-VNS).\n• The lack of cross-validation or presentation of variation throughout the manuscript is\nconcerning. Please apply these techniques where possible and report outcomes with quantified\nvariation.\nIn addition to the holdout validation performed originally, we have now revised the manuscript\nto include repeated hold-out cross validation. The mean AUC across models was significantly\ngreater than that expected under the null distribution (mean AUC = 0.61; t = 28.599, p = 0.001;\nFig. 3f), suggesting that the latent representations carried reproducible predictive signal. Our\noriginal holdout AUC did not represent an outlier within the permutation distribution (p ≈ 0.2,\nbased on 1,000 random splits), indicating that the held-out test set performance is consistent with\nthe expected range of outcomes under random partitioning. These results support the stability of\nour model’s predictive performance and reinforce the validity of the originally reported holdout\nevaluation.\nWe would also like to note that cross-validation is generally employed when the available\nsample size is limited and it is not feasible to set aside a dedicated test set.7,8 When sufficient\ndata are available, evaluation on a fixed held-out test set is typically preferred as it best reflects\ngeneralization to unseen patients. In our study, we originally adopted this approach by holding\nout 39 patients as an independent test set. We now complement this with cross-validation results.\n“To assess the stability and significance of the learned features, we first trained 1000 SVM\nclassifiers on randomly partitioned subsets of the VNS cohort. The mean AUC across models\nwas significantly greater than that expected under the null distribution (mean AUC = 0.61; t =\n28.6, p = 0.001; Fig. 3f), suggesting that the latent representations carried reproducible\npredictive signal. We then trained a final SVM classifier using the training (n = 199) and\nvalidation (n = 25) subsets of the VNS cohort and tested on an independent held-out set (n =\n39).”\n• Figure 1c includes subjects over the age of 20, can you explain their inclusion?”\nParticipants up to age 25 were included in order to capture the full developmental trajectory\nrelevant to epilepsy in pediatric and young adult populations. White matter microstructure and\nthalamocortical pathways continue to mature through the mid-twenties. Including young adults\nallowed us to better represent the continuum of structural variation and seizure outcomes. The\nage range is now explicitly reported under Clinical data do not predict response after VNS.\n“We compiled the largest database of children to undergo VNS (n = 1046) from 11 North\nAmerican centers through the Connectomics and Vagus Nerve Stimulation Outcome Study\n(CONNECTiVOS, Fig. 1a).28 Participants were deemed responders if seizure frequency was\nreduced by ≥80% and non-responders if ≤20% (n = 694, Fig. 1b) to reflect the bimodal\ndistribution of seizure outcomes and mitigate the effects of ambiguous intermediate outcomes.29\nThe mean [SD; range] age was 10.6 [4.7; 0.8-24.2] years, 54% male and included 51%\nresponders (Fig. 1b-c). Participants up to age 25 were included to capture the full\ndevelopmental trajectory of brain maturation amongst pediatric and young adult populations\nthat undergo VNS.”\nThe manuscript could benefit from a deeper/more explicit discussion of the heterogenous age\ndistribution of subjects and how this variation could affect indications of atrophy and findings\nrelated to serotonin networks.\nWe thank the reviewer for this important comment. Given that our study involves a pediatric\npopulation with substantial developmental variability, we carefully accounted for age in every\naspect of our analyses. For this reason, the VQ-VAE representation learning model was trained\non a dataset of T1w MRI images across the lifespan, ranging from 1 to 89 years old. In\nSupplementary Figure 8, we demonstrate that variability in cortical thickness across the cohort\nwas primarily driven by age, rather than site or outcome, with a significant age-related\ncomponent (r = -0.45, p = 6.7e-14). This confirms that age is a major contributor to variability in\ncortical thickness data, but these features are not predictive of response to VNS.\nTo directly assess potential effects of age-related cortical changes on VNS outcomes, we\nexamined the relationship between regional cortical thickness and treatment response.\nFreeSurfer-derived cortical thickness measures were entered into general linear models\npredicting VNS response, with age and intracranial volume included as covariates\n(Supplementary Table 3). To assess stability across age, we re-ran these analyses after\ndichotomizing the cohort by median age (Table R1 shown below). Across both the full and age-\nsplit samples, cortical thickness was not associated with VNS response in any region of the\nbrain, indicating that age-related variability does not explain group differences.\nIn addition, our cortical atrophy analysis explicitly accounted for age, as atrophy maps were\ncomputed relative to age- and sex-matched healthy controls within the lesion network mapping\nframework. This design ensures that apparent cortical thinning reflects pathological deviation\nfrom expected developmental trajectories rather than normative developmental effects.\nWe have now included a more explicit discussion on how age could affect indications of atrophy\nand the implication it has on our findings.\n“The VQ-VNS model demonstrated predictive accuracy across the pediatric age span. In the\nCONNECTiVOS cohort, we demonstrate that a child’s chronological age at implantation is\nnot associated with VNS outcomes. Rather, a pattern of brain atrophy characterized by the\nVQ-VNS salience map predicts response across development (Supplementary Fig. 7c). \nConvergence of the salience map with the LNM analysis – wherein brain atrophy is\nparametrized relative to normative age- and sex-matched developmental trajectories – supports\na framework where structural brain changes drive VNS response.”\nTable R1. Generalized linear model of FreeSurfer Cortical Thickness to VNS Response\nDichotomized by Age\nYounger children (≤ 9 years) Older children (> 9 years\nCoefficients Puncorrected PFDR Coefficients Puncorrected PFDR\nRight Right\nRegion Left Right Left Right Left Right Left Right Left Left\nBanks of\nsuperior sulcus 1.45 0.85 0.05 0.24 0.39 0.63 -1.01 -0.52 0.25 0.48 0.99 0.99\nCaudal anterior cingulate -0.36 0.3 0.56 0.65 0.8 0.82 0.97 0.84 0.22 0.26 0.99 0.99\nCaudal middle frontal 0.46 1.67 0.53 0.04 0.8 0.39 0.32 0.81 0.70 0.34 0.99 0.99\nCuneus 0.07 0.76 0.91 0.27 0.96 0.64 -0.63 -0.89 0.47 0.33 0.99 0.99\nEntorhinal 0.68 -0.1 0.16 0.82 0.56 0.92 0.24 -0.18 0.63 0.72 0.99 0.99\nFrontal pole 0.22 0.02 0.55 0.96 0.8 0.96 0.29 0.31 0.50 0.48 0.99 0.99\nFusiform 1.27 1.06 0.14 0.21 0.56 0.6 -0.07 -0.45 0.93 0.60 0.99 0.99\nInferior parietal 1.71 1.98 0.03 0.02 0.39 0.39 -0.96 -0.34 0.28 0.70 0.99 0.99\nInferior temporal 0.75 0.82 0.28 0.28 0.64 0.64 -0.01 -0.12 0.99 0.87 0.99 0.99\nInsula 0.6 0.72 0.53 0.43 0.8 0.79 0.26 -0.33 0.74 0.67 0.99 0.99\nIsthmus cingulate 1.06 -0.14 0.10 0.84 0.54 0.92 0.16 0.6 0.87 0.51 0.99 0.99\nLateral occipital 0.78 0.53 0.30 0.49 0.66 0.79 -0.09 -0.92 0.93 0.32 0.99 0.99\nLateral orbitofrontal 1.56 -0.27 0.05 0.72 0.39 0.87 -0.4 0.12 0.63 0.89 0.99 0.99\nLingual 0.16 0.49 0.82 0.48 0.92 0.79 -0.57 -0.91 0.52 0.33 0.99 0.99\nMedial orbitofrontal 0.29 0.51 0.66 0.47 0.82 0.79 -0.29 -0.48 0.72 0.58 0.99 0.99\nMiddle temporal 1.34 2.18 0.08 0.01 0.5 0.39 -0.8 0.26 0.33 0.75 0.99 0.99\nParacentral 0.73 0.07 0.32 0.93 0.66 0.96 -0.13 0.17 0.87 0.83 0.99 0.99\nParahippocampal 0.8 0.05 0.18 0.95 0.56 0.96 -0.41 -0.98 0.46 0.18 0.99 0.99\nParsopercularis 1.36 0.53 0.16 0.47 0.56 0.79 -0.6 0.15 0.51 0.87 0.99 0.99\nParsorbitalis 0.15 0.52 0.79 0.37 0.92 0.69 -0.38 0.52 0.56 0.48 0.99 0.99\nParstriangularis 1.73 0.09 0.03 0.90 0.39 0.96 0.49 0.14 0.56 0.87 0.99 0.99\nPericalcarine -0.15 -0.06 0.81 0.93 0.92 0.96 -0.84 -0.91 0.35 0.27 0.99 0.99\nPostcentral 0.38 0.87 0.65 0.32 0.82 0.66 -0.03 -1.13 0.97 0.24 0.99 0.99\nPosterior cingulate 1.63 1.16 0.06 0.20 0.44 0.59 0.88 0.34 0.35 0.70 0.99 0.99\nPrecentral 0.49 0.34 0.57 0.65 0.8 0.82 -0.16 0.63 0.85 0.46 0.99 0.99\nPrecuneus 1.15 1.4 0.17 0.10 0.56 0.54 -1.22 -1.22 0.23 0.23 0.99 0.99\nRostral anterior cingulate 2.34 0.58 0.00 0.34 0.29 0.69 -0.82 0.14 0.30 0.84 0.99 0.99\nRostral middle frontal 1.26 0.89 0.15 0.23 0.56 0.62 0.01 -0.1 0.99 0.92 0.99 0.99\nSuperior frontal 0.41 1.2 0.57 0.12 0.8 0.56 -0.24 0.17 0.75 0.83 0.99 0.99\nSuperior parietal 0.57 1.04 0.47 0.18 0.79 0.56 -1.45 -2.1 0.17 0.05 0.99 0.99\nSuperior temporal 0.89 1.05 0.27 0.15 0.64 0.56 0.21 -0.31 0.81 0.69 0.99 0.99\nSupramarginal 0.66 1.4 0.35 0.07 0.69 0.46 -0.05 -0.49 0.95 0.57 0.99 0.99\nTemporal pole 0.12 0.21 0.77 0.60 0.92 0.81 0.38 0.64 0.34 0.15 0.99 0.99\nTransverse temporal -0.27 0.36 0.62 0.51 0.82 0.8 -0.9 -0.09 0.21 0.89 0.99 0.99\n• Please revise the claims on lines 206-208. The model’s poor performance on out-of-distribution\ndata generally does not suggest any specificity to VNS, especially given that the model was pre\ntrained on healthy non-VNS patients. This likely reflects domain shift or lack of\nrepresentativeness in training data, and less likely biological specificity.\nWe thank the reviewer for raising this important point. To clarify, while the encoder-decoder was\npretrained on healthy subjects, the support vector machine (SVM) classifier itself was trained\nexclusively on the VNS cohort. As such, the features driving classification were deviations from\nhealthy anatomy that were specifically associated with VNS outcomes. The poor generalization\nto the temporal lobectomy (TLE) cohort therefore suggests that the latent features of T1w images\nthat predict VNS response do not directly translate to prediction of post-surgical seizure\noutcomes. We have revised the text to avoid implying biological specificity to VNS, instead\nemphasizing that the learned decision boundary is informative for VNS outcomes but not for\nresection outcomes. This has been made clearer in the text:\n“To evaluate the specificity of the SVM to predicting VNS outcomes, we tested its\nperformance on an independent cohort of children who underwent temporal lobectomy for\ntemporal lobe epilepsy39 (n = 102; mean age [SD] = 12.8 [3.9] years; 69.6% seizure freedom).\nThe SVM model did not generalize to this cohort (AUC = 0.57, p = 0.12; Fig. 3g), suggesting\nthat its decision boundary is not applicable to this population.”\n• Throughout the manuscript please be explicit about dimensions and sizes where possible. For\nexample, lines 194-195: “spatial dimension” could be made clearer with subscripts x,y,z.\nWe have specified the x, y, and z dimensions in lines 194-195, and applied the same throughout\nthe manuscript.\n“To evaluate the predictive utility of the latent space for response following VNS, we\ntrained a linear SVM classifier (Fig. 2e) using quantized latent vectors from the frozen VQ-VAE\nencoder. To reduce dimensionality, each four-dimensional latent tensor (x, y, and z spatial\ndimensions × feature channels) was spatially averaged across the x, y, and z axes, yielding a\nsingle 32-dimensional21 feature vector per subject. Each element of this vector therefore\ncaptures the global magnitude of the respective latent feature learned from the underlying\nT1w image.”\n• Explain the choice and dimensionality of the 32-dimension feature space on line 196.”\nThe 32-dimensional feature space corresponds to the latent embeddings produced by the VQ-\nVAE model, which were then aggregated across space for each subject. This dimensionality was\nchosen based on the VQ-VAE implementation in Tudosiu et. al, 2024, where this latent size was\nshown to balance reconstruction fidelity and representational compactness. In our case, the 32-\ndimensional representation provided strong reconstruction accuracy and yielded stable features\nfor downstream classification, so we retained this configuration. We have revised the Methods to\nclarify the chosen dimensionality of the feature space by citing the original work it was based on\nas well as what the dimensions represent. In addition, we have included the dimensions in Figure\n2.\n“To evaluate the predictive utility of the latent space for response following VNS, we\ntrained a linear SVM classifier (Fig. 2e) using quantized latent vectors from the frozen VQ-VAE\nencoder. To reduce dimensionality, each four-dimensional latent tensor (x, y, and z spatial\ndimensions × feature channels) was spatially averaged across the x, y, and z axes, yielding a\nsingle 32-dimensional21 feature vector per subject. Each element of this vector therefore\ncaptures the global magnitude of the respective latent feature learned from the underlying\nT1w image.”\nFigure 2: Training of the VQ-VNS Model. (a) Preprocessing pipeline of T1-w MRI for pretraining of the\nVQ-VNS Quantized Autoencoder Model for T1w reconstruction. Each raw T1w MRI image undergoes\nhybrid brain extraction (SynthStrip and SynthSeg, reorientation to standard RAS orientation, 1 mm\nisotropic resampling, bias field correction, and rigid registration to the MNI 152 template followed by\ncropping to brain margins. (b)-(c) VQ-VAE architecture for T1w MRI reconstruction, comprising a\nconvolutional encoder, vector quantization bottleneck, and convolutional decoder. (d) Composition of the\ndatasets used for VQ-VAE reconstruction training, which included healthy subjects and persons with\nepilepsy. (e) Schematic of the classification pipeline: the spatial mean of the VQ-VAE quantized latent is\ncomputed, feature selected, and used as input to a linear SVM for VNS response prediction.\nFigures\n• The manuscript would greatly benefit from a figure illustrating the full model development\nprocess: from pretraining, to feature extraction, to SVM training and prediction. Each step should\nalso include dimensionality of data flow, especially input, output, codebook, and latent volumes.\nWe have revised Figure 2 to demonstrate the full model pipeline, from the input image to the\noutput classification.\nFigure 2: Training of the VQ-VNS Model. (a) Preprocessing pipeline of T1-w MRI for pretraining of the\nVQ-VNS Quantized Autoencoder Model for T1w reconstruction. Each raw T1w MRI image undergoes\nhybrid brain extraction (SynthStrip and SynthSeg, reorientation to standard RAS orientation, 1 mm\nisotropic resampling, bias field correction, and rigid registration to the MNI 152 template followed by\ncropping to brain margins. (b)-(c) VQ-VAE architecture for T1w MRI reconstruction, comprising a\nconvolutional encoder, vector quantization bottleneck, and convolutional decoder. (d) Composition of the\ndatasets used for VQ-VAE reconstruction training, which included healthy subjects and persons with\nepilepsy. (e) Schematic of the classification pipeline: the spatial mean of the VQ-VAE quantized latent is\ncomputed, feature selected, and used as input to a linear SVM for VNS response prediction.\n• Figure 1g and Figure 3b both have labels for “VNS” but demonstrate modeling results for two\ndifferent models. Please use better labels in figure legends to better distinguish between different\ntypes of models used.\nWe thank the reviewer for pointing out this potential source of confusion. To improve clarity, we\nhave revised the label in Figure 1g to indicate chance performance and clinical model\nperformance. In Figure 3b, we have included a subtitle to indicate that we are demonstrating\nperformance of the VQ-VNS classifier. We have retained the labels to show the model’s\nperformance on patients who receive VNS vs. TLE patients, as the intent of this panel is to\ncompare across these groups.\n3\n1\n• Figure 2g borders on visual exaggeration. Please plot x-axis from zero and have a zoomed-in\nsubpanel or use a different plot type. Please also clarify how data from this figure were acquired.\nFor example if this is an aggregate measure across subjects, please show variations across\ngenerations for each region.\nIn this revised manuscript, Figure 2g corresponds to Figure 3e. We have revised Figure 3e to\ndisplay the differences in a zero-centered plot to enhance interpretability. The data in this panel\nare derived directly from the synthetic, archetypal responder and non-responder shown in panel\nd. The purpose of panel e is to highlight the volumetric differences between these archetypal\nreconstructions, which were generated by averaging the responder and non-responder groups,\nrespectively, as inputs to the VQ-VAE model.\nFigure 3: Performance of the VQ-VNS Quantized Autoencoder Model for T1w Reconstruction and\nClassifier for Predicting VNS Outcomes. (d) Synthetic images of the archetypal VNS responder and non-\nresponder reconstructed from the latent space using the VQ-VAE decoder. Differences are subtle and not\nvisually apparent. (e) The synthetic archetypal non-responder showed reduced grey matter, white matter,\nand subcortical grey matter volumes, and ventricular enlargement relative to the synthetic archetypal\nresponder.\nModeling\n• Explain the choice of using an SVM as the classifier for the embeddings rather than an\nensemble of models or other alternatives (random forest etc.). I strongly suggest comparing\ndifferent modeling strategies, where applicable.\nOur choice of SVM as the classifier for the learned embeddings was guided by both theoretical\nconsiderations and empirical evaluation. Among several standard classifiers tested, including\nlogistic regression, random forest, and k-nearest neighbors (kNN), SVM consistently provided\nthe best balance between model complexity, robustness to overfitting, and predictive accuracy on\nheld-out data. Logistic regression tended to underfit, random forests showed overfitting\ntendencies, and kNN performed less favorably overall.\nThese comparisons are now presented in Supplementary Figure 6, which confirms that the SVM\nachieved the highest classification accuracy for the learned latent representations. Together,\nthese results support its selection as the final classifier in the VQ-VNS framework.\nSupplementary Figure 6. Performance of support vector machine, logistic regression, k nearest neighbours (KNN),\nand random forest classifier in predicting response to VNS from latent feature representations of T1w MRI.\n• Please include a full suite of model assessment metrics where applicable (Accuracy, F1,\nPrecision, Recall, etc.)\nWe have included a full suite of performance metrics for each model assessed in this work –\nincluding accuracy, precision, recall, F1, and specificity – in Supplementary Table 9.\nSupplementary Table 8. Model Assessment Metrics\nModel Accuracy Precision Recall F1 Specificity\nClinical Data 0.53 0.53 0.51 0.52 0.54\nTest set 0.72 0.70 0.74 0.72 0.70\n50% 0.66 0.67 0.67 0.67 0.65\nBinarization\nTLE 0.56 0.69 0.62 0.65 0.44\n• For each subset of data (test, validation, test) please explicitly provide the distribution of\nresponders and non-responders and be explicit about data imbalances and highlight mitigation\nstrategies.\nTo minimize class imbalance, data were stratified at the time of splitting so that each subset was\nas close to 50/50 as possible. We have clarified this in the Methods Section “Training the SVM\nClassifier to Predict VNS Response.”\n“To assess the stability and significance of the learned features, we first trained 1000 SVM\nclassifiers with class balancing and a high regularization parameter (C=5000) on randomly\npartitioned subsets of the VNS cohort. To assess generalizability, a final SVM was trained with a\ntraining, validation set, and held-out test set. To minimize data imbalance, participants were\nsplit into training (n = 199, responders = 100, non-responders = 99), validation (n = 25,\nresponders = 13, non-responders = 12), and held-out test sets (n = 39, responders = 19,\nnon-responders = 20) using stratified sampling to preserve approximately equal proportions of\nresponders and non-responders in each subset.”\n• Add reconstruction scores to the images shown\nIn the revised manuscript, we have added MS-SSIM reconstruction scores to all displayed\nimages to quantitatively illustrate reconstruction quality.\nFigure 3: Performance of the VQ-VNS Quantized Autoencoder Model for T1w Reconstruction and\nClassifier for Predicting VNS Outcomes. (a) Representative reconstruction of a healthy subject from the\nheld-out test set of the VQ-VAE pretraining dataset. (b) Correlations between true and reconstruction\nvolumes across four tissue classes in the pretraining test set, demonstrating high reconstruction fidelity.\n(c) Representative reconstruction of a subject with a previous resection from the CONNECTiVOS cohort.\n(d) Synthetic images of the archetypal VNS responder and non-responder reconstructed from the latent\nspace using the VQ-VAE decoder. Differences are subtle and not visually apparent. (e) The synthetic\narchetypal non-responder showed reduced grey matter, white matter, and subcortical grey matter\nvolumes, and ventricular enlargement relative to the synthetic archetypal responder. (f) 1000 models\ntrained using hold-out cross validation perform significantly better than chance. (g) Receiver operating\ncharacteristic (ROC) curve showing the predictive performance of the SVM classifier on the held-out test\nset (AUC = 0.73, p = 0.007), which was specific to VNS and did not generalize to outcomes following\ntemporal lobectomy (AUC = 0.54, p > 0.99).\nSupplementary Figure 3: Reconstruction Performance of VQ-VAE Model. Representative\ninput images and corresponding reconstructions for ten randomly selected subjects from the test\nset (n = 744) of the pretraining dataset (n = 7443). Each row shows axial, sagittal, and coronal\nslices for the input image (left) and its reconstruction (right). The model preserves fine\nanatomical structure, captures intensity variation across individuals, and reduces motion-related\nartifacts visible in some original scans.\nSupplementary Figure 4: Reconstruction Performance Across Abnormal Imaging Data (e)\nVQ-VAE reconstructions and corresponding input images, demonstrating marginally accurate\nmodelling of ventricular enlargement and resection cavities, but reduced fidelity in replicating\nabnormal cortical folding, particularly in lissencephaly.\n• In supplemental figure 3, please double check that the cases shown in a-d do in fact correspond\nto the inputs in panel e\nThank you for this comment. We have revised Supplementary Figure 4 (previously 3) so that the\nexamples shown in a-d correspond to the inputs in panel e.\nSupplementary Figure 4: Reconstruction Performance Across Abnormal Imaging Data (a)\nTemporal lobectomy with prior surgical changes. (b) Subject with lissencephaly, characterized\nby absent sulcation. (c) Polymicrogyria and heterotopias. (d) Extensive periventricular\nleukomalacia. (e) VQ-VAE reconstructions and corresponding input images, demonstrating\nmarginally accurate modelling of ventricular enlargement and resection cavities, but reduced\nfidelity in replicating abnormal cortical folding, particularly in lissencephaly.\n• Ensure all mathematical equations are fully defined and explained with all parameters clearly\nlabeled.\nIn the revised manuscript, we have ensured that all mathematical equations are fully defined,\nwith each parameter explicitly described.\n“During training, we employed a composite loss function incorporating voxel-wise accuracy,\nspectral alignment, perceptual similarity, encoder-quantizer alignment (commitment), and\nadversarial objectives.3 Voxel-wise error (ℒ ) was quantified with mean absolute error (MAE)\nvox\nbetween the original 𝑋 and reconstructed 𝑋̂ images.\nℒ = 𝑀𝐴𝐸(𝑋,𝑋̂) (1)\nvox\nSpectral alignment was enforced using a frequency loss function (ℒ ) computed as the mean\nfreq\nsquared error (MSE) between the fast Fourier transforms9 of the original and reconstructed\nimages:\nℒ = 𝑀𝑆𝐸(𝐹𝐹𝑇(𝑋),𝐹𝐹𝑇(𝑋̂)) (2)\nfreq\nPerceptual similarity was assessed using MONAI's perceptual loss function (ℒ ), which\npercep\nleveraged the Learned Perceptual Image Patch Similarity (LPIPS)10 metric with AlexNet11 as the\nfeature extractor. Given that the 2D AlexNet model cannot directly process 3D volumes, we\ncomputed the LPIPS metric on individual slices extracted from the axial (ax), coronal (cor), and\nsagittal (sag) views, subsequently summing these per-slice scores as in:\nℒ = 𝐿𝑃𝐼𝑃𝑆(𝑋 ,𝑋̂ )+ 𝐿𝑃𝐼𝑃𝑆(𝑋 ,𝑋̂ )+ 𝐿𝑃𝐼𝑃𝑆(𝑋 ,𝑋̂ ) (3)\npercep ax ax cor cor sag sag\nCommitment loss (ℒ ) in the VQ-VAE encouraged the encoder outputs to commit to\ncommit\ndiscrete codebook entries. Since the codebook was updated using an exponential moving average\n(EMA) approach,16 gradients should not flow directly into the codebook vectors. Accordingly, we\napplied a stop-gradient operator to the codebook entry:\n2\nℒ = 𝽽‖𝑧 (𝑋)−𝑠𝑔(𝑒)‖ (4)\ncommit e 2\nWhere 𝽽 is a scaling constant, 𝑧 (𝑋) is the encoder’s output, 𝑒 is the closest codebook vector,\ne\nand 𝑠𝑔(∙) denotes the stop gradient operator.\nAdversarial loss (ℒ ) guided the generator to produce outputs that are indistinguishable from\nadv\nreal images by a patch-based discriminator 𝐷,12 thereby retaining local consistency and fine\nstructural detail. The loss employed a least-squares criterion with mean reduction, targeting\ndiscriminator outputs of 1 for generated samples. No activation was applied to the discriminator\noutputs prior to loss computation:\nℒ = (𝐷(𝑋̂)−1) 2 (5)\nadv\nTogether the expression representing the total loss term was,\nℒ = ℒ +ℒ +ℒ +ℒ +ℒ (6)\ntotal vox freq percep commit adv\nTraining stability across loss terms with differing magnitudes was maintained through the\napplication of gradient normalization, a method inspired by GradNorm.13 This involved\nrescaling gradients according to their norms, as detailed below:\n𝑤𝑔\n𝑔̃ = (7)\n||𝑔||+ 𝜀\nwhere 𝑔̃ represents the normalized gradient, 𝑔 is the original gradient, 𝑤 is a scaling factor, and\n𝜀 is a small constant to prevent numerical instability. This approach rescaled all gradients based\non their norms, effectively balancing their contribution to the training update and preventing\ndominance by large gradients or attenuation of smaller ones. The scaling factors were\nempirically selected.\nThe discriminator was trained using mean squared error (MSE) loss between its predicted\nscores for real and generated images, denoted as ℒ :\ndis\n1\nℒ = [𝑀𝑆𝐸(𝐷(𝑋))+𝑀𝑆𝐸(𝐷(𝑋̂))] (8)\ndis 2\nIn addition to gradient normalization, we applied LeCam regularization14 to further stabilize\nadversarial training. This method limits how much the discriminator’s outputs can fluctuate,\nhelping prevent mode collapse and improving generalization.\nLet 𝐷(𝑋) and 𝐷(𝑋̂) denote the discriminator logits for real and generated images, respectively.\nWe defined LeCam anchors 𝐴 and 𝐴 as exponential moving averages of real and fake logits:\nr f\n𝐴 ← 𝽽𝐴 +(1−𝽽)𝐸[𝐷(𝑋)], 𝐴 ← 𝽽𝐴 +(1−𝽽)𝐸[𝐷(𝑋̂)] (9)\nr r f f\nwhere 𝽽 is a momentum parameter and E represents the average of 𝐷(𝑋̂). The LeCam loss\n(ℒ ) penalized deviations of real logits from the fake anchor and fake logits from the real\nLeCam\nanchor:\nℒ = 𝐸[( 𝐷(𝑋)−𝐴 ) 2 ]+𝐸[( 𝐷(𝑋̂)−𝐴 ) 2 ] (10)\nLeCam f r\nThis regularization discourages abrupt changes in discriminator outputs over training and\nprevents overconfident separation of real and generated inputs. The final discriminator loss\ncombined the adversarial term with the LeCam penalty, scaled by a tunable coefficient:\nℒ = ℒ +𝜆 ∙ℒ (11)\nadv dis LeCam LeCam\n• In the UMAP embeddings the manifold shows some separation of data. Were those data points\ninspected further to determine which factors could have contributed to their separation? If so,\nplease explain or suggest possible separation factors.\nWe thank the reviewer for this insightful observation. To examine whether the possible manifold\nstructure in the UMAP embeddings reflected underlying biological or methodological factors, we\nperformed an additional principal component analysis (PCA) on the same dataset for greater\nexplainability. The top ten principal components collectively explained approximately 60% of\nthe total variance, with no single component dominating (the first principal component explained\nonly approximately 12% of the variance) (Supplementary Fig. 2). Importantly, visualization of\nthe first three components showed no discernible clustering by treatment outcome or age, and a\nslight separation by acquisition site (Supplementary Fig. 2). These findings suggest that the\napparent structure in the UMAP projection likely reflects the intrinsic variance of the dataset,\nwith a small contribution from acquisition site, but no systematic separation associated with\ntreatment outcome. We have included a summary of this analysis in the revised manuscript\n(Results, p. 5; Supplementary Fig. 2).\n“No clinical features were significantly associated with VNS outcome after false discovery rate\ncorrection (Supplementary Table 1). Nonlinear embedding of the clinical data using uniform\nmanifold approximation and projection31 (UMAP) also failed to separate responders and non-\nresponders (Fig. 1e; s = –2.7 × 10-4, p = 0.71), which was also corroborated by principal\ncomponent analysis which revealed low explained variance (<12% for PC1) and no systematic\nresponse associated differences (Supplementary Fig. 2). Labelling the UMAP embeddings by\ntreatment site revealed visually apparent but objectively weak (s = –0.04, p > 0.99) site-specific\ngroupings (Fig. 1f), indicative of subtle inter-site variability in patient selection criteria without\nstrong site-level effects.”\nSupplementary Figure 2: Principal component analysis of clinical dataset. (A) Boxplots of first\n10 principal component scores grouped by treatment outcome (NR = non-responder, R =\nresponder). None of the top 10 components demonstrated systematic separation by outcome. (B)\nScree plot of the percentage of variance explained by each component. (C) The first 10\ncomponents collectively explained ~60% of the total variance, with PC1 contributing ~12%.\nSupplementary Table 2. Loadings of demographic and clinical variables on the first ten\nprincipal components.\nPC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10\nAge 0.18 0.20 0.06 -0.19 -0.17 -0.12 -0.02 0.04 0.02 -0.03\nNumber of ASMs 0.15 0.05 0.14 0.02 0.10 0.03 0.16 0.08 0.22 -0.12\nDuration of epilepsy (yrs) 0.15 0.10 0.09 -0.24 0.04 -0.24 -0.02 0.16 0.00 0.01\nMale 0.03 -0.03 0.09 -0.07 0.22 -0.21 -0.25 -0.47 0.13 -0.23\nGenetic etiology 0.04 -0.17 0.00 0.01 0.13 0.18 0.35 -0.26 -0.04 0.24\nDevelopmental delay 0.08 -0.11 0.10 0.15 0.33 0.13 0.10 0.14 0.14 0.11\nFamily history of epilepsy -0.10 -0.09 0.09 -0.19 -0.18 -0.26 0.21 -0.11 0.35 0.24\nHistory of infantile spasms 0.11 -0.13 0.04 0.25 0.28 -0.21 0.00 0.13 0.03 -0.06\nCongenital etiology 0.08 0.04 -0.03 0.08 0.20 -0.26 0.18 -0.07 -0.34 -0.13\nAcquired etiology 0.06 0.12 0.04 0.08 0.00 0.12 -0.35 0.24 0.31 0.07\nPrevious surgery 0.05 0.17 -0.13 0.12 0.05 -0.32 0.00 0.05 -0.05 0.37\nGeneralized onset 0.12 -0.24 0.28 0.05 -0.20 -0.14 0.06 0.10 -0.02 -0.10\nNon-generalized onset -0.07 0.32 -0.14 -0.14 0.28 0.06 0.04 -0.04 0.13 0.06\nRight-lateralized EEG 0.04 0.18 0.00 0.28 -0.13 0.08 0.29 -0.01 0.13 -0.13\nLeft lateralized EEG 0.06 0.21 -0.03 0.27 -0.14 0.02 0.26 0.03 0.16 -0.14\nBilateral EEG anomalies 0.32 -0.07 0.15 0.07 0.07 0.11 -0.03 -0.03 0.06 0.15\nLennox-Gastaut syndrome 0.11 -0.01 0.04 0.17 0.08 -0.11 0.03 -0.07 0.29 0.08\nSCN1A mutation present -0.03 -0.09 0.00 -0.03 0.01 0.19 0.06 -0.05 -0.16 0.24\nTuberous sclerosis 0.06 -0.02 0.00 0.01 0.13 -0.12 0.13 0.05 -0.08 0.07\nFocal cortical dysplasia -0.01 0.09 -0.08 0.02 -0.02 -0.07 0.17 0.08 -0.03 -0.15\n• Methodological details tend to be lacking or missing for secondary comparisons. For example,\nin examining subjects with between 20%-80% seizure reduction, evaluating reconstruction with\nabnormal imaging etc. Please provide additional details and organize carefully.\nFor the secondary comparisons, the methodology was identical to that used for the primary\nanalyses. We have added a clarifying statement in the text to make this explicit.\n“To further assess robustness, we tested the VQ-VNS model’s performance in subjects with\nintermediate outcomes seizure reductions between 20% and 80%, without further model\ntraining. For this analysis response was binarized at the conventional >50% threshold.\nSimilarly, model performance was evaluated using AUC and confusion matrices.”\n“To evaluate the fidelity of VQ-VAE reconstructions, we computed three complementary image-\nlevel metrics – mean squared error (MSE), multiscale structural similarity (MS-SSIM),15 and\nFréchet Inception Distance (FID)16) between the original and the reconstructed 3D T1-weighted\nimages. MSE and MS-SSIM were computed using MONAI (MSEMetric, MultiScaleSSIMMetric),\nand FID was calculated using features extracted from a pretrained ResNet-50 model,17\nfollowing global average pooling across spatial dimensions. To assess the VQ-VAE model’s\nrobustness to anatomical abnormalities, we performed a secondary evaluation on 106 subjects\nwith abnormal imaging findings. Reconstruction fidelity was evaluated using the same three\nmetrics – MSE, MS-SSIM, and FID.\n• For lines 527-530, please reword for improved clarity.\nThank you for this suggestion. We have reworded lines 527-530 for improved clarity.\n“To further assess robustness, we tested the VQ-VNS model’s performance in subjects with\nintermediate outcomes seizure reductions between 20% and 80%, without further model\ntraining. For this analysis response was binarized at the conventional >50% threshold.\nSimilarly, model performance was evaluated using AUC and confusion matrices.”\n• Please clarify/detail assessment steps used for subjects between 20%-80% seizure reduction.\nFor the secondary comparisons, the methodology was identical to that used for the primary\nanalyses, and no new model was trained. Instead, the performance of the same model was tested\non the cohort of participants who experience a 20-80% reduction in seizures. We have added a\nclarifying statement in the text to make this explicit.\n“Similarly, model performance was evaluated using AUC and confusion matrices.”\n• In evaluating the cortical thickness estimation, explain the choice behind using t-SNE vs other\ntechniques. What are its benefits and shortcomings?\nWe used t-distributed stochastic neighbour embedding (t-SNE) as a dimensionality reduction\ntechnique to visualize cortical thickness patterns in a two-dimensional space24 . t-SNE is widely\nused for examining high-dimensional data because it preserves local structure and highlights\nneighbourhood relationships between subjects. This makes it ideal for visualizing subtle\ngrouping patterns in complex neuroimaging datasets. We recognize, however, that t-SNE also\nhas important limitations. The resulting embeddings can be sensitive to hyperparameters and\ninitializations, and global distances are not well preserved, which may complicate interpretation.\nNonetheless, for our goal of exploring relative clustering and qualitative separation of cortical\nthickness patterns across subjects, t-SNE provides an effective and widely accepted\ndimensionality-reduction tool that enabled effective visualization and quantitative analysis. We\nhave included a brief discussion in the revised manuscript.\n“To evaluate potential site-related batch effects in the cortical thickness data, we applied t-\ndistributed stochastic neighbor embedding24 (t-SNE) to all cortical vertex values and estimated\ntotal intracranial volumes (eTIV), projecting them into two dimensions. To preserve local\nneighbourhood structure in high-dimensional data, t-SNE was selected. Although t-SNE does\nnot preserve global distances well, it offers an established and effective visualization tool for\nour goal of exploring relative clustering and assessing potential batch effects. Embeddings\nwere visualized by site (Supplementary Fig. 9a), outcome, and age (Supplementary Fig. 9b–c).\nAge showed the strongest correlation with t-SNE component 1 (Supplementary Fig. 9d). As no\nconsistent clustering by site was observed, no further harmonization was applied.”\n• Explain the choice between permuting model outcomes versus other model comparison\ntechniques.\nThe use of permutation testing was not intended as an alternative, but rather as a complementary\nevaluation to provide an empirical significant estimate of the observed AUC. The purpose of the\npermutation analysis is to demonstrate that the model performs significantly better than chance,\nwhich reflects the baseline level of current clinical decision-making. By permuting the labels,\nwe demonstrate that mismatched labels yield a null distribution of performance, confirming that the\nobserved accuracy arises from the model learning a mapping between T1 features and outcome\nrather than chance. We have added clarification in the text.\nIn the revised manuscript, we now also report a full suite of conventional metrics (accuracy,\nprecision, recall, F1, and specificity) in addition to permutation testing (permuting the test/train\nsplit with refitting and retesting). Our goal is to provide a comprehensive assessment of model\nperformance, using both conventional measures and permutation-based evaluation for\nrobustness.\n“Robustness was assessed by permutation testing with 1,000 label shuffles to generate a null\ndistribution of AUC scores, enabling empirical significance estimation. This procedure tests\nwhether the VQ-VNS model’s observed performance exceeds chance expectations, thereby\nconfirming that predictions reflect a learned mapping between input features and outcome.”\n• In order to improve potential clinical applicability, consider framing predictions as probability\nscores of response (e.g. P(Response|Input)) rather than binary outputs.\nWe thank the reviewer for this insightful suggestion. We agree that probabilistic predictions\nwould improve clinical applicability. In our current implementation using a linear SVM,\nhowever, the decision boundary is very fine, resulting in most samples lying close to the\nseparating hyperplane. While this behavior is appropriate for the model’s internal decision-\nmaking, the distance from the decision boundary in our implementation does not yield\nconfidence scores with a wide enough range to support clinical decision making, given the fine-\ngrained separation in this dataset.\nIn our revised manuscript, we have presented results for a probabilistic Logistic Regression\nmodel (Supplementary Figure 7), but future work exploring probabilistic classifiers or\napproaches explicitly designed to calibrate output probabilities will be required to generate\nprobability scores that are better suited to support clinician decision making. We have added this\nas a future direction in the Discussion.\n“Finally, future work should incorporate probabilistic classifiers to yield confidence\nestimates that better support clinician decision-making.”\nSupplementary Figure 6. Performance of logistic regression, k nearest neighbours (KNN),\nand random forest classifier in predicting response to VNS from latent feature representations of\nT1w MRI.\nData and Code\n• I strongly suggest porting the model weights over to a more stable, permanent repository rather\nthan its current form, which is inaccessible. Some domain appropriate archives include\nOpenNeuro, Zenodo, etc.\n• This work would greatly benefit from a demonstrative Jupyter notebook. This could potentially\nbe shown with anonymized or synthetically generated responder/non-responder T1w MRI scans,\nwhich are then used to demonstrate feature extraction, visualization, and then prediction with the\ntrained classifiers.\nWe thank the reviewer for these suggestions. The model weights have now been uploaded to a\nstable and publicly accessible repository on Hugging Face. We have also added a demonstrative\nJupyter notebook to the project’s GitHub repository, which showcases the full inference pipeline\nusing a synthetically generated T1-weighted MRI.\nHugging face model: https://huggingface.co/hsuresh/vqvns\nGithub repository: https://github.com/gmilab/VQVNS\nMinor Concerns\nFigures\n• Please unionize all the font types across figures and sub panels\nIn the revised manuscript, we have unionized fonts across figures.\n• Double check that figures conform to the journal's guidelines (i.e. font sizes etc.)\nThank you for this comment. We have referred to the journal’s guidelines.\n• For figure 1d, discretize colors, space them out more and add the rounded percentages to the\ncorresponding block for clarity. For example, if “Development Delay” accounts for 90% of\nsubjects please label the bar with the corresponding percentage. Under each label add the type of\noutcome (i.e. binary, continuous, etc.)\n• Figure 1d please make axes square for clearer visual comparison\nThank you for these suggestions, we have incorporated them into the revised Figure 1d. All\noutcomes included in the figure are binary.\nFigure 1: Predictive Performance of Clinical Data (d) Clinical data, stratified by VNS response (p >\nFDR\n0.05 for difference between groups across all variables, Supplementary Table 1).\n• Figure 1c, please add the sample size for each subgroup (responders/non-responders). Can be\nadded inside the plot under each group or centered under the corresponding label.\nIn Figure 1C, 348 responders and 346 non-responders were included. These numbers have been\nadded to the figure.\nResults\n• Add min and max age ranges where applicable\nThe minimum age was 0.8 years, and the maximum age was 24.2 years. We have included this in\nthe manuscript.\n“The mean [SD; range] age was 10.6 [4.7; 0.8-24.2] years, 54% male and included 51%\nresponders (Fig. 1b-c).”\nReviewer #2 (Remarks on code availability):\n* I was not able to set up and run the code. I strongly suggest presenting verification of this work\nwithin a Jupyter notebook and plotting outcomes as interactive visualizations. Leverage Docker\nif necessary.\n* Model weights were not easily accessible, please port to an accessible platform similar to\nZenodo, OpenNeuro or others.\nWe appreciate the reviewer’s feedback and have addressed these concerns. The model weights\nhave now been ported to Hugging Face. In addition, we have created a Jupyter notebook\ndemonstration that reproduces the full inference workflow using a synthetic T1-weighted MRI.\nHugging face model: https://huggingface.co/hsuresh/vqvns\nGithub repository: https://github.com/gmilab/VQVNS\nReviewer #3 (Remarks to the Author):\nClinical\nResearchers developed VQ-VNS, a deep learning model that predicts the effectiveness of vagus\nnerve stimulation (VNS) for children with epilepsy using routine preoperative MRI scans.\nTraditional clinical data could not predict treatment response, but VQ-VNS achieved significant\npredictive accuracy and highlighted serotonin-rich brain regions and connectivity disruptions in\nnon-responders. These results are clinically relevant and highly translational, offering a\nbiologically interpretable tool to guide decisions and reduce unnecessary invasive procedures.\nMethods\nPatient grouping: Participants were divided into >80% seizure reduction vs. <20%. What about\npatients between these cut-offs? Excluding ~350 patients (694/1046 included) seems a large data\nloss. Since you collected seizure frequency reduction as a continuous variable, would correlating\nwith this continuous measure be more informative?\nWe thank the reviewer for raising this important point. We are keenly aware that patient and\nparent-reported seizure frequency reductions are less reliable in intermediate ranges, introducing\nuncertainty and potential misclassification in outcome labeling.20 To minimize this bias – and to\nrepresent the bimodal distribution of the data – we restricted the primary analysis to participants\nat the extremes of response (>80% or <20% seizure reduction), where outcome classification is\nmost reliable. While this approach necessarily reduced the sample size, it ensured greater label\nvalidity and outcome confidence for model training and evaluation.\nTo assess the robustness of this decision, we conducted a secondary analysis including\nparticipants with intermediate seizure reductions (20-80%). Even under this broader inclusion\ncriterion, the model maintained above-chance (albeit reduced) prediction accuracy, supporting\nthe generalizability of the learned representations despite increased label uncertainty.\nClinical outcomes: In the analysis section, you note that including patients with ambiguous\noutcomes did not improve SVM performance. This is an important observation — why not\ninclude them in the main analysis if results are unchanged? Please clarify the added value of\nexcluding them.\nAs noted, including participants with intermediate or ambiguous outcomes did not substantially\nchange overall SVM performance. However, we opted to retain the stricter outcome definition in\nthe primary analysis to ensure that the model was trained on the most reliable and biologically\nmeaningful labels. Incorporating uncertain outcomes can artificially stabilize accuracy metrics\nby introducing label noise that the model effectively “averages out,” but this does not guarantee\nthat the learned representations are clinically valid.\nBy training on well-defined extremes (>80% and <20% seizure reduction), we ensured that the\nlearned embeddings reflect robust neurobiological differences between clear responders and non-\nresponders. The secondary analysis including intermediate cases serves as a complementary\nrobustness check, confirming that model performance remains stable even when label\nuncertainty is increased.\nClinical variables: Have magnet effects and quality of life (QOL) been addressed? Both are well-\nestablished as clinically relevant (Danthine V., Boon P.), and could enrich the dataset. If\navailable, these should be included.\nWe completely agree that these are very important considerations. Unfortunately, these data were\nnot available for the majority of subjects in our cohort. Our study was not sufficiently powered to\nevaluate QoL outcomes in relation to imaging-derived predictors. Moreover, part of the novelty\nof the present work lies in demonstrating that treatment response can be predicted using only a\nstructural T1-weighted MRI, without requiring additional clinical variables. We agree that\nintegrating magnet effects and QoL would enrich future models when larger, multimodal\ndatasets become available.\nSVM classifier: Only 263 patients were included after excluding those with prior surgery or\nsignificant abnormalities. Please describe their clinical characteristics. Were they representative\nof the larger cohort (e.g., intellectual disability, DEEs, epilepsy duration)? The generalizability\nof the model is limited if large patient groups are excluded.\nWe thank the reviewer for this important comment. As expected, this exclusion led to a\nmarginally lower proportion of patients with developmental delay, which frequently occurs with\nstructural abnormalities. Conversely, a greater proportion of generalized-onset epilepsies were\nobserved in the analyzed subset, consistent with the exclusion of focal lesional cases. A smaller\nproportion of the 263 patients included in model training had a development delay and\ngeneralized onset seizures. These differences represent a limitation of our training set, which did\nnot include abnormal brains. Overall, the included patients remained demographically and\nclinically representative of the broader VNS cohort with respect to age, response, and other\nunderlying etiology. These characteristics are now summarized in Supplementary Table 12\nalongside those of the excluded patients.\nWe note in the Discussion that, because patients with abnormal brain anatomy were excluded\n(approximately 10% of the screened cohort), model generalizability to these cases has yet to be\ndetermined.\n“Third, images with significant anomalies or motion artifact were excluded from response\nprediction because the pretraining datasets did not include images with significant anatomical\nabnormalities. This may limit generalizability to grossly abnormal brains, which represented\n<10% of our datasets (45 out of 489 images screened, Supplementary Fig. 4).”\nFigure R2. Demographics of subset used for classifier training, validation, and testing.\nSupplementary Table 12. Demographic and clinical characteristics patients included in model\ntraining (n = 263) compared with full dataset (n = 1046).\nFull population count Subset used for classifier Chi squared\nCharacteristic (proportion) count (proportion) p-value\nGenetic\ncondition 451 (43.1%) 121 (46.0%) 0.391\nDevelopmental\ndelay 744 (71.1%) 150 (57.0%) 0.000\nFamily history 288 (27.5%) 79 (30.0%) 0.414\nInfantile spasm 153 (14.6%) 31 (11.8%) 0.238\nLennox-Gastaut\nsyndrome 105 (10.0%) 32 (12.2%) 0.311\nSCN1A 47 (4.5%) 14 (5.3%) 0.566\nCongenital\ncondition 267 (25.5%) 58 (22.1%) 0.247\nAcquired\ncondition 204 (19.5%) 43 (16.3%) 0.245\nPrevious surgery 136 (13.0%) 29 (11.0%) 0.391\nGeneralized\nonset 555 (53.0%) 163 (62.0%) 0.009\nNon-generalized\nonset 507 (48.4%) 113 (43.0%) 0.113\nRight sided IEDs 325 (31.0%) 75 (28.5%) 0.427\nLeft sided IEDs 377 (36.0%) 90 (34.2%) 0.589\nBilateral IEDs 786 (75.1%) 204 (77.6%) 0.400\nTLE cohort: It is unclear why specificity analyses were limited to a TLE subgroup. Were these\npatients outside the CONNECTIVOS trial? Were they seizure free? Why not include other\nsurgical subgroups (e.g., ETLE)?\nThank you for the opportunity to clarify. The CONNECTIVOS registry is a prospective study\nfocused specifically on VNS therapy; therefore, other surgical cohorts (e.g., those with TLE)\nwere not systematically enrolled. For the purpose of evaluating model specificity, we utilized an\nindependent, previously published, retrospective cohort of patients19 who underwent temporal\nlobectomy because high-quality preoperative imaging and clinical outcome data were available\nfor this group.\nNot all patients in the TLE cohort achieved seizure freedom (~60%), allowing us to evaluate\nmodel behavior across a clinically heterogeneous but well-characterized surgical population. We\nfocused on TLE because it represents the most common and well-defined surgical epilepsy\nsubgroup, enabling a more controlled assessment of whether the learned VNS-related imaging\nrepresentations generalize to a distinct treatment population.\n“To evaluate the specificity of the VQ-VNS model’s performance, we tested its ability to predict\nseizure freedom in a separate cohort of patients with temporal lobe epilepsy26–28 (TLE) who\nunderwent temporal lobectomy (n = 102; mean age [SD] = 12.8 [3.9] years; 69.6% seizure\nfreedom). This retrospective TLE cohort was selected because high-quality preoperative T1-\nweighted imaging and detailed postoperative outcome data were available, enabling an\nassessment of whether the VQ-VNS model, training exclusively on VNS recipients, generalizes\nto an independent surgical population. In this cohort, approximately 70% of patients achieved\nseizure freedom following surgery, which enabled the evaluation of model specificity. T1-\nweighted images were acquired using Philips and GE scanners (Achieva, Genesis Signa, Signa\nHDxt, Signa Excite) with field strengths of 1.5T or 3T. Imaging parameters were: flip angle =\n13.6° (range: 8.0°–30.0°), echo time = 4.3 ms (2.3–4.6 ms), repetition time = 8.7 ms (5.0–25.0\nms), and mean voxel size = 0.9 × 0.9 × 1.4 mm. Abnormal images were excluded in this cohort\nsimilar to the VNS cohort.”\nMethodological consistency: For the VNS analysis, patients with brain lesions were excluded. Is\nthis a methodological requirement, and if so, why would it not apply to the surgical cohort?\nWe thank the reviewer for noting this inconsistency. We excluded patients with brain lesions\nfrom the VNS analysis because there was insufficient representation of specific lesion types to\nallow the model to learn meaningful representations of these more heterogenous and\nanatomically variable abnormalities. As such, it performs optimally on brains that are\nmorphologically closer to normative anatomy.\nTo ensure methodological consistency, we have now applied the same exclusion criteria to the\nsurgical cohort to ensure methodological consistency across analyses. After this additional\nquality control step, the surgical test set comprised 102 subjects. The updated model for TLE did\nnot improve above chance, supporting specificity to VNS (ROC-AUC of 0.57, p = 0.12).\nFigure R3. Confusion matrix and AUC curve of SVM evaluation on TLE cohort.\nWould comparison to a DBS cohort (ANT, CM) be more appropriate than only TLE?\nWe completely agree that comparison to patients with DBS would provide valuable context. We\nalso have a prospective registry of DBS in children, however, at present, our dataset does not\ninclude a sufficient number of DBS cases to permit a meaningful inclusion.\nDefinitions: How was developmental delay defined, and does this overlap with DEE cases?\nDevelopmental delay was reported by each participating site according to their local clinical\nassessment, independent of whether or not the participant met criteria for a developmental\nepileptic encephalopathy (DEE). As such, the designation may reflect a heterogeneous set of\nclinical criteria, and we acknowledge this as a limitation.\nAs the reviewer suggests, there is likely an overlap between cases labeled as “developmental\ndelay” and patients who might also meet criteria for developmental and epileptic encephalopathy\n(n = 180).\nParameters: The discussion states VNS was “optimally programmed.” Please define this — were\nparameters such as rapid cycling or saturation cycling considered (e.g., Norpulse study)?\nWe thank the reviewer for this important point. In our study, VNS programming\nparameters were not standardized. They were left to the discretion of the treating clinicians at\nindividual sites. This variability is a limitation of the study and is discussed in the revised\nmanuscript. We did request stimulation parameter data, though they were only available for\nn=103. Within this cohort, the average output current was 1.8mA and the average duty cycle was\n20%. There was no difference in output current (1.90.7 vs. 1.70.7; t=-1.4; p=0.17; n=103) or\nduty cycle (2110% vs. 1910%; t=-1.3, p=0.2, n=103) between responders and non-responders.\nNo patient received rapid cycling or saturation cycling.\n“This study has several limitations. First, while VNS implantation is highly standardized with\nminimal inter-subject variability, stimulation parameters were not standardized across sites and\nwere adjusted at the discretion of the treating clinician, leading to potential differences in the\ncharge density delivered. However, among participants with available data (N=103), mean\nstimulation parameters (output current 1.8mA; duty cycle 20%) were consistent with optimally\nprogrammed populations.26”\nData splitting: The dataset was split 80/10/10 for training/validation/test. What is the rationale\nfor this choice? Did results vary with different splits (e.g., 50/25/25)? Please provide references\nor justification.\nYes, changing the split would change model performance. The dataset was divided into training,\nvalidation, and test sets using a stratified split that approximated an 75/10/15 distribution, to\nprovide greater power for out-of-sample evaluation on the test set. This strategy maximized the\nsize of the training set while retaining sufficient subjects for validation and robust assessment on\na held-out test cohort, which is commonly used in the context of clinical datasets. 27–29 While\nother splits are possible, an 80% training allocation is standard practice for clinical machine\nlearning applications.\nTo confirm that our results were not dependent on a single random partition, we repeated model\ntraining across multiple randomized splits through repeated-holdout cross validation. The\nresulting distribution from 1000 splits demonstrated that that latent features predict VNS\nresponse at AUCs significantly above chance (p < 0.001; mean = 0.613). Importantly, our\noriginal withheld dataset did not represent an outlier within this randomly permuted distribution\n(p ≈ 0.2, based on 1,000 random splits). These results support the stability of our model’s ability\nto predict VNS outcomes beyond what would be expected by chance alone.\n“To balance retaining a sufficient training set size with sufficient subjects for validation and\nrobust assessment on a held-out test cohort, these data were split into training, validation and\ntest sets in ratios of 75%, 10%, and 15% respectively, which is commonly used in the context of\nclinical datasets.30”\n“To assess the stability and significance of the learned features, we first trained 1000 SVM\nclassifiers on randomly partitioned subsets of the VNS cohort. The mean AUC across models\nwas significantly greater than that expected under the null distribution (mean AUC = 0.61; t =\n28.6, p = 0.001; Fig. 3f), suggesting that the latent representations carried reproducible\npredictive signal.”\nPET-derived neurotransmitter maps: The description of PET-based receptor/transporter profiles\nacross nine systems is intriguing. Were these data specifically acquired for this study, or part of a\npre-existing dataset? This is not fully described in the “pretraining data compilation” section.\nWe thank the reviewer for this comment. The PET-derived neurotransmitter receptor and\ntransporter maps were not acquired in our cohort, but come from previously published datasets\ncollated and made available through the neuromaps framework.30 We used these maps to provide\nneurochemical context for the model’s salience map, examining the associations across nine\nneurotransmitter systems.\n“To better understand the biological relevance of the regions in the VQ-VAE model’s salience\nmap, we studied the associations to the normative chemoarchitecture of the brain using PET-\nderived neurotransmitter receptor and transporter expression profiles across nine different\nneurotransmitter systems, as compiled from previously published datasets.31,32”\nResults / Discussion\nThe finding of serotonergic large-scale network disruptions is novel and important. However, the\nlack of results for the locus coeruleus–noradrenergic system is surprising, given its strong\nsupport in both human (P300, pupillometry) and animal studies (Raedt, Krahl). Could the authors\ncomment on this apparent discrepancy?\nIndeed, the locus coeruleus–noradrenergic (LC-NE) system has been strongly implicated in both\npreclinical and human studies of VNS, particularly through its role in arousal regulation, cortical\ngain modulation, and event related potential and pupillometric markers of noradrenergic tone.\nHowever, using PET-derived atlases, we did not find a significant alignment between our latent\nspace and norepinephrine transporter maps (p ~ 0.3). When visualizing and normalizing both\nPET maps to a common scale, we found that the NE distribution is concentrated largely within\nthe brainstem and thalamus (Supplementary Figure 13). Visually, there is some relation to the\nsalience map, however, this was not significant on spin permutation testing. This may be because\nthe serotonergic system exhibits a broader, more distributed projection pattern, spanning the\nthalamus, striatum, insula, and multiple cortical and subcortical targets. It is possible that the LC-\nNE contribution is under-represented in our latent space. The absence of LC-NE findings should\nnot be interpreted as evidence against noradrenergic involvement in VNS, but the evidence for\nserotonergic signalling should highlight its relevance and importance. In the revised manuscript,\nwe have expanded the discussion on these negative results.\n“Fourth, the VQ-VNS model was trained exclusively on T1w imaging. While more advanced\nimaging sequences could improve model performance, it would be at the expense of accessibility.\nFinally, while our salience maps aligned with the widespread serotonergic system, they did not\noverlap with the locus coeruleus-noradrenergic (LC-NE) system. This does not discount the\nestablished mechanistic role of the LC-NE pathway;33 which may be underrepresented within\nour latent space due to its highly focal nature (Supplementary Fig. 13). Thus, our findings\nhighlight the role of widespread serotonergic signaling without precluding the involvement of\nthe anatomically distinct LC-NE system.”\nSupplementary Figure 13. Comparison of PET-derived norepinephrine and serotonin\ntransporter maps. Norepinephrine transporter distribution shows a highly focal pattern largely\nrestricted to the brainstem and thalamus. In contrast, serotonin transporter exhibits a more\nspatially widespread projection architecture, encompassing striatal, insular, and other\nsubcortical and cortical regions. Min-max normalized maps highlight the relative difference in\nspatial extent. This distinction may contribute to an underrepresentation of noradrenergic\ncontributions to this current framework.\nThe statement that “beyond optimal programming, there are no associations between increased\ncharge density and outcome” (ref 59) may need qualification. Recent work, including Fahoum et\nal. (Brain Stimulation, 2022), suggests optimal intensity correlates with outcomes, especially in\npediatric cohorts. Please reconcile with this evidence.\nWe thank the reviewer for this comment. Within a subset of patients with stimulation parameters\navailable (n=114), we found no significant difference in stimulation parameters between\nresponders and non-responders (current: t = -1.2, p = 0.22; duty cycle: t = -1.30, p = 0.20).\nRecent analyses indicate that dose matters up to an evidence-based target, after which additional\nincreases in charge density do not consistently yield further benefit. Very interestingly, our\nfindings are aligned remarkably well with Fahoum et al. In their population-based study, a group\nmean of 1.6 mA and 17.1% duty cycle were associated with greater probability of response. Our\ngroup average current was 1.6 mA and an average duty cycle of 20%. Interestingly, we recently\nreported that strong pediatric responders (i.e. “super responders”) to VNS who become seizure-\nfree tend to do so at very low output currents (Tamura et al., Neuromodulation 2025).\n“This study has several limitations. First, while VNS implantation is highly standardized with\nminimal inter-subject variability, stimulation parameters were not standardized across sites and\nwere adjusted at the discretion of the treating clinician, leading to potential differences in the\ncharge density delivered. However, among participants with available data (N=103), mean\nstimulation parameters (output current 1.8mA; duty cycle 20%) were consistent with optimally\nprogrammed populations.26”\nDetails\nThe section on region-based morphometry would be clearer if you specify whether this refers to\nVBM, SBM, or other approaches.\nRegion-based morphometric analyses were performed using FreeSurfer-derived estimates of\ncortical thickness. We have clarified this detail in the Methods section entitled “Cortical\nThickness Estimation for Lesion Network Network Mapping.”\n“Cortical thickness was estimated from T1-weighted images using FreeSurfer as a surface-\nbased morphometry pipeline to derive individualized thickness maps used in all downstream\natrophy and network-based analyses.”\nThe statement “region-based morphometry is constrained by atlas-based parcellations…”\nrequires a proper reference. If reference 15 already supports this, please move it to the end of the\nsentence; if not, an additional reference is needed.\nWe thank the reviewer for this comment. This statement is supported by Guo et al. (2023). In\ntheir study of VNS patients, the authors note that common structural alterations such as cortical\nthickness and surface area could not be identified using traditional atlas-based morphometric\napproaches at the whole-brain level. They highlight this limitation as motivation for using\ndeformation-based morphometry, which provides voxel-wise sensitivity to distributed alterations\nacross the brain. We have clarified this in the revised text by moving the citation to the end of the\nsentence.\nConclusion:\nIn its current form, the manuscript is technically strong and demonstrates the use of highly\nspecialized methodologies.\nHowever, the presentation may be difficult for clinicians and readers outside the computational field to fully grasp. If the\nauthors aim for a high-impact journal and a broad audience, we encourage them to simplify and clarify the description of the\nmethods, and to place greater emphasis on the biological and clinical interpretation of their findings —\nparticularly the role of the implicated networks and neurotransmitter systems. This would\nconsiderably enhance the accessibility, readability, and overall impact of the work.\nWe thank the reviewer for this suggestion. In response, we have substantially revised the\nmanuscript to improve accessibility for a broader readership. Specifically, we now clearly state\nthe goal and rationale behind each technical step, define key methodological terms upon first use,\nand include concise synthesis statements throughout to guide interpretation. These revisions\nenhance readability and ensure that the biological and clinical relevance of our findings is\nemphasized alongside the technical innovation, resulting a more cohesive and broadly accessible\nmanuscript.\nSummary Methodological/Engeneering\nThe manuscript introduces the VQ-VNS framework, a methodology designed to predict and\nclassify responses to vagus nerve stimulation (VNS) in pediatric populations. The framework\nbuilds on the VQ-VAE architecture and incorporates interpretability mechanisms by adapting\nSHAP and Grad-CAM. Identifying likely responders has the potential to minimize the need for\ninvasive procedures.\nA key limitation of the study is that the proposed methodology is presented with substantial\ntechnical complexity, but with limited benchmarking against alternative approaches.\nMajor Comments\n- The authors reference several alternative approaches for responder identification [11–15]. A\ncomparative table would help quantify the accuracy improvements of the proposed framework,\nparticularly in relation to other T1-weighted MRI–based predictive models.\nTo our knowledge, there are currently no other predictive models of VNS response derived\nexclusively from routine T1-weighted MRI. Prior approaches have generally relied on\nspecialized acquisitions or advanced microstructural metrics. For example, Berger et al. reported\nstrong performance using diffusion MRI and tractography-based features, which required\nextended scan protocols and complex pre-processing pipelines not routinely available in clinical\nsettings. By contrast, our VQ-VNS framework operates directly on standard clinical T1-weighted\nscans, without any additional imaging modalities or advanced preprocessing. Arya et al.\nidentified clinical correlates of response, such as the absence of a lesion and a longer epilepsy\nduration but did not develop an imaging-based predictive model. Our model excluded images\nwith significant anatomical abnormalities. Other articles either did not report model performance\nor did not train a predictive model. We have included below a summary of the cited articles\nbelow.\nTable R1. Comparison of previously published studies on imaging-based VNS response.\nPublication Cohort Input features Model type AUC Comments\nBerger et al., 18 adults Diffusion SVM with 0.88 Small dataset;\nNeurotherapeu with DRE MRI (DTI recursive extensive\ntics 2024 and NODDI) feature preprocessing;\nelimination advanced imaging\ntechniques\nArya et al., 43 with Clinical and Logistic 0.80 No image-based\nJournal of childhood- MRI lesion regression predictive modeling;\nChild onset status lesions are rare\nNeurology epilepsy\nSeizure 2010 DRE EEG, and logistic reported multimodal clinical\n(pediatric and MRI regression and EEG data;\nadult) lesions are rare\nGuo et. al., 93 Clinical, MRI Logistic Not No image-based\nFrontiers in encephalomal lesion regression reported predictive modeling;\nNeurology acia-related laterality, and lesions are rare\nXie et al., 95 pediatric Clinical and Descriptive N/A Associations only;\nEpilepsy and DRE, MRI lesion statistics no predictive\nBehaviour structural characteristic (no classifier was\n2023 etiology s predictive trained\nmodel)\n- The manuscript states: “The resulting archetypal non-responder showed reduced grey matter,\nwhite matter, and subcortical volumes, together with ventricular enlargement.” If such\nmorphological alterations are discriminative, why would automated segmentation of T1-\nweighted images not provide an equally effective classification of responders? A comparison\nwith standard volumetric and morphometric approaches would strengthen the contribution of the\nproposed framework.\nWhile the archetypal reconstructions highlight group-level morphological differences, these\nalterations are not consistently present at the level of individual subjects. Indeed, it is the failure\nof these models and metrics that inspired the current representation learning approach to VNS\noutcome prediction. In the revised manuscript, we have included our analysis of basic imaging-\nrelated metrics using conventional modelling techniques, the inadequacy of which underscores\nthe superiority of more sophisticated, state-of-the-art techniques like VQ-VNS. These findings\nare concordant with the central thesis of the paper, which is that latent space representations\nbased on generative models designed to learn compact representations of brain patterns are\nrequired to understand complex associations outcomes to VNS.\nInformed by differences highlighted by the archetypal reconstructions, we trained both a logistic\nregression model and a support vector machine using automated segmentation of T1-weighted\nimages alone. These models, trained exclusively on FreeSurfer-derived gray matter, white\nmatter, subcortical gray matter, and ventricular volumes, demonstrated near chance performance\n(logistic regression: AUC = 0.49; SVM: AUC = 0.55) in predicting VNS response. This suggests\nthat, individually, automatically segmented features on T1w imaging are not on their own\npredictive of response. The strength of the proposed framework lies in its ability to learn internal\nlatent representations that capture subtle, distributed, and non-linear structural features across the\ncohort. These features are not easily detectable through conventional volumetric or\nmorphometric pipelines. While the differences between archetypal responders and non-\nresponders helps to explain some of the characteristics that may drive classification, they are\ninsufficient to enable classification on an individual subject level.\n- The original spin permutation procedure rotated cortical surface maps while preserving\ncontralateral symmetry. Was this constraint maintained in the present experiments? Furthermore,\nthe original test statistic excluded non-cortical regions. Were such regions incorporated in the\ncurrent analysis? These methodological details require clarification.\nOur approach was inspired by the originally described spin permutation framework35 but was not\nidentical in implementation. Instead of applying surface-based rotations, we implemented a 3D\nrotational affine transformation directly to the latent brain volumes. This procedure preserves the\nrelative spatial arrangement of all voxels (i.e., no voxel-wise shuffling is performed), thereby\nmaintaining bilateral or contralateral symmetries present in the data. Additionally, our\nimplementation of spin permutation includes both cortical and subcortical structures. This\nbroader inclusion was intentional, as our latent representation captures patterns beyond just the\ncerebral cortex.\nMinor Comments\n- The statement that prognostic signatures are “likely embedded in subtle, distributed, and\nnonlinear anatomical patterns” would benefit from supporting references.\nThank you for allowing us to clarify this claim. Prior studies that have successfully predicted\nVNS outcomes using advanced neuroimaging (e.g. diffusion tensor imaging, functional MRI,\nmagneotencephalography, etc.), as in (Berger et al., 2024; Ibrahim et al., 2017; Mithani et al.,\n2019), have all identified subtle, distributed features predictive of response. Furthermore, most of\nthem have used non-linear modelling for classification. Conversely, studies that have tried to use\ngross anatomical features for VNS prediction have had limited success. As such, we propose\nthat, based on these prior works, predictive features on anatomical imaging are also likely to be\nsubtle, distributed, and non-linear. We have added citations to this claim now.\n- The rationale for selecting an SVM classifier should be expanded, in comparison with other\nmachine learning methods (e.g., random forests, k-nearest neighbors).\nOur choice of SVM as the classifier for the learned embeddings was guided by both theoretical\nconsiderations and empirical evaluation. Among several standard classifiers tested, including\nlogistic regression, random forest, and k-nearest neighbors (kNN), SVM consistently provided\nthe best balance between model complexity, robustness to overfitting, and predictive accuracy on\nheld-out data. Logistic regression tended to underfit, random forests showed overfitting\ntendencies, and kNN performed less favorably overall.\nThese comparisons are now presented in Supplementary Figure 7 which confirms that SVM\nachieved the highest classification accuracy for the learned latent representations. Together,\nthese results support its selection as the final classifier in the VQ-VNS framework.\nSupplementary Figure 6. Performance of support vector machine, logistic regression, k nearest\nneighbours (KNN), and random forest classifier in predicting response to VNS from latent\nfeature representations of T1w MRI.\n- The contribution of individual components within the composite loss function could be clearer.\nWe thank the reviewer for this suggestion. In response, we have added a supplementary table to\nbetter describe the contributions of each individual component in the composite loss function.\nSupplementary Table 11. Summary of composite loss function terms.\nLoss term Description Formula\nVoxel-wise Mean absolute error ℒ = 𝑀𝐴𝐸(𝑋,𝑋̂)\nvox\naccuracy (MAE) between the\noriginal and\nreconstructed images,\nensuring voxel-level\nfidelity\nSpectral Mean squared error ℒ = 𝑀𝑆𝐸(𝐹𝐹𝑇(𝑋),𝐹𝐹𝑇(𝑋̂))\nfreq\nalignment (MSE) between the\nFourier transforms of\nthe original and\nreconstructed images,\npromoting frequency-\ndomain consistency\nPerceptual Learned Perceptual ℒ = 𝐿𝑃𝐼𝑃𝑆(𝑋 ,𝑋̂ )+ 𝐿𝑃𝐼𝑃𝑆(𝑋 ,𝑋̂ )\npercep ax ax cor cor\nsimilarity Image Patch Similarity + 𝐿𝑃𝐼𝑃𝑆(𝑋 ,𝑋̂ )\nsag sag\n(LPIPS) computed on\n2D axial, coronal, and\nsagittal slices using\nAlexNet as the feature\nextractor\nEncoder- Enforces alignment 2\nℒ = 𝽽‖𝑧 (𝑋)−𝑠𝑔(𝑒)‖\nquantizer between encoder 𝑐𝑜𝑚𝑚𝑖𝑡 𝑒 2\nalignment outputs and discrete\n(commitment) codebook vectors using\na stop-gradient\noperator; 𝽽 controls\nscaling\nAdversarial Least squares loss ℒ = (𝐷(𝑋̂)−1) 2\nadv\nobjectives guiding the generator to\nproduce realistic\nreconstructions\nindistinguishable from\nreal images by a patch-\nbased discriminator\nAblation experiments comparing performance with and without each term could justify their\ninclusion.\nThe composite loss function we adopted was based on Tudosiu et al., 2024, where the authors\nperformed extensive ablation experiments to evaluate the contribution of each loss component.\nAs this prior work has already rigorously demonstrated the necessity and impact of each term,\nwe did not repeat these ablations here.\nWe have included this detail under “Self-supervised representation learning captures anatomical\npatterns associated with VNS outcome” in the Results:\n“The composite loss function we adopted was based on prior experiments43 where extensive\nablations were performed to evaluate the contribution of each loss component.”\n- The claim that “the model’s reconstructions preserved fine anatomical details” appears\nsomewhat contradictory, since VQ-VAE models are known to lose fine-scale structural\ninformation. A more precise description of which features are preserved or lost would enhance\ninterpretability.\nThank you for allowing us to clarify this point. Our implementation was explicitly evaluated for\nanatomical fidelity using a multi-modal set of metrics. As detailed in the section entitled “Self-\nsupervised representation learning captures anatomical patterns associated with VNS outcome,”\nwe quantified reconstruction fidelity with MS-SSIM, FID, and MSE, and further assessed\nanatomical plausibility through structure-wise correlations after SynthSeg segmentation. These\nanalyses demonstrated that the model preserved key anatomical details, with correlations across\ntissue classes ranging from 0.96 to 1.00. Although reconstruction fidelity was somewhat lower in\nsubjects with pronounced structural abnormalities, overall performance remained high,\nconsistent with training on predominantly anatomically typical brains. Our claim is therefore\nspecific to the present implementation, supported by the reported metrics.\n- The observation that “reconstruction fidelity was lower in subjects with significant structural\nabnormalities” raises the question of whether this reduction affects classification performance.\nAre response predictions less reliable for patients with pronounced abnormalities?\nThe VQ-VAE was pretrained exclusively on brains without significant structural abnormalities,\nwhich explains the lower reconstruction fidelity observed in cases with pronounced structural\nabnormalities. To avoid introducing systematic biases, patients with major abnormalities were\nexcluded from the classification analyses. As a result, the reduction in reconstruction fidelity in\nsubjects with significant structural abnormalities would not have biased classifier performance.\n- Additionally, the reconstruction fidelity of the CONNECTiVOS pediatric epilepsy cohort\nshould be reported.\nWe have now included a supplementary table reporting the MS-SSIM and MSE scores for the\nCONNECTiVOS dataset.\nSupplementary Table 5. Reconstruction scores for the CONNECTiVOS dataset.\nImage MS-SSIM MSE\n1 0.962 0.0035\n2 0.966 0.0026\n3 0.971 0.0025\n4 0.975 0.0023\n5 0.966 0.0029\n6 0.967 0.0033\n7 0.967 0.0032\n8 0.946 0.0064\n9 0.970 0.0027\n10 0.960 0.0041\n11 0.954 0.0039\n- Gradient-based saliency mapping seems to implicate more regions than those emphasized in the\ntext (thalami, brainstem, insular cortex, peri-Rolandic regions). How were the reported regions\nidentified and selected (e.g., thresholding criteria, local maxima)?\nFor interpretability, we focused our description on regions corresponding to the vagus afferent\nnetwork, which has been implicated in mediating vagus nerve stimulation effects and\nresponsiveness43. All regions highlighted meet a minimum threshold of 0.30 on a 0-0.5 log scale,\nrepresenting approximately the top 10% of voxels. We have provided a thresholded map of\nsaliency values below.\nFigure R1. Model salience map thresholded at 0.3, representing the top 10% of values.\n- Differences between synthetic responders and non-responders in Figure 2F are visually subtle.\nEnhancing or highlighting these differences would improve interpretability.\nGiven the subtle differences between the images, we provided a quantitative comparison in the\naccompanying bar plot, which more clearly highlight the volumetric distinctions between the\narchetypal responder and archetypal non-responder. To avoid confusion, we have revised the\nfigure caption to explicitly state that the differences are subtle and reconfigured panel e to better\nemphasize them.\nFigure 3: Performance of the VQ-VNS Quantized Autoencoder Model for T1w Reconstruction and\nClassifier for Predicting VNS Outcomes. (a) Representative reconstruction of a healthy subject from the\nheld-out test set of the VQ-VAE pretraining dataset. (b) Correlations between true and reconstruction\nvolumes across four tissue classes in the pretraining test set, demonstrating high reconstruction fidelity.\n(c) Representative reconstruction of a subject with a previous resection from the CONNECTiVOS cohort.\n(d) Synthetic images of the archetypal VNS responder and non-responder reconstructed from the latent\nspace using the VQ-VAE decoder. Differences are subtle and not visually apparent. (e) The synthetic\narchetypal non-responder showed reduced grey matter, white matter, and subcortical grey matter\nvolumes, and ventricular enlargement relative to the synthetic archetypal responder. (f) 1000 models\ntrained using hold-out cross validation perform significantly better than chance. (g) Receiver operating\ncharacteristic (ROC) curve showing the predictive performance of the SVM classifier on the held-out test\nset (AUC = 0.73, p = 0.007), which was specific to VNS and did not generalize to outcomes following\ntemporal lobectomy (AUC = 0.57, p =0.12).\n- “The model’s attention was concentrated within brain regions enriched in serotonin transporter\n(5-HTT) expression”. Can you quantify this result?Clarification is needed.\nWe found that the concentration of model attention within the brain was enriched in serotonin\ntransporter (5-HTT) expression based on normative PET data. This was quantified by computing\nthe correlation between attention weights and regional 5-HTT density. The correlation\ncoefficient and corresponding p-value are more explicitly denoted in the revised results.\n“We generated 1,000 such null maps per comparison to compute empirical p-values and account\nfor spatial autocorrelation. The VQ-VAE’s attention was concentrated within brain regions\nenriched in serotonin transporter (5-HTT) expression (r = 0.26, p = 0.04; Fig. 4b).”\n- “To assess the predictive value of clinical data, a support vector machine (SVM) was trained on\n80% of the dataset (discovery cohort, n = 550) to predict VNS outcome, with the remaining 20%\nheld out for validation (n = 144)”. Was there an independent test set? Clarification is needed.\nIn our analysis of the predictive value of clinical features, the SVM was trained exclusively on\nthe 80% training set. The remaining 20% was never used during model development and was\nheld out as an independent test set to evaluate performance after training was complete. We have\nrevised the text to clarify this distinction.\n“To assess the predictive value of clinical data, a support vector machine (SVM) was trained on\n80% of the dataset (n = 550) to predict VNS outcome. The remaining 20% (n = 144) was held\nout as an independent test set to evaluate model performance after training.”\n- “Structure-wise agreement was evaluated across four tissue classes – grey matter, white matter,\nsubcortical grey matter, and ventricles – by computing the Pearson correlation between\nreconstructed and original volumes across subjects”. This only accounts for the volume size of\neach tissue type and not their location. Adding the DICE scores or an equivalent would\nprovide a more comprehensive assessment.\nWe have now computed Dice similarity coefficients between the reconstructed and original\ntissue maps for each class. The results show high agreement across tissue classes, consistent with\nthe strong volumetric correlations reported previously.\nSupplementary Table 6: Dice-Sørensen coefficients\nTissue Class Mean Dice Score 95 % CI\nCSF 0.85 0.0035\nGrey Matter 0.78 0.0032\nSubcortical Grey 0.91 0.0024\nWhite Matter 0.80 0.0036\n- “Atrophy networks were defined by combining endpoint maps across all HCP subjects”. Why\nwere endpoints chosen instead of complete streamlines?\nEndpoints were chosen instead of complete streamlines because our objective was to map direct\nstructural disconnections resulting from regional atrophy, rather than to characterize entire\ntrajectories. Full streamlines often traverse multiple white-matter regions and are susceptible to\ngyral bias and false continuations, making them less specific to true cortical terminations.44\nTerminations. Structural connectivity derived from streamline endpoints provides more accurate\ncortical coverage, reduces gyral bias, and yields higher reproducibility compared to full-tract\nmeasures.28 Thus, using streamline endpoints provides a biologically grounded estimate of the\nfirst-order disconnection pattern associated with each atrophic node. This endpoint-based\napproach has been used in prior structural connectivity studies29.\n- A visualization of the cortical atrophy map should be included for completeness.\nThank you for this suggestion. A visualization of the cortical atrophy map was provided in\nSupplementary Figure 10c-d, copied below for reference.\nSupplementary Figure 10: Size and Spatial Distribution of Atrophy. (a) Distribution of atrophy size across all\nsubjects, relative to the 95th percentile and the 5% exclusion threshold. (b) Comparison of atrophy size between\nresponders and non-responders shows no significant group difference. (c) Spatial distribution of atrophy across the\nfull cohort. (d) A small focal cluster of atrophy is preferentially observed in non-responders, with no reciprocal\npattern in responders."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Authors have addressed my comments well. Congratulations on the nice work!", "reviewer_2": "The authors have adequately addressed concerns raised and incorporated suggestions to significantly improve the\nmanuscript.\nMinor suggestions\n- Lines 237 and 303. When using the term “predicted” please add qualifiers to be clearer in specifying the level of predictive\nperformance (predictions can be correct, incorrect, above chance, at chance, or below chance). For example:\n+ “The SVM model did not predict outcomes above chance level after…”\n+ “...model demonstrated predictive accuracy above chance levels across”\n-Line 220, Change “reproducible” to “measurable” or “detectable” as it better aligns with the scope of the outcomes (AUC ~ 0.61 for N=1000 repeats).", "reviewer_3": "All comments have bene adequately answered", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We thank the reviewer for their careful evaluation of our manuscript. We are pleased to hear that\nall comments have been satisfactorily addressed. We greatly appreciate the reviewer’s feedback\nthroughout the process.\nWe thank the reviewer for their careful comments. We have revised the text at lines 237 and 303\nto explicitly qualify predictive performance relative to chance level.\nThank you for your comments. The HuggingFace link has been added with the GitHub link in\nthe manuscript.\nWe thank the reviewers for their thoughtful feedback. We appreciate the reviewers’\ncontributions, which have helped us improve the manuscript."}]

Decision letter