Novel multimodal neurostimulation protocol (MNS protocol)

As shown in Fig. 1, the novel MNS protocol is a sequential stimulation approach that combines personalized tACS and robot-assisted navigated rTMS. The personalized tACS intervention is administered at 40 Hz for a duration of 30 min. Immediately following the tACS, neuro-navigated rTMS is performed. This protocol leverages individualized DTI data to precisely identify fiber tracts between the DLPFC and dorsal anterior cingulate cortex (dACC), selecting the point of highest structural connectivity as the stimulation site. By optimizing stimulation parameters for both tACS and rTMS based on individualized head model, the MNS protocol aims to enhance cognitive function in BD patients by utilizing advanced neuroimaging and targeted stimulation techniques. The specific method is described as follows.

Fig. 1

Schematic diagram of multimodal neurostimulation protocol. Individual DTI data are analyzed to determine the neural fiber tract density between the DLPFC and dACC. The point with the highest connectivity density is identified as the stimulation site. Using individual T1-weighted structural images, personalized induced EF models for tACS and rTMS are constructed. The tACS parameters are then optimized ensuring precise stimulation in the stimulation site, followed by the application of navigated rTMS

Individualized stimulation site localization based on DTI imaging

For individualized stimulation site localization based on DTI imaging, as shown in Fig. 1, T1-weighted structural images and DTI data were acquired from patients and subsequently preprocessed. The detailed scanning parameters can be found in the “Imaging” section. The processing of T1 images utilized Freesurfer software for automated region-of-interest segmentation, including image grayscale normalization, gray and white matter segmentation, and fine cortex segmentation using the Desikan-Killiany template. The DTI images were processed using MRtrix3 and FSL software, including image noise reduction, Gibbs artifact removal, and eddy current correction. Secondly, after calculating the response function and CSD fiber directional distribution function, a fiber tracking algorithm based on the second-order FOD was used to reconstruct one million fibers in the whole brain. Subsequently, the two pre-segmented regions of interest DLPFC and dACC in the T1 structural image were extracted, registered into the DTI space, the fiber bundles passing through the two regions at the same time were screened, and the fiber density in each voxel in the fiber path was calculated as a measure of the structure, an indicator of the strength of the connection. Finally, the voxel position corresponding to the highest point of structural connection strength in the individual DTI space was converted to the MNI152 (Montreal Neurological Institute) standard T1 template space after nonlinear registration as the stimulation site. Detailed information on the stimulation target location for each individual can be found in the Supplementary Material.

Treatment protocol of personalized tACS

After determining the individualized stimulation site, we proceeded to optimize the treatment parameters of tACS for each patient as depicted in Fig. 1. We optimize the electric field intensity and focus of the target stimulation point. While maximizing the electric field intensity in the target area, it minimizes the electric field intensity in non-target areas and improves the focality of the electric field. In addition, some current constraints are set for safety reasons, including the size limit of the total electrode current and the size limit of a single electrode current. The specific steps are as follows.

First, the SIMNIBS tool software was applied for constructing individualized head models. The selected stimulation electrode employed in this study is a 10-mm diameter disc electrode. The electrical conductivity characteristics of each corresponding tissue and material are set to 0.465 S/m (scalp), 0.01 S/m (skull), 1.65 S/m (cerebrospinal fluid), 0.276 S/m (gray matter), 0.126 S/m (white matter), 2.5e − 14 S/m (air cavity), 5.9e7 S/m (electrode), and 0.3 S/m (conductive paste).

The proposed electrode optimization approach for tES involves the following sequential steps. Primarily, the initial step involves the calculation of the leadfield matrix, denoted as A:

where:

$$e=\begin\begine\left(r_1\right)\\e\left(r_2\right)\\:\end\\e\left(r_N\right)\end,\;A=\left[\begin\alpha_1\left(r_1\right)\\\alpha_1\left(r_2\right)\\:\\\alpha_1\left(r_N\right)\end\;\begin\alpha_2\left(r_1\right)\\\alpha_2\left(r_2\right)\\:\\\alpha_2\left(r_N\right)\end\;\cdots\;\begin\alpha_M\left(r_1\right)\\\alpha_M\left(r_2\right)\\:\\\alpha_M\left(r_N\right)\end\right],s=\begins_1\\s_2\\:\\s_M\end$$

rN represents the node position on the head model, N denotes the number of nodes in the head model, and M signifies the total number of candidate electrodes.

To enhance computational efficiency, the EF intensity of the minimized non-target region is quantified as the EF intensity on the voxel point situated at a distance of 1/4* the maximum distance from the target point. The objective of parameter optimization is expressed through the optimization function:

