Variability in cTBS Aftereffects Attributed to the Interaction of Stimulus Intensity With BDNF Val66Met Polymorphism

Objective: To evaluate whether a common polymorphism (Val66Met) in the gene for brain-derived neurotrophic factor (BDNF)—a gene thought to influence plasticity—contributes to inter-individual variability in responses to continuous theta-burst stimulation (cTBS), and explore whether variability in stimulation-induced plasticity among Val66Met carriers relates to differences in stimulation intensity (SI) used to probe plasticity. Methods: Motor evoked potentials (MEPs) were collected from 33 healthy individuals (11 Val66Met) prior to cTBS (baseline) and in 10 min intervals immediately following cTBS for a total of 30 min post-cTBS (0 min post-cTBS, 10 min post-cTBS, 20 min post cTBS, and 30 min post-cTBS) of the left primary motor cortex. Analyses assessed changes in cortical excitability as a function of BDNF (Val66Val vs. Val66Met) and SI. Results: For both BDNF groups, MEP-suppression from baseline to post-cTBS time points decreased as a function of increasing SI. However, the effect of SI on MEPs was more pronounced for Val66Met vs. Val66Val carriers, whereby individuals probed with higher vs. lower SIs resulted in paradoxical cTBS aftereffects (MEP-facilitation), which persisted at least 30 min post-cTBS administration. Conclusions: cTBS aftereffects among BDNF Met allele carriers are more variable depending on the SI used to probe cortical excitability when compared to homozygous Val allele carriers, which could, to some extent, account for the inconsistency of previously reported cTBS effects. Significance: These data provide insight into the sources of cTBS response variability, which can inform how best to stratify and optimize its use in investigational and clinical contexts.


