Personalizing repetitive transcranial magnetic stimulation for precision depression treatment based on functional brain network controllability and optimal control analysis

and psychiatric disorders such as Depression. How- ever, due to patient heterogeneity, neuromodulation treatment are often highly variable, requiring patient-speciﬁc stimulation protocols throughout the recovery stages to optimize treatment outcomes. Therefore, it is critical to personalize neuromodulation protocol to optimize the patient-speciﬁc stimulation targets and parameters by accommodating inherent interpatient variability and intersession alteration during treatments. The study aims to develop a personalized repetitive transcranial magnetic stimulation (rTMS) protocol and evaluate its feasibility in optimizing the treatment eﬃciency using an existing dataset from an antidepressant experimental imaging study in depression. The personalization of the rTMS treatment protocol was achieved by personalizing both stimulation targets and parameters via a novel approach integrating the functional brain network control- lability analysis and optimal control analysis. First, the functional brain network controllability analysis was performed to identify the optimal rTMS stimulation target from the eﬀective connectivity network constructed from patient-speciﬁc resting-state functional magnetic resonance imaging data. The optimal control algorithm was then applied to optimize the rTMS stimulation parameters based on the optimized target. The performance of the proposed personalized rTMS technique was evaluated using datasets collected from a longitudinal antidepressant experimental imaging study in depression ( n = 20). Simulation models demonstrated that the proposed personalized rTMS protocol outperformed the standard rTMS treatment by eﬃciently steering a depressive resting brain state to a healthy resting brain state, indicated by the signiﬁcantly less control energy needed and higher model ﬁtting accuracy achieved. The node with the maximum average controllability of each patient was designated as the optimal target region for the personalized rTMS protocol. Our results also demonstrated the theoretical feasibility of achieving comparable neuromodulation eﬃcacy by stimulating a single node compared to stimulating multiple driver nodes. The ﬁndings support the feasibility of developing personalized neuromod- ulation protocols to more eﬃciently treat depression and other neurological diseases.


Introduction
Brain neuromodulation has shown efficacy in treating numerous psychiatric diseases such as Parkinson's disease (PD) and Depression ( Lewis et al., 2016 ;Fenoy et al., 2018 ;Fang et al., 2021 ). For example, repetitive transcranial magnetic stimulation (rTMS) is a noninvasive neuromodulation tool approved by the US FDA to treat depres-sider patient heterogeneity and patient-specific recovery curve during treatment sessions in the stimulation montage. A recent Stanford Neuromodulation Therapy (SNT) protocol considered the patient heterogeneity by employing the resting-state functional magnetic resonance imaging (fMRI) signal to personalize the stimulation target located within the left dorsolateral prefrontal cortex (DLPFC) and adjust the stimulation intensities by the cortical depth of the stimulation target to ensure the delivery of 90% RMT ( Cole et al., 2021 ). However, the stimulation target was still restricted within the DLPFC and the stimulation parameters were not optimized across treatment sessions ( Cole et al., 2021 ). In order to maximize the rTMS efficacy in depression treatment, it is critical to personalize the rTMS protocol for each depressive patient by simultaneously optimizing the stimulation target(s) and parameters under different treatment sessions ( R.F. Cash et al., 2021 ;Figee and Mayberg, 2021 ;Fitzgerald, 2021 ;Klooster et al., 2022 ;Modak and Fitzgerald, 2021 ). Unfortunately, such a personalized rTMS technique, in which both stimulation target(s) and parameters will need to be personalized over treatment sessions for precision depression treatment, is not currently available because of the incomplete understanding of the neural control mechanisms and its alterations in depression and during the interventional process.
Personalized rTMS treatment could, in principle, involve the stimulation target(s) and parameters being optimized iteratively for each specific patient and each specific treatment session. A potential framework for such optimization is Network Control Theory, which has been developed to characterize the brain's dynamic properties, assuming the brain system is deliberately shifted or guided along a particular trajectory to support specific goals ( Medaglia et al., 2017 ). As such, it provides a potentially powerful tool for modeling the neuromodulation process in humans ( Tang and Bassett, 2018 ). Prior studies demonstrate that network control diagnostics can predict neural activations in response to exogenous stimulations, where the brain dynamics are constrained by a static, structural brain network built from the white matter tracts characterized using diffusion tensor imaging (DTI) data ( Bassett and Sporns, 2017 ;Gu et al., 2015 ). For example, the average controllability is defined as the controllability diagnostic evaluating a node's ability to drive the network system into easy-to-reach states. In contrast, modal controllability is defined as the controllability diagnostic assessing a node's capability to steer the network system into difficult-to-reach states ( Gu et al., 2015 ). A previous simulation study demonstrates that stimulating brain regions with high average controllability contributes to easily moving the network system into nearby brain states (brain states with high similarity such as resting state vs. an alternate resting state, or resting state vs. sleeping state) with little control energy ( Muldoon et al., 2016 ).
Moreover, the nodes with high average controllability are shown to be predominantly located in the default mode network (DMN), while high modal controllability nodes being mainly located in the frontoparietal network (FPN) ( Muldoon et al., 2016 ). These results provide a control framework for developing stimulation protocols by selecting optimal target nodes for specific brain state transition (for example, nearby states transition) based on the average and modal controllability diagnostics. While powerful, these approaches usually disregard other nonstructural components of the brain dynamics ( Deng and Gu, 2020 ), limiting their ability to accommodate the short-term changes of the underlying neural control patterns associated with the brain reorganization and their applications in personalizing the brain time-variant neuromodulation protocols to treat psychiatric diseases like depression.
Recently, a functional brain network controllability analysis (fBNC) method was developed to overcome this limitation ( Deng and Gu, 2020 ). The fBNC method characterizes brain controllability properties by studying network dynamics based on the effective connectivity matrices ( Ljung, 1999 ;M. Gilson et al., 2016 ), which can be constructed from functional neuroimaging data such as functional magnetic resonance imaging (fMRI) ( Deng and Gu, 2020 ). Conventional methods for effective connectivity estimation, including dynamic causal modeling ( Stephan and Roebroeck, 2012 ), granger causality ( Seth, 2007 ), and partial correlation-based models ( Marrelec and Benali, 2009 ), do not estimate the dynamics that occur on top of the network. In contrast, the fBNC approach allows us to consider the effective relations between brain regions and the brain network dynamics simultaneously, based on the system's input and output signals ( Ljung, 1999 ). Prior studies have demonstrated accurate estimation of the effective connectivity network and a good fitness of the fBNC approach in constructing the brain network dynamics ( Deng and Gu, 2020 ;M. Gilson et al., 2016 ). The fBNC analyses showed a consistent relationship between the controllability diagnostics, average and modal controllability ( Deng and Gu, 2020 ), as demonstrated in the structural brain controllability analysis ( Gu et al., 2015 ). These findings indicate the potential use of the functional time courses to explore the brain's underlying neural control mechanisms and to locate the optimal target regions based on the fBNC-derived controllability diagnostics. This also represents a significant enhancement of the approach through incorporation of functional properties and the shortterm changes of the underlying neural control patterns in contrast with prior methods that only took into account structural brain network information ( Muldoon et al., 2016 ;Beynel et al., 2020 ;Karrer et al., 2020 ).
Stimulation parameters could also be optimized for each individual patient and specific treatment session to personalize the rTMS technique for maximal neuromodulation efficacy ( Mohan et al., 2020 ). Optimal control theory is a branch of mathematical optimization that deals with finding optimal control parameters for a dynamic system over a period of time such that an expected network state can be reached ( Ross, 2009 ). By minimizing both the control energy and the difference of the initial state from the target state, optimal control theory can be employed to predict the optimal stimulation parameters required to drive the brain state from an initial state to the desired target state through a series of state transitions along an optimal trajectory ( Stisoet al., 2019 ). These predictions have been quantified for direct electrical stimulation to the structural brain networks to improve memory ( Stisoet al., 2019 ). However, it is clear that formating the controllability estimates with static structural networks fails to capture dynamic brain changes throughout the treatment sessions ( Battaglia et al., 2012 ). Optimal control theory has recently been extended to explore the optimal control energy required for seizure control in epilepsy patients within an effective local network ( Scheid et al., 2021 ). However, optimal control analysis was limited to the temporal lobe with implanted electrodes ( Stisoet al., 2019 ;Scheid et al., 2021 ). To date, generalization to the control of global brain dynamics and the consideration of the underlying dynamic neural alterations across the whole brain has not yet been thoroughly investigated.
The current study proposes a personalized rTMS protocol by fusing fBNC and optimal control analysis for patient-specific precision depression treatment. The performance of the proposed personalized rTMS technique was evaluated with the resting-state fMRI datasets collected in an antidepressant experimental imaging study in depression. The resting state fMRI (rs-fMRI) data were recorded at pre-treatment (visit 1), 1-week post-treatment (visit 2), and 6-weeks post-treatment (visit 3) and utilized to construct the effective brain networks across the entire brain regions and perform fBNC-based optimal control analysis. This work offers empirical support for developing and utilizing functional brain controllability analysis and optimal control approach in personalizing any neuromodulation protocols for neurological diseases' intervention and understanding the neural control impairments underlying neurological diseases.