$$\arg\;\max\nolimits_s\left(e_0^TCs\mathit-\lambda E__-\mathrm\mathrm x}\right)$$

(2)

where s denotes the magnitude of the stimulation current, e0TC signifies the distribution of the required EF intensity at the target, and E0.25d_max denotes the EF intensity at a node located at 1/4* the maximum distance from the target point. The weight parameter λ is used to balance the two objectives of focality and EF intensity.

Constraint conditions are integral to the optimization method. Firstly, the total input current of the electrode must equal the output current, ensuring a balanced current sum.

Additionally, the sum of the absolute values of the current is not larger than Itotal.

$$\sum\left|s_i\right|\leq I_$$

(4)

and the absolute value of the current of a single electrode is not greater than Imax.

$$\left|s_i\right|\leq I_$$

(5)

for clinical applicability, constraints limiting the number of electrodes (denoted as \(n\)) are introduced:

and the absolute value of the current of a single electrode is not greater than Imax.

In this study, we set the Itotal at 4 mA, with a maximum current Imax of 2 mA per electrode. The number of electrodes \(n\) was selected to be either 4 or 5.

Treatment protocol of navigated rTMS

In this study, we combined the acquired individual structural MRI images to construct personalized head models using SIMNIBS software. Our head models encompassed essential components such as the scalp, skull, cerebrospinal fluid, gray matter, white matter, and other brain tissues. The figure-of-eight coil, featuring an inner diameter of 20 mm and an outer diameter of 70 mm, was integrated into the simulation model. The coil was positioned directly above the stimulation site and simulated as current sources. We conducted simulation modeling of the electric fields induced by TMS to evaluate how varying coil orientations affect the distribution of the electric field within the brain. Specifically, the TMS coil was rotated in 15-degree increments to systematically simulate the induced electric field at different angles. Based on the simulation results, the optimal coil orientation was selected for each participant to ensure that the target brain region received sufficiently strong electric field.

A neuro-navigation system was employed for rTMS. Firstly, the individualized head model was co-registered with a standard brain using AC-PC alignment. The MRI images were converted into three-dimensional skin and brain models to replicate the head characteristics of the patient. Using a visual localization system, a head positioning mask, and appropriate computer software, the treatment site was identified on the patient’s scalp. Subsequently, the treatment was administered by aligning the coil at the specified angle with the determined treatment site.

Sequential personalized tACS and navigated rTMS

The MNS protocol proposed in this study represents a departure from traditional single-modal electrical or magnetic stimulation techniques. This protocol employs a dual-modal approach, combining personalized tACS with navigated rTMS. The personalized tACS intervention is administered at 40 Hz for a duration of 30 min, including a 30-s ramp up at the start and a 30-s ramp down to 0 mA at the end. Immediately following the tACS, neuro-navigated rTMS is performed. Each daily rTMS treatment consisted of 60 five-second 10-Hz trains delivered at 100% of the resting motor threshold (RMT) with inter-train intervals of 15 s (i.e., 3000 pulses per session, 45,000 pulses overall). For the sham treatment, the coil was directed to the same stimulation site of scalp sensations with the production of the same sounds. Each participant received 15 days of tACS-rTMS treatment. This sequential application is designed to maximize the synergistic effects of electromagnetic stimulation on cognitive function.

The MNS protocol was conducted by a skilled technician. The transcranial electrical stimulation device (Starstim, Neuroelectrics, Barcelona, Spain) and Circular Ag/AgCl electrodes (PISTIM, Neuroelectrics, Barcelona, Spain) were utilized. Conductive gel (HD-GEL, Soterix, NJ, USA) was applied as the conducting medium. The rTMS procedures utilized the Magstim Rapid2 rTMS device (The Magstim Company, Whitland, UK). We employed the Black Dolphin Navigation Robot (S-50, a sub-millimeter smart positioning system, Solide Medical Sci. & Tech.Co., Ltd., Xi’an, China) with a figure-of-8 coil (Solide Brain Control Medical Technology, Xi’an, China) to perform the rTMS.

Trial design

As illustrated in Fig. 2A, to validate the efficacy of the proposed MNS protocol, a randomized controlled double-blind clinical trial was conducted. The stringent experimental design ensures that both participants and investigators remain blinded to the treatment allocation, thereby minimizing potential biases and enhancing the validity and reliability of the findings.

Fig. 2