INTRODUCTION
Transcranial magnetic stimulation (TMS) has received considerable attention in both research and clinical settings due to its ability to probe and transiently modulate cortical activity. In clinical contexts, TMS is routinely employed as a treatment for depression (Connolly et al., 2012), and recently as a treatment for migraine and obsessive-compulsive disorder (Schwedt and Vargas, 2015;Voelker, 2018) and numerous other neurologic and psychiatric conditions (Schlaepfer et al., 2010). In research settings, it has proven to be a powerful tool for interrogating brain structure-function relationships, pertaining to a wide range of motor and cognitive abilities (Devlin and Watkins, 2006;Lowe et al., 2018;Medaglia et al., 2018) as well as social Era et al., 2018Era et al., , 2020 and emotional processes (Moors et al., 2019;Fini et al., 2020), and to characterize and index fundamental neurophysiologic properties, including but not limited to cortical excitability (Pascual-Leone et al., 1998), interhemispheric interactions (Mochizuki et al., 2004), and activity-induced neuroplasticity (Bolognini et al., 2009). However, the findings from these studies are often muddled by the amount of inter-and intra-individual variability that is typically observed in response to stimulation (Ridding and Ziemann, 2010;Chung et al., 2016). Identifying and controlling for factors that contribute to this variability, both in the physiologic and behavioral aftereffects, is an important strategy for optimizing the investigative and therapeutic utility of TMS protocols.
In recent years, there is a growing interest in theta-burst stimulation (TBS), which is a modified form of repetitive TMS (rTMS). TBS reportedly interferes with cortical excitability and produces aftereffects to the same degree as the conventional rTMS protocols but in a fraction of the time; the application of TBS can take 20-190 s compared to 15-30 min of conventional rTMS protocols. While the short implementation time has made TBS a very attractive tool for research and clinical investigations, it is not immune to the observed response variability.
A TBS pattern consists of 50 Hz bursts of stimulation pulses delivered in triplets every 200 ms (at 5 Hz). Depending on the number of repetitions of the TBS pattern and the interval between repetitions, TBS can increase or decrease cortical excitability. Across studies, TBS has typically been found to be excitatory with intermittent delivery of TBS pattern (iTBS) and inhibitory with continuous TBS (cTBS) pattern (Huang et al., 2005). While the induced changes in cortical excitability with iTBS and cTBS are found at the group level, further evaluations of the data indicate that there is considerable variability in response to TBS, both within and between individuals (Ridding and Ziemann, 2010;Goetz et al., 2014;Hordacre et al., 2017;Corp et al., 2020). Evidence indicates that extrinsic factors related to stimulation parameters and experimental design such as the degree of pre-activation of targeted muscles (Iezzi et al., 2008;Goldsworthy et al., 2014), the stimulation intensity (SI) of singlepulse TMS for probing cortical excitability with TBS (Vallence et al., 2015;Goldsworthy et al., 2016b), and the SI of TBS pulse triplets (Jannati et al., 2017;Sasaki et al., 2018) contribute to this variability. In addition, properties that are intrinsic to each individual such as age (Freitas et al., 2011) and genetic factors can also contribute to the variability.
One of the most notable genetic factors that has been purported to affect learning, memory, and neuroplasticity is the brain-derived neurotrophic factor (BDNF; e.g., Cheeran et al., 2008;Jannati et al., 2017). BDNF is a protein encoded by the BDNF gene, supporting the survival, growth, and differentiation of neurons and synapses (Huang and Reichardt, 2001). BDNF plays an important role in synaptic plasticity as its release is thought to aid long-term potentiation (LTP) and long-term depression (LTD) processes (Lu, 2003). A relatively common polymorphism in the BDNF gene (Val66Met allele) is associated with a decrease in activity-dependent release of BDNF (McHughen et al., 2009) and diminished synaptic plasticity in animal models (Ninan et al., 2010). Evidence suggests that this polymorphism is also associated with impairments in learning and memory in humans (Bath and Lee, 2006;Soliman et al., 2010). Given that the persistent effects of rTMS-including TBS-are believed to be mediated by LTP-or LTD-like effects on synaptic plasticity, individuals with the Val66Met allele may be less responsive to TBS. However, evidence to support this hypothesis has been mixed. While some studies have found that Val66Met allele carriers exhibit little-to-no aftereffects of TBS when compared to their homozygous (Val66Val) counterparts (Cheeran et al., 2008;Antal et al., 2010;Lee et al., 2013;Jannati et al., 2017), others report no difference in susceptibility to TBS as a function of BDNF genotype status (Li Voti et al., 2011;Mastroeni et al., 2013). One possible explanation for this discordance between studies is that Val66Met allele carriers may be more variable than Val66Val homozygotes with respect to response to brain stimulation. A recent meta-analysis supports this view: Chung et al. (2016) found that Val66Val carriers show more consistent response across studies, with greater effect sizes for iTBS. By contrast, Val66Met carriers were more variable (also see Jannati et al., 2019).
Here, we ask whether increased variability in cTBS response among Val66Met allele carriers relates to the SI of single-pulse TMS used to probe cortical excitability before and immediately after cTBS; we refer to the single-pulses as test pulses in this study. Prior work has shown that the SI of test pulses impacts cTBS-induced suppression of motor evoked potentials (MEPs; Vallence et al., 2015;Goldsworthy et al., 2016b). This is directly relevant to studies involving TBS because different approaches have been employed by different investigators to establish SI. For example, some studies determine SI based on the resting motor threshold (rMT), while others stimulate at an SI that is empirically determined to be sufficient to generate an MEP of a certain amplitude, often 1 mV (SI 1mV ). These methodological differences in SI determination may introduce noise, especially if SI differentially affects cTBS response in Val66Val vs. Val66Met carriers (Cheeran et al., 2008;Antal et al., 2010;Lee et al., 2013;Jannati et al., 2017). In the current study, we re-examine the impact of BDNF genotype status on cTBS-induced suppression of motor excitability. We hypothesized that cTBS aftereffects for BDNF Met allele carriers would be less reliable, and investigated whether differences in SI [as determined using the percentage of maximum stimulator output (MSO) required to produce MEPs with a peak-to-peak amplitude of ∼1 mV] account for variable cTBS responses as a function of BDNF genotype status.

Overview
The experiment consisted of a single session. Participants were seated in a comfortable chair with their right arm resting on a pillow and were instructed to keep their right hand relaxed throughout the duration of the experiment. Single TMS pulses were delivered over the left primary motor cortex. We first determined the rMT, and then gradually increased SI to the percentage MSO required to produce MEPs with peak-to-peak amplitudes of approximately 1 mV (SI 1mV ). Prior to cTBS (Baseline) and in 10 min intervals immediately following cTBS for a total of 30 min post-cTBS (0 min post-cTBS, 10 min post-cTBS, 20 min post cTBS, and 30 min post-cTBS), we collected 30 MEPs from the first dorsal interosseous (FDI) of the dominant (right) hand at SI 1mV , following the recommendation from prior research (Goldsworthy et al., 2016a). Test pulses were delivered with an inter-stimulus interval of 6 s with a random jitter of ±6%. Approximately 10-15 min following baseline MEP collection and 10-15 min prior to cTBS administration, we obtained active motor threshold (aMT). Following post-cTBS MEP collection, we obtained saliva samples for BDNF genotyping (Figure 1).