Study design
The study was approved by Oxford Research Ethics Committee and was conducted as per the Declaration of Helsinki. Participant recruitment is detailed in ( Godlewska et al., 2016 ). Structured Clinical Interview for DSM-IV ( First et al., 1995 ) was used to ascertain the diagnostic status and for the presence of current and past psychiatric disor-ders. Depressed patients met DSM-IV criteria for a primary diagnosis of major depressive disorder; exclusion criteria for both groups were: past/current DSM-IV diagnosis of axis I disorder (other than MDD in patient group), substance dependence as defined by DSM-IV, a clinically significant risk of suicidal behavior, major somatic or neurological disorders, pregnancy or breast-feeding, any contra-indications to MR imaging, or concurrent medication that could alter emotional processing. All participants were right-handed. Patients with MDD received escitalopram 10 mg each morning for six weeks without dose adjustement and had MRI scans at baseline (before treatment), one week, and six weeks post-treatment. Healthy volunteers were tested once and did not receive treatment. In total, written informed consent to participate in the study was obtained from 39 patients with MDD (24F:15 M) and 32 healthy individuals (18F:14 M).
Depression severity was measured using the 17-item Hamilton Depression Rating Scale (HAMD) ( Hamilton, 1960 ) and Beck Depression Inventory (BDI) ( Beck et al., 1961 ) at each testing point. After completion of the study, patients were offered treatment openly with escitalopram according to usual clinical practice. Clinical response to the SSRI was defined as a reduction in HAMD of 50% or more from baseline after six weeks of treatment.

Data acquisition
fMRI data were acquired on a 3T Siemens Magnetom TIM TRIO scanner (Siemens AG) at the University of Oxford center for Clinical Magnetic Resonance Research (OCMR). A total of 180 vol of resting-state fMRI (rs-fMRI) were acquired with a voxel resolution of 3 × 3 × 3.5 mm, repetition time (TR)/echo time (TE)/flip angle (FA) = 2000 ms/28 ms/89°The experiment lasted 6.04 min. While averaging over long time-series might intuitively provide a better estimator of real connectivity, long time-series scan may be changeling for participants especially for patients with depression. It has been shown that a relatively small number of data points (ranging from 180 to 512 data points) are enough to measure sufficient BOLD activity for identifying resting state brain networks ( Cole et al., 2010 ;Braun et al., 2012 ). Therefore, if network metrics derived from short time-series provide equally reliable results as it was shown for conventional resting-state analysis ( Van Dijk et al., 2010 ), it could also reduce the burden on patients. T1-weighted structural images were acquired using a magnetization prepared rapid acquisition by gradient-echo sequence (MPRAGE) with a voxel resolution 1.0 × 1.0 × 1.0 mm on a 208 ×256 ×200 grid, TR/TE/inversion time (TI) = 2040 ms/4.68 ms/900 ms. Gradient echo phase and magnitude field maps to correct distortion were also acquired with a voxel resolution of 3.5 × 3.5 × 3.0 mm, TR/TE1/TE2/FA = 488 ms/5.19 ms/7.65 ms/60°Resting-state data acquisition was performed with participants having their eyes closed.