Clinical trial design and intervention paradigm. A Participants underwent a 3-week multimodal intervention combining transcranial electrical stimulation and transcranial magnetic stimulation. Each week consisted of daily sessions, totaling 15 sessions. Cognitive assessments and imaging examinations were conducted before and after the intervention. B Group A received 30 min of active tACS followed by 20 min of active rTMS. Group B received 30 min of sham tACS followed by 20 min of active rTMS. Group C received 30 min of active tACS and sham rTMS. Group D received 30 min of sham tACS and sham rTMS. Each participant with recorded scores at three time points, W0: week 0, W3: week 3

Participants were randomly assigned to any one of the four groups: A, B, C, or D. In group A (active tACS and active rTMS (MNS protocol)), tACS is the active stimulation, and rTMS is the active stimulation group. In group B (sham tACS and active rTMS), tACS is the sham stimulation, and rTMS is the active stimulation group. In group C (active tACS and sham rTMS), tACS is the active stimulation, and rTMS is the sham stimulation group. In group D (sham tACS and sham rTMS), both tACS and rTMS are sham stimulation. The personalized localization was conducted for the four groups of patients, All the participants received a consecutive 3-week treatment, including five times per week.

As illustrated in Fig. 2B, the active rTMS group received rTMS treatment at 100% RMT, while the sham stimulation group received stimulation at 20% RMT. The active tACS group underwent 30 min of alternating current stimulation, whereas the sham tACS group had current passing only during the initial and final 30 s of the stimulation. Cognitive and imaging assessments were conducted before (week 0, W0) and after (week 3, W3) the treatment. Psychological scale evaluations were performed during the treatment, and adverse reactions were documented.

Participants