Participants
Thirty-three neurologically healthy individuals (16 females) aged 18-45 [mean (M) ± standard deviation (SD) = 24.6 ± 6.2 years] with no contraindications to TMS participated in the study. All participants provided informed consent in accordance with the Institutional Review Board at the University of Pennsylvania. See Supplementary Table 1 for demographic and stimulation parameter data for each participant.

Electromyography
Electromyographic (EMG) activity was recorded using surface electrodes overlying the belly of the right FDI. Signals were amplified and band-pass-filtered between 20 and 2,000 Hz, digitized (sample-rate 5 kHz), and stored for offline analysis using SIGNAL software (Cambridge Electronic Devices, Cambridge, UK).

Transcranial Magnetic Stimulation
Single-pulse TMS with monophasic waveform was administered using a hand-held figure-of-eight coil (Magstim Company, Whitland, Dyfield, UK). The coil was positioned over the left motor cortex (M1) to a site that reliably elicited an MEP in the right FDI muscle (i.e., the motor hotspot). The Brainsight (Rogue Inc., Montreal, QC, Canada) neuronavigational system was used to mark the motor hotspot on native-space magnetic resonance image volumes collected prior to the experiment. In line with widely accepted methods, rMT was defined as the minimum pulse intensity required to elicit MEPs with peak-topeak amplitudes of at least 50 µV in 5 of 10 consecutive trials with the FDI at rest Rothwell et al., 1999;Schlaepfer et al., 2010).

Continuous Theta-Burst Stimulation
CTBS was administered with a biphasic waveform using a Magstim SuperRapid 2 Stimulator (Whitland, UK). Following the standard procedure (Huang et al., 2005), cTBS entailed continuous delivery of 50 Hz triplets of TMS pulses at 5 Hz for a total of 600 pulses (∼40 s). SI of cTBS pulses was set to 80% of aMT, defined as the minimum pulse intensity required to produce MEPs with peak-to-peak amplitudes of at least 200 µV in 5 of 10 consecutive pulses while participants contracted the right FDI muscle at approximately 20% of the maximal voluntary contraction. EMG activity was displayed to participants in real-time using SIGNAL software (Cambridge Electronic Devices, Cambridge, UK) in order to ensure they maintained approximately 20% of maximal voluntary contraction of the right FDI immediately prior to and during the delivery of single-pulse TMS. Participants were instructed to relax the right FDI in between test pulses. The same biphasic stimulator was used to determine aMT and administer cTBS.