Resting state fMRI data processing
Resting state fMRI data preprocessing was performed in Statistical Parametric Mapping 12 (SPM12)-based DPABI software (http://rfmri.org/dpabi) running in MATLAB R2017a ( Brett et al., 2002 ;Yan et al., 2016 ). The first ten time points were discarded to allow magnetization equilibrium. Preprocessing of the image volumes included slice-time correction, head motion correction, spatial normalization to Montreal Neurological Institute (MNI) space with a resampling resolution of 3 × 3 × 3 mm, spatial smoothing, and linear trend removal. All images were then filtered using a typical temporal bandpass filter (0.01 -0.08 T) to reduce the low-frequency drift, physiological respiratory artifact, and cardiac noise. The Desikan-Killiany-Tourville (DKT) atlas was used for the parcellation of the brain into 62 regions of interest (ROIs) ( Klein and Tourville, 2012 ).

Functional brain network controllability analysis
The schematic of the methodology is shown in Fig. 1 . For details of parameters estimation for the network dynamic model, refer to ( Deng and Gu, 2020 ;M. Gilson et al., 2016 ). Briefly, fBNC analysis was first performed by fitting an advanced network dynamic model, where the rate of state changes were determined by the current state and a random diffusion process ( M. Gilson et al., 2016 ). Mathematically, the network dynamic model can be described as: where is the brain state of region , is the time constant indicating state decay over time, represents the effective connectivity between region and , and ( * ) is formally a Wiener process with covariance ( Stirzaker, 2005 ). To fit the dynamics, the unknown parameters , , and were estimated by minimizing the loss between the modelderived and empirical covariance matrices. Let = −1∕ + ∑ ≠ be the Jacobian matrix of Eq. (1) ( Chen et al., 2017 ). Lyapunov optimization ( Zheng and Cai, 2014 ), a machine learning-based approach, was applied to estimate the unknown parameters and in Eq. (1) . Analogous to the classic control equation ( Gu et al., 2015 ;Friedland, 2012 ): where is the brain connectivity matrix with dimension x ( is the number of brain regions), is the x input matrix, and is the x 1 external stimulation, the fitted dynamics of Eq.
(1) will be reformulated as: Based on the classic control Eq.
(2) and the fitted dynamics in Eq.
(3) , the linear control system was built by setting up the transition matrix as the estimated Jacobian matrix ( ̂ , x ), formating the control matrix as the estimated covariance matrix of the noise structure ( ̂ * , x ) , and configuring the external stimulation as ( ) with dimension x 1, following the model construction in ( Deng and Gu, 2020 ).

Average controllability
The average controllability reflects the capability of a node to steer the network system into easy-to-reach states with minimal effort, i.e., input energy. ( Gu et al., 2015 ). It was calculated as the 2 -norm of the network system ( Bunse-Gerstner et al., 2010 ). Mathematically, it was defined as: Cognitively, the brain areas with high average controllability are important in allowing the brain to move efficiently between different cognitive states that require little cognitive effort ( Gu et al., 2015 ).

Modal controllability
The modal controllability reflects the ease of a node to push the brain network system into difficult-to-reach states ( Gu et al., 2015 ). Mathematically, it was defined as: is the element of the eigenvectors matrix of and λ is the th eigenvalue. From a cognitive perspective, the brain areas with high modal controllability are important in switching between different cognitive states requiring significant cognitive effort ( Gu et al., 2015 ). using Desikan-Killiany-Tourville atlas. The resting state fMRI signals were extracted from each ROI and utilized to construct the effective brain network and calculate the functional controllability of each ROI. (B) Antidepressant treatment periods. The resting state fMRI signals were recorded at baseline (before treatment), one week, and six weeks post-treatment. (C) Schematic of the optimal control paradigm. In the optimal control design, the initial brain state (0) has some position in space that evolves over time toward a predefined target state ( ) . At every time point, the optimal energy ( ( ) ) required at the stimulating site to propel the system to the target state was calculated.

Personalization of the rTMS targets via fBNC
Average and modal controllability are two systematic diagnostics that characterize the ease of driving the brain into new states as time evolves ( Gu et al., 2015 ). We focused on these two controllability measures to determine and personalize the optimal rTMS intervention target to steer the resting brain state of depressed patients from visit 1 to visit 3. The resting brain state was defined as the z-scored mean of rs-fMRI signal amplitudes over times at each region of interest (ROI) ( Retter et al., 2021 ;Weaver et al., 2016 ;Larsson and Smith, 2012 ). As reported in ( Schalk et al., 2017 ), the physiological signal amplitude is better predictor of cortical excitability and inhibition than either oscillatory power or phase employed in other study ( Stisoet al., 2019 ), which is more appropriate for detecting the stimulation effects of rTMS. The resting brain state of patients at pre-treatment was defined as the initial state (visit 1) with the highest clinical assessment score for depression (HAMD scores: 23.3 ± 4.7 (mean ± standard deviation)). The resting brain state after the antidepressant treatment is defined as the target state (visit 3) with the lowest clinical assessment score for depression (HAMD scores: 9.5 ± 9.0 (mean ± standard deviation)). To decide which controllability measure was relatively more critical for the resting-to-resting state's transition, we calculated the control energy required for each node to steer each patient to the target state. The node with less control energy demand indicates that it is more capable and efficient in steering this specific brain state transition; thus, it is chosen as the optimal stimulation target. Comparisons were made for the energy required for the nodes with the maximum average controllability, maximum modal controllability, and the left DLPFC. The left DLPFC was included to compare the energy consumption between our selected optimal target node and the standard rTMS stimulation target. Once the controllability diagnostic was determined (the one with less control energy required to steer the resting brain states transition), for each patient, the node with the maximum controllability value was adopted as the optimal target region for the subsequent optimal control simulations. The energy calculation for the brain states transition is introduced in Section 2.5 .
To compare the stimulation effects by single node vs. multiple nodes targeting, we also determined the driver nodes on top of the brain network for each patient. Driver nodes are the nodes in the network that can provide complete control over the network system once the appropriate external inputs are exerted on these driver nodes. We used the maximum matching algorithm to find the driver nodes from the functional brain networks ( Liu et al., 2011 ;Tahmassebi et al., 2018 ;Leitold et al., 2019 ). Briefly, given a graph = ( , ) , where is the vertices and represents the edges, a matching in is a set of pairwise nonadjacent edges, none of which are loops, which means no two edges share a common vertex. A maximum matching is a matching that contains the largest possible number of edges ( Glover, 1967 ).

Personalization of the rTMS parameters via the optimal control analysis
Upon identifying the optimal rTMS target in Section 2.3 , an optimal control analysis was then employed to iteratively identify the optimal control energy required to steer the resting brain state at visit 1 to visit 3 of depression patients. The details of parameters estimation of optimal control process can be found in ( Stisoet al., 2019 ). Briefly, the optimization problem was first defined as: where is the target state, is the control horizon, a free parameter that defines the finite amount of time to reach the target state, and is a free parameter that weights the input constraint.
is the set of nodes to constrain. With these definitions, two constraints emerged from the optimization problem. First, ( − ( ) ) ′ ( − ( ) ) constrains the trajectories of a subset of nodes by preventing the system from traveling too far from the target state. Second, ( ) ′ constrains the amount of input utilized to reach the target state, a requirement for biological systems, which are limited by metabolic demands and tissue sensitivities.
Once the optimal control energy was calculated for each patient, it was then applied to the designated target regions to steer the initial resting brain state of the depression patient at visit 1 to the simulated final state by performing the open-loop control ( Stisoet al., 2019 ). Then, the maximum model correlation between the simulated final state and the actual final state (the resting brain state at visit 3) was computed to quantify the stimulation effects of the targeting strategies. The maximum model fitting accuracy was computed by selecting the control horizon from 0 to 0.75, with an increased step of 0.05 to the (so = 0.05, 0.10, 0.15, 0.20, …), and calculating the maximum model correlation within this period. The total energy required for the state's transition was calculated. The measure of total energy can be described as: To clarify, in this study, the term "parameter " refers to the control energy to steer the brain states transition, but does not include the stimulation targets. The control energy can further be decomposed into different parameters such as stimulation intensities, frequencies, and stimulation durations for the external stimulation devices in clinical practice.

Performance evaluation
The performance of the proposed personalized rTMS technique was evaluated by comparing it to the standard rTMS protocol in treating depression in terms of the treatment efficacy. Three treatment strategies were compared: rTMS with personalized stimulation node, rTMS with multiple driver nodes, and the standard rTMS with the lDLPFC as the stimulation target. The stimulation parameters in the first two strategies were optimized as described in Section 2.5 to personalize the rTMS treatment. The stimulation parameters in the last strategy, i.e., standard rTMS, were fixed with the parameters reported in the literatures ( Khokhar et al., 2021 ;McClintock et al., 2017 ). The treatment efficacy was evaluated by the maximum model fitting accuracy achieved by exerting the stimulation parameters to the initial state (visit 1) of each depression patient and then computing the maximum model correlation between the simulated final state and the actual final state (visit 3). The first treatment strategy was to design the stimulating region as the optimal target by choosing either of the controllability diagnostics (average controllability or modal controllability). Then the next step was to exert the optimal stimulation parameters inferred from the optimal control algorithm to the node to guide the resting brain states transition. The second treatment strategy was similar procedure to the first, but setting the driver nodes as the stimulation targets. The last strategy was following the standard rTMS protocol by selecting the lDLPFC as the stimulation target and exerting the empirical parameters on the target node. The empirical stimulation amplitude in Tesla was in the range from 1.5 T to 2 T, the stimulation frequency in Hertz was between 1 Hz and 20 Hz, and the stimulation duration was the number of sample points (126 in this study) divided by 1000 ( Stisoet al., 2019 ;Khokhar et al., 2021 ;McClintock et al., 2017 ). The control energy exerted on the lDLPFC for standard rTMS treatment was then calculated as log(frequency) x amplitude x duration, as outlined in ( Stisoet al., 2019 ).

Statistical analysis
The average and modal controllability of DMN and FPN were compared in healthy controls using Wilcoxon signed-rank test ( Woolson, 2007 ). The energy required for the node with the maximum 1.4 ± 1.20 18.6 ± 2.6 * * Significant difference between depression patients and healthy subjects ( P < 0.05).
average controllability, the node with the maximum modal controllability, and the standard rTMS target, lDLPFC, to transition the resting brain state at visit 1 to visit 3, were also statistically compared. The Friedman one-way ANOVA was first utilized to perform the group comparisons and then the Wilcoxon signed-rank test was employed for group pairwise comparisons ( Zimmerman and Zumbo, 1993 ). A linear regression model was employed to identify the correlation between the simulated final state and actual final state (resting brain state at visit 3) to reflect the model fitting accuracy ( Aalen, 1989 ). Then, the maximum model fitting accuracy achieved by the three targeting strategies (targeting the node with the maximum average controllability, the driver nodes, and the lDLPFC) were statistically compared. The Friedman one-way ANOVA was also applied to perform group comparisons and then the Wilcoxon signed-rank test was employed to compare each two groups. Clinical scores were compared with healthy controls by independent t-test and between different visits using paired t tests ( De Winter, 2013 ). Outliers, defined by mean plus/minus three standard deviations, were excluded ( Howell et al., 1998 ). A critical P-value of 0.05 was used and the false discovery rate (FDR) method was employed to correct for multiple comparisons ( Laird et al., 2005 ). In this study, the term of "group comparisons " refers to the comparisons between different targeting strategies, instead of the comparison between groups of depression patients and healthy subjects.

Demographic and clinical behavior data
Five patients were excluded from the analysis due to the incompletion of the antidepressant experimental imaging study in depression. Two patients decided not to continue beyond the first scan and three patients had two scans but could not attend the final appointment. Nine healthy controls and nine patients were excluded from the analysis due to the poor quality of resting state fMRI data contaminated by large motion artifact recorded at the early stage of the study. Five patients were excluded due to the inconsistent and lesser resting state fMRI recording time compared to the other patients, which may induce higher variability and lower reliability of the functional brain network estimation and the network dynamic model fitting in this study ( M. Gilson et al., 2016 ;Birn et al., 2013 ). Table 1 summarizes the demographic information and clinical behavioral scores of the remaining 20 patients with depression and 23 healthy controls, including gender, age at onset of depression patients, age at the baseline appointment, and baseline clinical scores of HAMD and BDI and QIDS. Statistical analysis showed that there were no significant differences between healthy controls and depression patients in terms of gender ( P = 0.72, 2 ( 1 , = 43 ) = 0.124, chi-square test) and age ( P = 0.43, t -value = − 0.8, df = 41). As expected, clinical scores of healthy controls were significantly lower than depression patients in terms of the HAMD ( P = 1.43 x 10 −25 , t -value = 23.3, df = 41), BDI ( P = 8.93 x 10 −27 , t -value = 25.1, df = 41) and QIDS ( P = 3.60 x 10 −28 , t -value = 27.2, df = 41) scores. Fig. 2. Controllability of cognitive systems in healthy controls. The average and modal controllability on the two cognitive systems, DMN, and FPN, were plotted and compared. The values on the y-axes represent the mean of controllability ranks. " * " indicates significant difference.

Fig. 3.
Energy required for the nodes with the maximum average controllability (AC_max), maximum modal controllability (MC_max), and standard rTMS target in lDLPFC to steer the resting brain states transition. The values on the y-axes represent the energy ranked over all brain regions. So larger ranks indicated more energy required for the specific node to transit the depression patients' resting brain state at visit 1 to visit 3. " * " indicates significant difference.

Personalized rTMS target(s)
The average and modal controllability diagnostics were characterized by the fBNC method and served as potential indices to determine the optimal stimulation target(s). Based on the physiological properties of DMN and FPN, we assumed that the DMN controls the "easy-to-reach " states transition and, thus, had relatively higher average controllability than the FPN ( Gu et al., 2015 ). Conversely, the FPN controlled the "difficult-to-reach " state, thus, had relatively higher modal controllability than the DMN ( Gu et al., 2015 ). The average and modal controllability of DMN and FPN in healthy controls were calculated, respectively, and compared in Fig. 2 . The results showed a significantly higher average controllability in the DMN compared to the FPN ( P = 0.007, zvalue = 2.447), and a significantly higher modal controllability in the FPN compared to the DMN ( P = 0.008, z -value = − 2.410). This finding is consistent with the findings in a previous study ( Gu et al., 2015 ).
We then computed the energy required for the node with the maximum average and modal controllability to steer the resting brain state from visit 1 to visit 3. Fig. 3 shows the energy needed for the node with the maximum average controllability, maximum modal controllability, and the lDLPFC. The results show that the mean energy rank required for the node with the maximum average controllability was 1.81, which was almost the minimum rank across all 62 brain regions (ranked from 1 to 62, higher ranks represent more energy required). These findings indicate that the least energy necessitated for the node with the maximum average controllability for this specific resting-toresting brain states transition. The Friedman one-way ANOVA revealed significant group differences in the energy consumption between the three targeting strategies ( 2 (2) = 27.11, P = 1.29 x 10 −6 ). Meanwhile, the mean energy rank required for the node with the maximum average controllability was significantly lower than the node with the maximum modal controllability ( P = 1.10 x 10 −3 , z -value = − 3.704) and the lDLPFC ( P = 1.60 x 10 −3 , z -value = − 3.598), suggesting that the node with the maximum average controllability was more efficient than the node with the maximum modal controllability as well as the standard rTMS target in lDLPFC in being able to steer the resting-to-resting state's transition of depression. It further indicated that the resting-to-resting state's transition of depression patients was an "easy-to-reach " process; thus, the node with the maximum average controllability was the node that was most optimal. Therefore, the ROI with the maximum average controllability was selected as the optimized stimulation target in the development of the proposed personalized TMS technique since it required much less energy than the other nodes, meaning it was easier to modulate the resting brain states transition from depression (visit 1) to 'healthy' (visit 3, post-treatment). The node with the maximum average controllability of each patient was selected as the optimal target for the following simulation and comparisons with performance of multiple driver nodes. Fig. 4 shows the results of the optimal control analysis from three typical patients (HAMD scores of pre vs. post are Patient #14: 28 vs. 1, Patient # 13: 13 vs. 1, and Patient # 05: 31 vs. 8, respectively). The node with the maximum average controllability for each patient was selected as the stimulation target. The optimal control analysis was applied to estimate the optimal control energy required for this node to steer the brain state from initial depression (visit 1) to final healthy (visit 3). The initial state was defined as the resting brain state of patients at visit 1 (pre-treatment), while the final state was defined as the resting brain state at visit 3 (post-treatment) as shown in Fig. 4 A. Fig. 4 B shows the fitted curves for the simulated final state and the actual final state; Fig. 4 C demonstrates the correlation between the Simulation and the Observation results. The Simulation represents the simulated final state predicted from the initial state by delivering the estimated optimal control energy to the stimulation target. The actual final state (Observation) is defined as the resting brain state characterized from resting state fMRI data collected at visit 3. Results showed that the observation and the simulation curves were well-fitted and all significantly correlated (Patient #14: r = 0.50, P = 3.80 x 10 −5 , Patient # 13: r = 0.51, P = 2.00 x 10 −5 , Patient # 05: r = 0.32, P = 0.012), indicating the theoretical feasibility of utilizing the proposed personalized rTMS technique to efficiently steer the resting brain states transition from depression at visit 1 to healthy at visit 3.

Comparison between the standard and personalized rTMS treatment protocol
We compared the performance of the proposed personalized rTMS protocol and the standard rTMS protocol in terms of the maximum model fitting accuracy, as shown in Fig. 5 . The stimulation target was personalized as the node with maximum average controllability ( Section 2.4 ) and the rTMS parameters were personalized using the optimal control analysis ( Section 2.5 ). The lDLPFC was used as the stimulation target and the empirical stimulation parameters were applied in the standard rTMS protocol ( Section 2.5 ). The Friedman one-way ANOVA demonstrated significant group differences of the maximum model fitting accuracies achieved by the three targeting strategies ( 2 (2) = 6.26, P = 0.044). From the results in Fig. 5 , we can see that the maximum model fitting accuracies achieved by targeting the node with the maximum average controllability or the driver nodes were not significantly different ( P = 0.18, z -value = 0.915), indicating the potential to obtain the comparable level of model fitting accuracy by targeting either the optimal single target or multiple driver nodes. Nevertheless, the standard rTMS targeting strategy showed significantly lower maximum model fitting accuracy than both personalized rTMS protocols by targeting either the node with the maximum average controllabil- Fig. 4. Optimal control simulation of three typical patients. The node with the maximum average controllability was selected as the target node. (A) The initial states (blue lines) and the final states (red lines). The initial states were defined as the depression patients' resting brain state at visit 1 and the final states were defined as the depression patients' resting brain state at visit 3. (B) Model fitting. The red lines represent the actual final state (resting brain state at visit 3, also called Observation) and the dark blue lines represent the simulated final state (Simulation). (C) Regression analysis between the Observation and the Simulation results to reflect the model fitting accuracy. ity ( P = 0.018, z -value = 2.090) or the driver nodes ( P = 0.0038, zvalue = 2.670).
As shown in Fig. 6 , our results indicated that the main brain areas with high average controllability at visit 1 included the left temporal cortex, left precuneus, left/right ventrolateral prefrontal cortex (VLPFC), left precentral cortex, left/right ventromedial prefrontal cortex (VMPFC), left/right anterior cingulate cortex (ACC), right inferior parietal cortex, and right fusiform. The brain regions with high likelihood to be driver nodes were mainly located in the left/right dorsolateral prefrontal cortex, left temporal cortex, left caudal middle frontal cortex, left anterior cingulate cortex, left parahippocampal cortex, right posterior cingulate cortex, and right lingual cortex.

Discussion
While brain neuromodulation has shown great therapeutic potential, such as the FDA-approved rTMS intervention treatment for depression ( Schalk et al., 2017 ), its optimization and personalization remain challenging. This is partly due to limited knowledge of the underlying neural control mechanisms and the accommodation of inherent interpatient variability and intersession alterations during treatments in depression ( Cole et al., 2021 ;Liu et al., 2011 ). This study proposed a personalized rTMS intervention protocol to treat patients with depression by integrating functional brain network controllability analysis and optimal control analysis. By employing a novel fBNC analysis strategy ( Medaglia et al.,  Fig. 5. Comparison of the maximum model fitting accuracies achieved by three targeting strategies (stimulating driver nodes, node with the maximum average controllability, and the standard rTMS target in lDLPFC), respectively. 2017 ; Muldoon et al., 2016 ), we adopted a control model for exact brain states transition in which the optimal control energy demanded the transition of the state is minimized ( Modak and Fitzgerald, 2021 ). Our results suggested that the personalized rTMS protocol was more efficient than the standard rTMS treatment in steering the resting brain states transition of depression patients indicated by the significantly lower control energy required and higher maximum model fitting accuracy achieved by the personalized protocol. We also demonstrated that the node with the maximum average controllability outperformed the other nodes in steering the resting brain state of depression patients at an early diseased stage (visit 1) to the later healthier stage (visit 3) with less control energy needed. These findings indicate the preference for selecting the node with the maximum average controllability for the transition of the resting-to-resting state on depression patients. Finally, our results further illustrated that by targeting the single node with the maximum average controllability, the model fitting accuracy could reach the level of multiple driver nodes targeting, indicating the theoretical feasibility of minimizing the number of stimulation targets while not compromising efficacy in clinical application.
Network control theory predicts the controllability of large-scale neural circuits ( Gu et al., 2015 ). Conventional controllability analysis on structural brain networks predicts the ability to alter large-scale neural circuitry, assuming the structural brain networks can modulate the transition between brain states ( Gu et al., 2015 ;Beynel et al., 2020 ). However, recent advances in physics and engineering have promoted control system analysis and switched the focus from structural architecture to functional dynamics ( Deng and Gu, 2020 ;Baggio et al., 2021 ). Consequently, instead of building upon the limited models based on structural brain network that ignore differences between functional states and the short-term alterations, it is critical and clinically important to develop and optimize theoretical control analysis based on functional brain networks ( Deng and Gu, 2020 ). In practice, it is also more attainable, cost efficient, and clinically practicable to build a control framework from the functional signals constructed from the functional imaging modalities such as electroencephalogram (EEG), magnetoencephalography (MEG), and the blood-oxygen-level-dependent (BOLD) signals ( Verplaetse et al., 2016 ;Foldes et al., 2015 ;Li et al., 2018 ). The current work employed the functional controllability strategy to evaluate a set of nodes' ability in moving the brain into disparate states, inferring the brain states' transition and co-varying patterns in the functional signal space ( M. Gilson et al., 2016 ). Previous studies have shown that the signal space inference can reveal the intrinsic interactions between different regions that mediate the underlying control mechanism, indicating the feasibility of considering the functional controllability in the signal space ( Deng and Gu, 2020 ;Baggio et al., 2021 ).
Average and modal controllability are two controllability diagnostics utilized to evaluate the ease of the brain to enter different cognitive states ( Gu et al., 2015 ). Average controllability is posited to quantify the nodes' roles in moving the brain system into many easily reachable states, while the modal controllability measure defines a node's ability to steer the brain system into many difficult-to-reach states ( Gu et al., 2015 ). If control energy can be likened to cognitive effort and if brain states can be likened to cognitive functions, then the easy-to-reach state refers to the brain state that requires little cognitive effort to reach. At the same time, difficult-to-reach state refers to the brain state that requires significant cognitive effort to reach from the initial brain cognitive state ( Gu et al., 2015 ). The easy-to-reach and difficult-to-reach states have been explored using the brain controllability before ( Gu et al., 2015 ), however, they have not been applied to depression brain studies. To the best of our knowledge, this study represents the first effort to employ the brain controllability to study depressed brains, and further developed into a personalized neuromodulation technique for precision depression treatment. We showed that the DMN has significantly higher average controllability than the FPN, while the FPN has significantly higher modal controllability than the DMN cognitive system in healthy controls ( Fig. 2 ). Prior work has suggested that the DMN is active when a person is not focused on the outside world and the brain is at wakeful rest, such as during daydreaming and mind-wandering ( Kumar, 2019 ). In contrast, FPN is thought to support distinct functional roles such as task-switching and task-set maintenance ( Marek and Dosenbach, 2018 ). These findings indicate that the DMN is usually involved in transiting the brain network into many effortless brain states, with low energy requirements. In contrast, the FPN cognitive system is more crucial in steering the brain system into cognitively demanding states associated with active engagement in complex tasks. These interpretations are consistent with the underlying neural control driver nodes distribution averaged over all patients at baseline (visit 1). In the average controllability distribution maps, the value of each brain region represents the z-score of the average controllability value. In the driver nodes distribution maps, the value of each brain region indicates the z-score of the node's possibility to be driver node. The z-scores are thresholded to be above 0 so only the main brain regions are shown in the distribution maps. mechanism explained by the two cognitive systems' available average and modal controllability diagnostics, linking the engineering control concepts to the psychological construct of cognitive control.
Upon corroborating the control diagnostics' properties, by correlating the engineering control concepts with the brain cognitive control systems, we further assessed the ability of various nodes to steer a specific brain state transition, from the resting brain state at the early depressive stage to the later healthier stage, directly relating the control theory to specific psychiatric disease treatment. This portion of the investigation was made possible by an important modeling advance addressing the challenge of stimulating a trajectory whose control is dominated by a single node ( Klooster et al., 2022 ;Modak and Fitzgerald, 2021 ;Medaglia et al., 2017 ). Our results revealed that the average controllability played a relatively more important role in determining the target region in steering the resting-to-resting functional states transitions between different treatment sessions and preferred to be designated as the optimal target node. This conclusion is inferred from our results that the node with the maximum average controllability required the least control energy across all 62 brain regions for this specific resting-to-resting states transition of depression and was much less than the node with the maximum model controllability and the standard rTMS target in lDLPFC ( Fig. 3 ). The significantly lower control energy required for the node with the maximum average controllability than the node with the maximum model controllability further supported our assumption that the resting-to-resting state's transition in depression was an easy-to-reach process. Thus, the node with the maximum average controllability, which is the node most capable of steering the easily reachable states with less control energy needed ( Gu et al., 2015 ), demanded less control energy than the node with maximum model controllability and other nodes. Together, our results indicate the potential to utilize our proposed fBNC-based optimal control approach to determine the optimal stimulation target for specific brain state transitions based on the functional average and modal controllability diagnostics.
Prior studies have demonstrated a reduced blood flow and metabolism in the left DLPFC in patients suffering from major depression disorder, suggesting the hypoactivity of the left DLPFC tends to cause depressive symptoms ( Rajkowska et al., 1999 ). For this reason, the stimulation is applied to the left DLPFC for standard rTMS treatment to increase the releasing of neurotransmitters and local brain activity of depression patients ( Baeken and De Raedt, 2011 ). However, rTMS treatment outcomes are still highly variable for depression treatment. Early studies in a large cohort of data showed that about 30% of patients responded to the rTMS treatment, with 18.6% achieving remission ( Berlim et al., 2014 ). The variability caused by the standard rTMS treatment may be due to a lack of understanding of the underlying neural control mechanisms and the heterogeneity of stimulation montage for each depression patient undergoing rTMS treatment. To overcome the limitation of the heterogeneity in rTMS treatment protocols, we proposed a personalized neuromodulation protocol in this study. We first performed the fBNC analysis to locate the optimal target region for each patient, and then applied the optimal control algorithm to obtain the optimal stimulation parameters demanded for the optimal target to steer the specific brain state transitions. Our results revealed that the personalized rTMS protocol showed higher efficiency than the standard rTMS treatment as implied by the significantly lower control energy required and higher maximum model fitting accuracy achieved by the personalized protocol ( Fig. 5 ). Furthermore, our results also demonstrated the theoretical feasibility of achieving comparable model fitting accuracy by targeting a single node with the maximum average controllability or multiple driver nodes, evidenced by the insignificant differences in the maximum model fitting accuracies achieved by the two targeting strategies, respectively. Our findings may indicate the potential to propose a more straightforward and more practical neuromodulation protocol by targeting the node with the maximum average controllability for depression intervention treatments.
Although the current study proposed a personalized rTMS treatment protocol, which requires personalizing both stimulation targets and parameters at each visit for each patient, our results also revealed some critical brain regions that have a higher chance for being selected as the optimal stimulation target at the early stages of depression, such as the ACC, VLPFC, and VMPFC. Previous studies have demonstrated these brain regions as potential rTMS stimulation targets for depression treatment ( Downar and Daskalakis, 2013 ;Kreuzer et al., 2015 ;Tan et al., 2020 ;MO et al., 2021 ). For example, a prior study that employed rTMS with double cone coil revealed that the stimulation of the ACC was significantly superior to conventional rTMS stimulation at the end of the 3 week treatment period ( Kreuzer et al., 2015 ). The VLPFC has been reported to draw input from external sensory modalities and control the other brain regions such as the basolateral amygdala, hypothalamus, ventral striatum, and other limbic structures to mediate the success of cognitive arousal in depression ( Downar and Daskalakis, 2013 ). A recent rTMS stimulation study has further demonstrated that the stimulation to the VLPFC could effectively improve the ability to explicitly regulate cognition in depression patients ( MO et al., 2021 ). A previous study also found that rTMS targeting the VMPFC could alleviate depression symptoms and was more sensitive in patients with severe emotional symptoms ( Tan et al., 2020 ). Overall, the critical stimulation targets in this study, combined with the functional brain controllability measure, may help us further interpret rTMS treatments' effectiveness and provide some promising alternative stimulation targets for future rTMS treatments.
To the best of our knowledge, this study represents the first effort to personalize rTMS technique in treating depression by personalizing the stimulation target using fBNC and personalizing the stimulation parameters using the optimal control analysis. It is worth noting that the current work still has some specific limitations. The nonlinear effect still exists in our model and remains to be solved in the future. Moreover, due to the constraints of unpredictable noise in the measurement, the approximation to the observed trajectories is still unsatisfactory, especially for the actual values ( Deng and Gu, 2020 ). Meanwhile, in this study, we utilized the DKT atlas with 62 regions of interest, a successor to the DK adult cortical parcellated atlas with 68 regions of interest employed in the previous study ( M. Gilson et al., 2016 ). Some regions with unclear or arbitrary boundaries were removed, and many existing boundaries were revised to conform to sulcal fundi in the DKT atlas. This provides greater anatomical consistency across individuals due to clearer and more reproducible landmarks. However, we should notice that a more fine-grained parcellation scheme that respects both functional and structural features of the brain (e.g., the Glasser, Schaefer, or Brainnetome atlas) may be more appropriate for future improvement of the studies ( Schaefer et al., 2018 ). By utilizing a more fine-grained parcellation scheme, we may reduce the variability within each single brain region, such as the variability within the DLPFC. Finally, the stimulation effect of the standard rTMS treatment protocol was evaluated by fitting the empirical parameters to each patient and calculating the maximum model fitting accuracy. While the parameters are commonly used in clinical rTMS treatment ( Khokhar et al., 2021 ;McClintock et al., 2017 ), they do not strictly follow the standard protocol by determining the stimulation amplitude based on the motor-evoked potentials (MEP) ( Fitzgerald et al., 2006 ). Immediate next steps will include testing the validity of our personalized neuromodulation protocol with rTMS intervention treatment on depression patients and other psychiatric diseases.
In addition, as reported in ( Fitzgerald, 2021 ), there are mainly two issues that need to be solved in developing an optimal rTMS therapy. The first issue is how to specifically delineate the optimal brain target for the best neuromodulation efficacy, and the second one is how to position the rTMS coil precisely over the identified stimulation brain target in practice. In our study, we aim to solve the first issue, by delineating the optimal cortical target for stimulation. We went even further by optimizing stimulation parameters for the optimized stimulation target(s). The second issue on how to position the rTMS coil in practice beyond the scope of this current study and will be deferred to a future TMS navigation study. Meanwhile, conventional rTMS protocol can only directly stimulate superficial cortex effectively, with an approximately 3 cm maximum depth from the surface of the head ( Deng et al., 2013 ). Recent studies have demonstrated the downstream effects from rTMS stimulation on deep brain structures and showed the clinical efficacy of indirectly targeting sub-cortical brain regions (like subgenual anterior cingulate cortex (sgACC)) by stimulating the cortical brain areas (like DLPFC) ( R.F. Cash et al., 2021 ;Cash et al., 2019 ;Luber et al., 2022 ). Those findings support the feasibility of utilizing our personalized rTMS technique to better modulate the cortical and sub-cortical brain areas such as the sgACC, which are known associated with the depression pathology, to further maximize the treatment efficacy. However, we should also be aware that the deep rTMS also has some flaws such as less flexibility for personalized treatment due to the shape of the stimulator ( Lonergan et al., 2018 ).
In our study, we optimized both the stimulation targets and parameters via the brain controllability analysis and optimal control analysis in order to most efficiently complete a specific brain states transition. Compared to the recent SNT trials and other personalized rTMS protocols that limited their target optimization to the DLPFC only ( Cole et al., 2021 ;Fitzgerald, 2021 ;Klooster et al., 2022 ;Modak and Fitzgerald, 2021 ), our study optimized the stimulation target(s) within the whole brain network to avoid local optimization limitation. Besides the stimulation target optimization, we also further optimized the stimulation parameters (such as control energy, intensity, frequency, duration). Therefore, our approach allows for a global optimization for both the stimulation target(s) and the corresponding stimulation parameters, which can most efficiently steer a depressive brain state to the healthy brain state for precision treatment. The personalization framework developed in this study can be easily transferred to other neuromodulation modalities such as direct electrical stimulation, deep brain stimulation and other neurmodulation techniques for treating other neurological diseases and psychiatric disorders including Stroke, Parkinson' disease, Alzheimer's disease, Depression, Anxiety disorder, and others. The fBNC assessment and monitoring techniques developed in this study can also be used to understand the controlling mechanisms of the brain in disease states or underlying neuromodulation or neurorehabilitation treatments.

Conclusion
Our study demonstrated the preliminary feasibility of personalizing the rTMS technique in treating depression with optimized treatment efficiency and outcome. The stimulation target and stimulation parameters were personalized for specific patients using the fBNC and optimal control analysis. The developed personalized rTMS was evaluated by comparing against the standard rTMS in terms of the model fitting accuracy and correlation with clinical assessment. Our results showed that the developed personalized rTMS protocol was more efficient than the standard rTMS treatment. In conclusion, our work offers empirical support for integrating functional brain networks controllability and optimal control approach in designing a personalized neuromodulation protocol and interpreting the psychiatric diseases' underlying neural control mechanisms like depression. The personalization framework developed in this study can be transferred to other neuromodulation modalities for treating other neurological diseases and psychiatric disorders.

Data and code availability statement
The datasets used in this study and the code developed in this study are available per reasonable request.
FF developed the algorithms, conducted data analysis, draft the manuscript and intemperate the results. BRG conducted the antidepressant clinical trial, collected resting state fMRI data and performed clinical assessment. RYC helped with the result interpretation and manuscript revision. SIS helped with the result interpretation and manuscript revision. SS helped with the antidepressant clinical trial design, resting state fMRI data collection, clinical assessment, also helped with the result interpretation and manuscript revision. YZ initiated the idea, initiated the project, supervised the study design, algorithm development and data analysis, and also helped with result interpretation and manuscript revision.

Declaration of competing interest
None