This study recruited patients with stable-phase BD and cognitive impairment from the outpatient Psychiatric Department of the First Affiliated Hospital, Zhejiang University School of Medicine, from July 2023 to February 2024. All participants enrolled in this study were clinically diagnosed with bipolar disorder type II. The research was approved by the Clinical Research Ethics Committee of the First Affiliated Hospital, Zhejiang University School of Medicine (Approval number: IIT20230058C-R1). The study has been registered in the Clinical Trials Registry (http://www.clinicaltrials.gov) with the registration number NCT05964777. Before participation, all participants provided informed consent, ensuring their understanding of the objectives, procedures, and potential risks of this study. Participant confidentiality and well-being were upheld throughout the study.

The inclusion criteria included: (i) aged 14–45, regardless of gender; (ii) diagnosis of BD in accordance with the Mini-International Neuropsychiatric Interview (M.I.N.I.); (iii) stable medication treatment for at least 3 months; (iv) clinical remission for at least 3 months before entering the randomization phase, with a Young Mania Rating Scale (YMRS) score ≤ 6 and 17-item Hamilton Depression Rating Scale (HAMD-17) score ≤ 7; (v) cognitive impairment indicated by a Perceived Deficits Questionnaire (PDQ) score ≥ 17; (vi) right-handed; (vii) education level ≥ 9 years.

The exclusion criteria included: (i) comorbidities of any other mental disorder in the Diagnostic and Statistical Manual of Mental Disorders, Fifth Edition (DSM-5); (ii) any contraindications for MRI or the discovery of structural abnormalities in brain MRI; (iii) pregnant, lactating women, breastfeeding or those planning to become pregnant; (iv) patients who have undergone ECT, rTMS, or tACS treatment in the last 3 months; (v) history of serious neurologic illnesses, such as epilepsy or traumatic brain injury; (vi) significant and unstable medical conditions, including diabetes or cardiovascular, hematological, endocrine, liver, or kidney disease; (vii) history of substance abuse and dependence; (viii) achromatopsia, hypochromatopsia, or dysaudia; (ix) medication with antidepressant or anticholinergic agents.

Imaging

Imaging was conducted using a 3.0 Tesla scanner (GE SIGNA) with a standard whole-head coil at the First Affiliated Hospital, Zhejiang University School of Medicine. Participants were instructed to lie on the scanner bed with their eyes closed. Foam cushions were placed on either side of the head to minimize head motion, and ear plugs were provided to reduce noise interference. High-resolution 3D anatomical images were acquired using a T1-weighted magnetization prepared rapid acquisition gradient-echo (MPRAGE) sequence. The scanning parameters were as follows: 3D T1 structural image: TR/TE = 7.1/2.9 ms; 146 slices, slice thickness 1 mm; flip angle = 8°; matrix size = 256 × 256; field of view = 260 × 260 mm2. For the DTI diffusion image, a planar echo sequence was used with the following parameters: TR/TE = 11,000/89.4 ms; matrix size = 256 × 256; field of view = 240 × 240 mm2; number of slices = 44; slice thickness = 3 mm. The b-value for DTI scanning was 1000 s/mm2, and 25 diffusion directions were included. During scanning, participants were instructed to keep their eyes closed. To minimize head motion, foam cushions were strategically placed on either side of the head, and earplugs were utilized to mitigate noise.

Imaging analysisfMRI analysis

fMRI data preprocessing was conducted using Data Processing Assistant for Resting-State fMRI (DPARSF) Advanced Edition V5.4 in MATLAB 2020. The key preprocessing steps included segmenting the MRI data into 10 periods, discarding the first 10 noisy time points, and synchronizing slice timing. Participants with head movements exceeding 3 mm or 3° were excluded from the analysis. Functional images were normalized to the MNI space using T1-based segmentation and registration. Spatial smoothing was applied using an 8 mm × 8 mm × 8 mm Gaussian kernel. Linear drift was removed through linear regression to minimize the effects of thermal noise. Band-pass filtering (0.01–0.1 Hz) was employed to eliminate physiological noise. The brain was parcellated into regions based on the AAL116 template, with each region treated as a network node. Pearson correlation coefficients of the time series were computed to construct brain network matrices, which were Fisher Z-transformed for subsequent statistical analysis. Graph theory analysis was conducted using GRETNA software to calculate global and nodal metrics of the brain network matrices for both the baseline and post-treatment groups over a sparsity range of 0.01 < S < 0.35 with a step size of 0.01. Global metrics included the small-world index, normalized characteristic path length, characteristic path length, clustering coefficient, global efficiency, local efficiency, and network efficiency. Nodal metrics comprised nodal efficiency, nodal local efficiency, nodal clustering coefficient, nodal shortest path, betweenness centrality, and degree centrality.

DTI analysis

The DTI data were processed using FSL software (FMRIB Software Library, Oxford, UK), which included motion and eddy current correction, brain extraction, tensor model fitting, and normalization to the MNI space using Tract-Based Spatial Statistics. First, motion and eddy current corrections were applied to each participant’s diffusion-weighted imaging data, and the diffusion gradient directions were rotated accordingly using affine registration. Subsequently, the diffusion tensor model was fitted using a weighted least squares approach in each participant’s native space. Spatial normalization involved mapping each participant’s FA image to a standard space (2-mm isotropic MNI152 template) using nonlinear registration. The FA reference and subsequent region of interest (ROI) definitions were based on the ICBM-DTI-81 white-matter atlas from Johns Hopkins University. Finally, scalar metrics from 40 ROIs were extracted from each participant’s projected skeleton.

Clinical assessment

The PDQ was utilized to screen for cognitive impairment in participants at baseline. Symptom severity was evaluated using the HAMD-17 and YMRS at different time points. Assessors underwent training for consistency prior to assessments, and the adverse events was recorded.

The THINC integrated tool (THINC-it) is a self-administered iPad-based assessment tool that has been demonstrated to be reliable and valid for evaluating cognitive symptoms in BD patients. THINC-it consists of Spotter, SymbolCheck, Codebreaker, Trails, and Perceived Deficits Questionnaire—5-item version for Depression (PDQ-5-D), corresponding to 5 traditional neuropsychological tests: the recognition task (CRT), Stroop task (1 back), the Digit Symbol Substitution Test (DSST), Digit Matching Line (TMT), and PDQ. These tests assess attention, processing speed, working memory, and cognitive and executive functions. Cognitive function in all participants was assessed using THINC-it at week 0 (W0) and week 3 (W3) for each participant. Any adverse events encountered were meticulously documented.

The primary outcome measure was the change in scores on the five THINC-it tests from week 0 and week 3. Secondary outcomes included changes in resting brain activity at the end of the 3rd week.

Clinical outcome analysis

The collected scores from the assessment tools were subjected to statistical analysis using factorial repeated measures analysis of covariance (ANOVA) and paired t-tests. Repeated measures ANOVA was employed to assess within-group and between-group differences across the two time points (W0 and W3). Paired t-tests were used to investigate specific changes between baseline and post-treatment or follow-up assessments.

This study employed MATLAB 2020b, SPM12, and DPABI for statistical analysis. SPM12 was used to perform paired t-tests to examine variations in fALFF, ReHo, and FC between groups. The DPABI software, in conjunction with Gaussian Random Field (GRF) theory, was utilized to correct for multiple comparisons in the paired t-test results. This correction involved a voxel threshold of P < 0.01 and a clump threshold of P < 0.05.