BDNF Genotyping
Genomic DNA from human saliva samples was collected in Oragene DNA collection kits and was then isolated using the prepIT.L2Preagent (cat # PT-L2P-5, DNA Genotek Inc., Canada) and precipitated with ethanol according to manufacturer's instructions. The DNA samples were genotyped for BDNF (the single nucleotide polymorphism rs6265) using the TaqMan SNP Genotyping Assay (C__11592758_10) designed by Thermo Fisher Scientific. Primers and probes were mixed with TaqMan Universal PCR Master Mix (Thermo Fisher Scientific). 4.5 µl of genomic DNA (2.5 ng/µl) was transferred in triplicate to a 384-well plate, with each well containing 5.5 µl of the PCR mixture. The PCR reaction was performed following a protocol provided by ABI. The allele was discriminated by post-PCR plate reading on the ViiA TM 7 System. Data were processed using the ViiA TM 7 Software (Thermo Fisher Scientific).
under investigation impact cTBS-induced changes in motor excitability above and beyond contributions from other factors that may affect MEP amplitudes and cTBS responsiveness. To account for age-related changes in neuroplasticity (Freitas et al., 2011), we first fit a base model with age as a covariate and sequentially added to the base model the fixed effect timepoint (baseline vs. 0, 10, 20, and 30 min post-cTBS). We then sequentially added the remaining fixed effects of interest (the categorical predictor representing BDNF status [Val66Val vs. Val66Met] and the continuous predictor representing the range of SI 1mV values) and tested whether the inclusion of each factor and/or interaction term significantly improved the model fit using the change in the deviance statistic (−2 times the log-likelihood), which follows a chi-square distribution with degrees of freedom equal to the number of parameters added. The main advantage of this analytic approach is that it improves the interpretability and generalizability of findings by accounting for the random effects that can influence model parameter estimates, resulting in more precise estimates of the fixed effects of interest (e.g., Warrington et al., 2014;Harrison et al., 2018). Moreover, it is robust to violations of normality (e.g., Zuur et al., 2010) and unequal group sizes (here, unequal n per BDNF genotype group) when compared to regression analyses conducted on mean values per subject (e.g., Field and Wright, 2011). Nonetheless, graphical tools for assessing linear mixed model assumption violations [e.g., normal quantile-quantile (QQ) plot of residuals] were evaluated (following the guidelines of Zuur et al., 2010; see also Harrison et al., 2018), revealing that the assumptions had not been violated.
All models included by-participant random intercepts to capture the inherent correlation among multiple measurements within a participant and adjust for individual differences present prior to the intervention. Following the recommendation of Barr et al. (2013), we attempted to implement a maximal random effects structure reflecting by-BDNF genotype status random slopes for the fixed effect of timepoint, as other fixed effects of interest were not fully crossed within participants (i.e., SI and BDNF genotype status) and the inclusion of by-participant random slopes of timepoint are not warranted theoretically given our prediction that the slope of the timepoint effect differs as a function of BDNF genotype status. However, this maximal random effects structure triggered convergence warnings (indicative of model overfitting), and did not significantly improve model fit (p = 1).
We conducted post hoc comparisons using the emtrends function in the Estimated Marginal Means R package (Lenth, 2020). We first assessed whether the change in MEPs from baseline to each post-cTBS timepoint differed as a function of the slope of SI within each BDNF genotype group separately, which allowed us to evaluate whether MEPs significantly change from baseline within each group and if that change differs as a function of increasing SI. We then evaluated between-group differences in MEPs across the slope of timepoint for the median of the upper and lower tercile of SI values (SI 65 and SI 44 , respectively) in order to assess whether the two groups differ in terms of how higher/lower SI affects change in MEPs from baseline to post-cTBS timepoints. The Tukey method was used to correct for multiple comparisons within each family of estimates.

RESULTS
Among the 33 participants, 20 were Val66Val carriers, 11 were Val66Met carriers and two were Met66Met carriers, consistent with the known prevalence of BDNF Met allele polymorphism (Shimizu et al., 2004). However, because there were only two homozygous Met allele carriers in the current sample and the extent to which Met66Met carriers behave similarly to heterozygous Met allele carriers remains to be clarified (e.g., Egan et al., 2003), data from these two participants were excluded from analyses (rather than collapsed into the Val66Met group; see e.g., Cheeran et al., 2008). Comparisons between the two groups revealed that there were no significant differences in age Thus, the two groups did not differ with respect to SI 1mV used to collect MEPs ( Table 1). All subjects tolerated the cTBS with no adverse effects. Note. Independent samples t-tests were performed unless data violated the assumption of normality and/or equal variances, in which the Mann-Whitney (U) tests were performed. Effect sizes for t-tests represent Cohen's d, whereas effect sizes for Mann-Whitney tests represent the rank biserial correlation. SI 1mV (defined as the percentage of maximum stimulator output; MSO) and % rMT reflects the SI used to probe plasticity prior to and immediately following cTBS. BDNF, brain-derived neurotropic factor; MEP, motor evoked potential; rMT, resting motor threshold; SI, stimulation intensity; M, mean; SD, standard deviation; mV, millivolt.  First, we assessed whether cTBS decreased motor excitability while accounting for potential age-related differences in cTBS response by adding to the base model the fixed effect of timepoint. Including the timepoint significantly improved model fit (χ 2 (1) = 27.97, p < 0.0001). We then evaluated the impact of BDNF genotype status on cTBS response by testing improvements in model fit (over that which included age and timepoint as single terms) with the addition of BDNF genotype status and the timepoint × BDNF interaction. Neither BDNF genotype status (p > 0.61) nor the interaction between BDNF genotype status and timepoint improved model fit (p > 0.45). We then explored whether test pulse SI (defined as MSO) impacted MEPs overall and/or cTBS-induced changes in cortical excitability as a function of BDNF genotype status. Adding SI did not improve model fit (p > 0.20), but including the interaction between BDNF, timepoint, and SI did significantly improve model fit (χ 2 (6) = 206.10, p < 0.0001). Model comparison results are reported in Table 2, and full model results including factor estimates and associated significance values are provided in Table 3.
Post hoc comparisons evaluating the effect of SI within each BDNF genotype group revealed that SI affected the degree to which MEP amplitudes differed from baseline to each post-cTBS timepoint for both Val66Val and Val66Met carriers. Specifically, the extent to which cTBS suppressed MEPs decreased as a function of increasing SI up to 30 min post-cTBS. Moreover, comparisons between BDNF genotype groups at the median of the upper and lower tercile of the range of SI values revealed no difference in the change in MEPs across timepoint for lower SI values (SI 44 ; p > 0.21), but a significant difference across timepoint for higher SI values (SI 65 ; p < 0.0001), with Val66Met carriers exhibiting (paradoxical) MEP-facilitation with test pulses   Figure 2 displays the full model predicted estimates of MEPs for each timepoint as a function of SI for each BDNF genotype group. Estimates were extracted from the model using the effects package in R (Fox and Weisberg, 2018), which yielded separate estimates per 10% increments of total percentage MSO across the range SI 1mV values for the two BDNF genotype groups and each timepoint. The model-estimated MEPs pre-and post-cTBS reflect predicted values after accounting for the covariate age and inter-individual variability in MEPs prior to cTBS (represented in the by-participant random intercept). As shown in Figure 2, the present findings demonstrate that test pulse SI has a differential impact on cTBS aftereffects among BDNF Val66Val carriers when compared to Val66Met carriers. While increasing SI attenuated cTBS-induced MEP suppression for both BDNF genotype groups, this effect was more pronounced for Val66Met carriers, as higher test pulse SIs yielded paradoxical (excitatory) cTBS responses from baseline to post-cTBS timepoints raw mean MEPs per SI interval for Val66Val vs. Val66Met carriers (corresponding to the predicted values depicted in Figure 2) can be found in Supplementary Table 2.

DISCUSSION
This study reports two key findings important for understanding potential sources of variability in cTBS response. First, SI impacted cTBS aftereffects, with both BDNF genotype groups exhibiting less MEP suppression from baseline to post-cTBS timepoints with higher SIs. Second, SI used to probe cortical excitability before and after cTBS accounted for variable cTBS aftereffects among Met allele carriers to a greater extent than their homozygous (Val66Val) counterparts. Thus, our data confirm and extend existing results demonstrating that inter-subject variability in response to cTBS arises due to intrinsic (genetic) factors (Cheeran et al., 2008;Mastroeni et al., 2013;Jannati et al., 2017) and extrinsic (methodologic) factors (Vallence et al., 2015;Goldsworthy et al., 2016b), and indicates that these two sources of variability interact in meaningful ways.
To date, only a few studies have systematically investigated whether BDNF Val66Met genotype status affects the response to cTBS (Cheeran et al., 2008;Mastroeni et al., 2013;Jannati et al., 2017). However, the findings from these studies have been mixed, precluding any definitive conclusions regarding whether BDNF genotype status contributes to inter-individual differences in cTBS response more generally. This suggests that cTBS response is more variable among Met allele carriers when compared to their homozygous (Val66Val) counterparts (Chung et al., 2016). Here, we demonstrated that the SI used to probe cortical excitability predicted both the magnitude and direction of cTBS-induced changes in plasticity for Met allele carriers to a greater extent than Val allele carriers, which provides a framework for understanding prior equivocal results and variability in cTBS response. Across the studies using different methods for determining SI to probe cTBS-induced plasticity, SIs evoking 0.5 mV yielded null effects of BDNF genotype status (Mastroeni et al., 2013), SI 1mV yielded little-to-no changes in cortical excitability for Met allele carriers (Cheeran et al., 2008), and SI at 120% rMT (evoking on average MEPs of 1.4 mV) yielded paradoxical effects of cTBS (MEP facilitation) for Met allele carriers (Jannati et al., 2017). When considered alongside the findings from the current study, this suggests that for Met allele carriers, cTBS becomes increasingly facilitatory with increasing test pulse intensity, especially when compared to homozygous Val allele carriers, which is consistent with the hypothesis that polymorphisms in the BDNF genotype constitute an important source of inter-individual variability in cTBS responsiveness.
Inconsistencies regarding the impact of BDNF polymorphisms on stimulation-induced neuroplasticity are not unique to cTBS protocols, as both significant and null effects of BDNF genotype status have been obtained in studies using several different noninvasive brain stimulation techniques, such as iTBS (Antal et al., 2010;Lee et al., 2013), paired associative stimulation (PAS; Cheeran et al., 2008;Witte et al., 2012), and cathodal tDCS (Cheeran et al., 2008;Antal et al., 2010). However, whether these discrepant findings also reflect methodological differences in the SI used to probe plasticity remains unclear. As in the current research, these studies assessed changes in cortical excitability using SI 1mV -a method that necessarily introduces variability as it pertains to the degree to which SI exceeds a given individual's resting motor threshold in order to attain MEPs of ∼1 mV. However, to our knowledge, no prior study has explored whether variability in the intensity required to achieve MEPs of 1 mV at baseline, i.e., prior to repetitive stimulation, interacts with BDNF genotype status. If true, our results predict that variability in SI among each study's sample would explain why Met allele carriers as a group either exhibit the expected pattern of response or reduced and/or paradoxical responses to plasticity-inducing noninvasive brain stimulation techniques. Future research is needed to fully explore this hypothesis.
It is important to note, however, that our finding of decreased MEP-suppression (and increasing MEP-facilitation among Met allele carriers) with increasing SI is not consistent with prior research exploring cTBS aftereffects across individuals' inputoutput (IO) curves. That is, evidence elsewhere demonstrates that cTBS-induced MEP-suppression is greatest with SIs ∼150% of rMT (Vallence et al., 2015;Goldsworthy et al., 2016b). However, these studies differ from the current research in two potentially important ways. First, cTBS intensity was determined using rMT, whereas the present research used aMT to determine cTBS intensity. Given that prior activation of the targeted muscle has been shown to affect both the magnitude and direction of cTBS aftereffects (Gentner et al., 2007;Iezzi et al., 2008), it is possible that this methodological difference also contributes to variability in cTBS response. Second, these studies did not investigate potential differences in SI as a function of BDNF genotype status, leaving it an open question as to whether or not pre-active vs. pre-relax muscle contraction would represent an additional source of variability depending on this genetic factor. Thus, it will be important to explore other potential interactions between genetic and experimentally-imposed sources of cTBS response variability in future work.
Why is there facilitation of MEPs, rather than the expected MEP-suppression, following cTBS among Val66Met carriers at higher SIs? While the precise mechanisms underlying this effect are beyond the scope of the current study, findings from prior work provide future avenues of study. It has been shown that the increase in the percentage of MSO required to produce the full range of possible MEP amplitudes can vary considerably across individuals and that both intrinsic and external factors contribute to this variability (Goetz et al., 2014). Evidence elsewhere suggests that test pulse SI along the IO curve impacts response to plasticity-inducing protocols (Vallence et al., 2015;Goldsworthy et al., 2016b). Moreover, some evidence suggests that, depending on whether the initial response to cTBS is in the predicted direction, increasing or decreasing cTBS intensity can reverse the direction of the initial response (Sasaki et al., 2018). In other words, the aftereffects of plasticity-inducing protocols like cTBS may vary as a function of both the intensity at which stimulation is administered and the intensity at which the cortex is probed before and following modulation. Our findings add further to this complexity by suggesting that the IO curve of Met allele carriers may inherently differ from that of homozygous Val allele carriers, which in turn affects the direction of their response to plasticity-inducing protocols like cTBS. Partial support for this hypothesis comes from studies that have explored whether rTMS intensity interacts with BDNF genotype status (Hwang et al., 2015; also see Jannati et al., 2017 for evidence indicative of a potential interaction between cTBS intensity and BDNF genotype status). We suggest that future studies investigate the ways in which stimulation protocols may influence response to plasticity-inducing protocols, and more importantly the potential interaction between methodological choices and BDNF genotype status. This may yield additional promising insights into the sources of variable rTMS response and in turn, the neurophysiology that dictates said variability.
There are a few limitations to our study. First, the sample size is relatively small, which may affect the generalizability of our results. However, several factors may militate against this concern. Our sample size is: (1) comparable to prior studies investigating the effects of BDNF polymorphisms on stimulation-induced changes in neuroplasticity-whether reporting positive (Cheeran et al., 2008) or negative results (Mastroeni et al., 2013); and (2) representative of the occurrence of BDNF polymorphisms in the population (Shimizu et al., 2004). Critically, the analysis approach adopted here is relatively robust to biases in estimated effects for between-subject factors with unequal sample sizes (Maas and Hox, 2005;Bell et al., 2014). This is because linear mixed effects modeling can use trial-level data, which captures variability between subjects prior to the intervention. We argue that this represents a major advantage over prior statistical methods that have obscured inter-and intra-individual differences by analyzing changes in an individual's mean MEPs across a block of trials, which may be contaminated by cumulative effects (Pellicciari et al., 2016) and/or transforming MEP values to reflect a percentage change in mean/maximum MEP values obtained at baseline (e.g., Sasaki et al., 2018). Second, SI used to probe plasticity varied across individuals, rather than within individuals. Thus, it will be important for future work to confirm and extend these findings by comparing the IO curves of individuals with and without the BDNF polymorphism as it relates to neuroplastic responses to rTMS protocols. Finally, it could be argued that the findings reported here were not specific to cTBS, as the current study did not include a sham-control arm-as is the case in the majority of studies investigating stimulation-induced motor excitability (e.g., Cheeran et al., 2008;Vallence et al., 2015;Goldsworthy et al., 2016b;Jannati et al., 2017Jannati et al., , 2019. However, we believe the lack of a sham-control condition does not undermine the main findings. Prior work that has included a sham-control arm, but did not investigate genetic moderators of rTMS response, demonstrates that cortical excitability changes in the expected direction following real, but not sham, stimulation (e.g., Todd et al., 2009;Albuquerque et al., 2018), suggesting that the effects of rTMS on MEP amplitudes are not due to order effects (i.e., prior elicitation of MEPs via single-pulse TMS) and/or placebo effects. Regarding the BDNF polymorphism, studies have shown that MEPs decrease to a greater extent for homozygous Val carriers when compared to Val66Met carriers in the absence of repetitive NIBS protocols (e.g., motor learning tasks; e.g., Kleim et al., 2006); however, the current study did not include a voluntary motor (learning) component. Thus, we would expect similar changes in MEPs from preto post-cTBS for both BDNF genotype groups if the findings were not specific to the stimulation itself. Nonetheless, it is still important to replicate and extend these findings in a sham-controlled study with a larger sample size sampled across the IO curve.

Conclusions
Overall, our data provide novel insight into the sources of variability in response to rTMS protocols, which has important implications for optimizing the clinical utility of this neurorehabilitative tool. The finding that SI used to probe plasticity modulated cTBS response for Met allele carriers to a greater extent than their homozygous (Val66Val) counterparts indicates that genetic factors interact with methodological sources of variability. Given that these two sources of variability interact in a meaningful way, it is important that future work take this into consideration. Further refinement of our understanding of the complex interplay between stimulation parameters that determine the effects of TMS and biologically-based individual factors that influence neuroplasticity may allow for further optimization and personalization of both experimental and therapeutic brain stimulation protocols.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Institutional Review Board of the University of Pennsylvania. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
RH, PS-B, and OF contributed to the design of the study. DH, LD, PS-B, RW, DS, and OF collected the data. DH, LD, PS-B, RW, DS, AA, AT, and OF assisted with pre-processing and analysis of neurophysiological data. FL performed genotyping of saliva samples. DH and LD analyzed and interpreted the full dataset and drafted the manuscript. All authors revised the manuscript, approved the final version, and agreed to be accountable for the content of the work. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by the National Institute on Deafness and Other Communication Disorders, National Institutes of Health (grant number R01 DC012780) and the Dana Foundation Brain and Immunoimaging Award, New York, NY awarded to RH. DH was supported by the National Institutes of Health Eunice Kennedy Shriver National Institute of Child Health & Human Development (grant number T32 HD071844). The study sponsors had no involvement in the collection, analysis and interpretation of data and in the writing of the manuscript.