A Pilot Study on the Functional Stability of Phonation in EEG Bands After Repetitive Transcranial Magnetic Stimulation in Parkinson’s Disease

1. Introduction

Parkinson’s disease (PD) is a neurodegenerative disease first described by James Parkinson in 1817.¹ Recent studies quantify its incidence in 15 cases per 100,000 people, with a prevalence ranging from 100 to 200 cases per 100,000.² The impact of PD on the quality of life of persons with PD (PwP) is associated with motor symptoms such as in effective neuromotor control, difficulty in walking and handling objects, resting tremor, facial rigidity, etc.,³ and other nonmotor symptoms which are also affecting PwPs’ ability to carry an independent life. PwPs also experience alterations in respiration, phonation, articulation, and prosody, collectively referred to as hypokinetic dysarthria (HD), which is a complex motor speech disorder characterized by manifestations such as mono-pitch and monoloudness, imprecise articulation, impaired speech rate and rhythm, and irregular pitch fluctuations.⁴ PD symptoms may be mitigated by medication, neurostimulation, and rehabilitation, and stabilized temporarily to facilitate PwPs’ motor functions.⁵ Transcranial magnetic stimulation (TMS) is a noninvasive method for the stimulation of neural tissue, including the cerebral cortex, spinal roots, and cranial and peripheral nerves. In this sense, repetitive TMS (rTMS) has been proposed as a therapy to improve patients’ neuromotor conditions.⁶

The main objective of this study is to explore the potential of assessing the functional improvement of phonation conditions of PwPs after rTMS treatment, using electroencephalography (EEG)-related frequency bands of glottal neuromechanical correlates. Phonation symptoms such as thyroarytenoid muscle hypotonia, vocal fold imbalance, and tremor in voice (altered neuromotor feedback) are some manifestations of PD-related neurodegeneration on speech.^7,8 The process of phonation is driven by lung air pressure forcing flow through the vocal folds. These structures are composed of muscle (thyroarytenoid) and connective tissues, fixed to the thyroarytenoid cartilage structure.⁹ The thyroarytenoid muscle, also known as musculusvocalis together with the cricothyroid, and transverse and oblique arytenoid muscles, is responsible for vocal production. Especially, the cricothyroid, and thyroarytenoid muscles are involved in contracting and stretching the vocal folds, thereby controlling phonation pitch. The inferior laryngeal nerve (recurrent) is the principal pathway in neuromotor innervation of the mentioned muscles in the larynx, except for the cricothyroid, which is innervated by the superior laryngeal nerve. The function of these nerves is to induce contraction and adduction of vocal folds to induce vibration under lung pressure.¹⁰ Therefore, it can be hypothesized that vocal fold tension (a biomechanical correlate of phonation) might be the direct consequence of the combined NMA of the inferior and superior laryngeal nerves (neuromechanical correlates), both being branches of the $X$ cranial nerve (vagus). The control of larynx neuromechanics can be directly related to the activity of premotor (PM), motor (M), and supplementary motor (SM) areas of the brain, which is ultimately related to the nature of speech being produced (voiced or voiceless). Therefore, meaningful cognitive traits determining prosody (energy, intonation, rhythm) will intervene in the control of vocal fold adduction, abduction, and tension. There are many brain areas implied in the laryngeal motor control, both in phonation and nonphonation tasks. Primary vocalM areas include the laryngeal motor cortex,¹¹ PM cortex, SM area, and cerebellar lobule VI. Secondary areas comprehend the cingulate M area, the ventral nuclei of the thalamus, the putamen, the frontal operculum, and the anterior insula.¹² The structural and functional relationship among them is still a matter of ongoing research, out of the reach of the present contribution, although there is a consensus in that as far as vocalization is concerned, the direct pathway is mainly controlled by the ventromedial central sulcus peak (VmCSP) corresponding to Broadmann area 4$p$, and the dorsolateral peak (DsP) in area 6, both on the precentral gyrus next to articulator control areas. In close relationship with this paper’s primary objective, Rödel et al.¹³ reported selective stimulation of the vocal fold tensor (cricothyroid) or relaxer (thyroarytenoid) using TMS. In this sense, Brown et al.¹¹ suggested that both areas mentioned before may have functional roles representing different muscles in the larynx because both muscular systems are innervated differently, specifically by the external laryngeal (cricothyroid) and the recurrent laryngeal nerves (thyroarytenoid).⁴ The working assumption in this study will consider that the main neuromotor activity (NMA) projected to the lower neuromotor units through the direct activation pathway of laryngeal neuromotor control is mainly related to VmCSP and DsP activity, given the fact that the exercises examined are solely the phonation (sustained vocalization) of a maintained vowel [a:], therefore, the role of secondary areas (indirect activation pathway) more relevant in other dynamic speech tasks would play a subsidiaryrole. The vocal fold stiffness (VFS) resulting from the tension and stretching of larynx muscles will condition the patterns of the glottal source (GS), a correlate of the pressure in the supraglottal ridge of the vocal folds, and the GS will induce the vocal sound pressure (VSP) at the radiation extreme when propagated through the oro-nasopharyngeal tract (ONPT). A simplified view of the information flow would be $\mathrm{\text{NMA}}\to \mathrm{\text{VFS}}\to \mathrm{\text{GS}}\to \mathrm{\text{VSP}}$. Specifically, the relationship $\mathrm{\text{NMA}}\to \mathrm{\text{VFS}}$ is known as cortical muscle coupling (CMC).¹⁴ As it is well known, alterations in NMA coherence seen in the EEG may be present in motor-related disorder treatment and rehabilitation,^15,16,17 a consequence of the working hypothesis is that these alterations might be in part modeled from VSP as frequency bands related to classical EEG activity using signal model inversion to reconstruct the way-back path as $\mathrm{\text{VSP}}\to \mathrm{\text{GS}}\to \mathrm{\text{VFS}}\to \mathrm{\text{NMA}}$. This possibility would benefit from ongoing studies on CMC, as it can be used to explore the coupling relationship of different functional frequency bands.^18,19 Therefore, this study proposes to conduct the reconstruction of the way back path from VSP to EEG-related frequency bands using long-lasting utterances of an open vowel as [a:] to produce functional phonation descriptions in terms of NMA and to use them in detecting functional phonation improvements after rTMS treatment. The use of maintained phonations of [a:] in speech pathology studies is well suited for functional phonation evaluation, this fact being recognized by its wide application in clinical practice, given that slight variations in NMA will be immediately reflected in GS, and thus, easily tracked and monitored.

Having fixed the conducting narrative justifying the reconstruction process from VSP to NMA, it should be decided which frequency bands would be of higher interest for being explored further. It would seem reasonable to focus on the activity in the $\vartheta$- and $\gamma$-bands following the description of the nonlinear character of motor unit recruitment in muscular agonist–antagonist activation in Darbin and Montgomery,²⁰ taking into account the direct neuromotor pathways involving the cortex-thalamus-basal ganglia (Ctx-Th-BG) circuitry activating the cricothyroid and thyroarytenoid muscles, and the projection of the organized oscillators that explain the ultimate nonlinear character of motor unit activity involved in the laryngeal function, because the $\vartheta$-band is strongly related to unstable nonlinear NMA,²¹ whereas, on the other hand, the $\vartheta$–$\gamma$ coupling should be related to the Ctx-Th-BG activity.²²

Thus, a primary objective of the study would be to evaluate the functional competence of phonation in TMS participants. A secondary objective is based on the existence of strong relationships between muscular contraction under biomechanical drive and neuromotor EEG activity on the brain areas responsible for PM and motor control^18,23,24 in the sense that laryngeal motor activity is transferred by larynx muscles, inducing the contraction of the thyroarytenoid muscle, measured by the unbiased vocal fold stiffness (UVFS) through nonlinear projection operators transforming neural discharges into muscle contraction. This modeling would allow reverse system inversion, provided that adequate inverse operators could be designed based on system identification methodologies. Therefore, it would be possible to advance in the projection of the NMA estimated from phonation biomechanics over the brain area activity, modeled as EEG frequency bands.

According to the primary objective, this work is intended to assess the validity of EEG-aligned frequency contents of glottal neuromechanics, but using only indirect estimations from vowel phonation acoustic recordings, to characterize the stability of phonation in pre-stimulus and post-stimulus vocal emissions from a limited set of PD participants, using a methodology aligned with the secondary objective. This paper is structured as follows. Section 2 gives a description of the participants in the experimental setup, the stimulation protocol, and the signal processing methods to estimate the signal correlates used in the study, with special attention to EEG-related frequency bands of the vocal fold biomechanical stiffness, and the statistical methods to evaluate phonation stability conditions from each participant, evaluation session, and frequency band. Section 3 presents a descriptive example of the results of processing a segment from a pre-stimulus phonation, the statistical distributions of frequency band estimates, the results from comparing pre- and post-stimuli distributions using log-likelihood ratios (LLRs) and hypothesis tests, and the global scores helping to assess the effects of stimulation on phonation stability conditions of participants. Section 4 discusses the meaning of results in terms of treatment efficiency and methodological accuracy, highlighting achievements and limitations. Section 5 summarizes the findings, contributions, and conclusions derived from the study. An Appendix has been added to include supplementary data tables.

2. Materials and Methods

2.1. Sample description

This study has been conducted on PwP participants showing mild to moderate HD directly related to PD, all of them right-handed, and native speakers of Czech. All were on stable dopaminergic medication for the duration of the whole study. The patients were tested in the ON state, which means that they had received the adequate dosage of dopaminergic medication according to their respective prescriptions 2$$h before the evaluations were conducted (patients showing levodopa-induced dyskinesias were excluded from the study). They were informed of the nature of the research and gave their written consent. The study considered only those cases presenting a pre-stimulus and four post-stimulus evaluations spaced in time. Eighteen cases were selected from the 33 participants included in the original study database⁵ because those were the ones fulfilling the before-mentioned condition. All other cases completed less than four post-stimulus examinations, and therefore, were not considered for this first study. The demographic description of the participants is given in Table 1. The cohort distributions are broadly similar in terms of UPDRS grade (females: $16.6\pm 4.1$; males: $12.3\pm 3.9$) and age (females: $74.6\pm 3.0$; males: $69.7\pm 8.4$).

2.2. Pre- and post-tTMS audio recordings

Each participant went through a baseline assessment (pre-stimulus evaluation at T0) before being submitted to 10 stimulation sessions (stimulation process) within the interval of two weeks; a follow-up evaluation session two weeks after stimulation (post-stimulus at T1); additional follow-up evaluations around six weeks (post-stimulus at T2), and around 10 weeks (post-stimulus at T3). The 18 participants of the subset included in this study were submitted also to a fourth post-treatment evaluation session around fourteen weeks after the stimulation process (post-stimulus at T4). The evaluation dates are listed in Appendix A in Table A.1. Each participant in the study was randomly assigned to active or sham stimulation. A perceptual assessment was conducted by a speech therapist rating speech performance (evaluation of dysarthria profile), faciokinesis, phonorespiration, and phonetic competence at each evaluation step. Audio recordings of the following utterances from each participant were taken before and after the stimulation process (pre- and post-stimuls): one free-topic monologue; one short neutral reading; one short emission of vowels [a:], [i:], [u:] (of $\sim 1.5$$$s); one long emission of a sustained [a:] (of around 15$$s); one long emission of a diadochokinetic exercise consisting in the repetition of [pataka] ($\mathrm{\text{lasting}}>10$$$s), and one single emission of 10 different selected tri-syllabic words. The trial was registered on clinicaltrials.gov (Number NCT04203615).

2.3. Speech processing methods

The study was conducted on 4$$s fragments of the long emissions of a sustained [a:], selected 2$$s after the vocal onset, to avoid phonation start transients and also fatigue effects on phonation.²⁵ The stimulation protocol and speech recording conditions are described in detail in Brabenec et al.⁵ Participants were submitted to rTMS (DuoMAG${}^{\text{TM}}$XT-100, Deymed Diagnostic) in 10 stimulation sessions over two weeks at CEITEC, Masaryk University. Each stimulation session took 40$$min. to complete, during which an eight-shaped coil over the right posterior superior temporal gyrus (STG, MNI coordinates $X=40$, $Y=-38$, $Z=14$) applied trains of pulses at a frequency of 1$$Hz, 100% intensity of the pre-estimated resting motor threshold (1800 pulses per stimulation session). In the case of sham stimulation sessions, the same coil was used in quite a similar fixture, producing the same sounds as in the active case described, but no magnetic field was applied. The stimulation thresholds and settings were established in a preliminary study, and because of their complexity, the interested reader will be referred to it for further clarification.⁵ Audio recordings were taken using a large capsule cardioid microphone M-AUDIO Nova mounted on a boom arm RODE PSA1 at a distance of approximately 20$$cm from the patient’s mouth. Acoustic signals were digitized by the M-AUDIO Fast Track Pro audio interface with a 48$$kHz sampling frequency and 16-bit resolution. In what follows the processing methods for the estimation of NMA descriptions consisting of EEG-aligned frequency bands of the biomechanical VFS estimated from speech recordings will be commented on

•	Fragments of 4$$s long from the recordings of the vowel [a:], were down-sampled at a sampling rate of 16$$kHz, selected between the time instants at 2$$s and 6$$s, to skip potential onset and decay effects. This sampling rate preserves the frequency contents of glottal signals.
•	The ONPT transfer function was evaluated by a 24-pole inverse adaptive lattice-ladder filter26 based on the iterative adaptive inverse filtering (IAIF) algorithm.27 The size of the filter has been fixed according to Akaike’s criterion,28 using a factor of 1.5 over the sampling frequency divided by 1000 for an optimal filter size on linear predictive error estimation, therefore, for a sampling frequency of 16$$kHz, the size of the filter has been set to order 24. The adaptive lattice-ladder inverse filter estimates a prediction-error polynomial reducing the speech segment being analyzed $s(n)$ to a residual $r(n)$ by classical deconvolution as $r(n)=h_T(n)\ast s(n)$, where $h_T(n)$ is the impulse response of the prediction-error polynomial emulating the inverse transfer function of the ONPT, such that $H_T(\omega )={\sum }_m{h}_{Tm}{e}^{-jm\omega }$, where $j=\sqrt{-1}$, $m$ is the filter coefficient order, and $\omega$ is the angular frequency, therefore, the lattice structure reconstructs the tube chain structure of the ONPT, and its associated transfer function is removed from the spectral contents of the speech signal. A description of the inversion filter details can be found in Gómez et al.29
•	The GS was estimated in pitch-synchronous cycles30 by numerically integrating the inverse filter residual.
•	The VFS ($\xi$), given in $\mathrm{\text{N}}\cdot {\mathrm{\text{m}}}^{-1}$ was estimated from the spectral tilt of the GS adjusted on a 2-mass model of the vocal fold biomechanics.29,31
•	The VFS was de-biased and de-trended by a moving-average filter.
•	The working hypothesis established that VFS is the direct consequence of the neuromotor activation of the cricothyroid and thyroarytenoid muscles. To sustain a given stable phonation frequency $f0$, a delicate equilibrium between both activations is necessary.11 This equilibrium is represented by an average baseline value of VFS (trend). Oscillations around this trend would reproduce small-signal alterations of the neuromotor activation, therefore, a de-trended VFS would produce a good estimate of neuromotor instability due to agonist–antagonist balance misadjustment. To obtain a frequency-band description of neuromotor instability, the de-trended VFS was processed by a bank of fifth-order Butterworth band-pass filters tuned at the respective EEG-related frequency bands, ($\delta$: $f\le 4$$$Hz; $\vartheta$: $4\mathrm{\text{Hz}}<f\le 8$$$Hz; $\alpha$: $8\mathrm{\text{Hz}}\le f\le 16$$$Hz; $\beta$: $16\mathrm{\text{Hz}}<f\le 32$$$Hz; $\gamma$: $f>32$$$Hz; $\mu$: $8\mathrm{\text{Hz}}<f\le 12$$$Hz). As a result, a set of de-trended VFS frequency-band time signals ${\xi }_{ij}^k(n)$ is produced, where $i=(0,\dots ,I)$ is the evaluation session index ($I=4$), $j=(1,\dots ,J)$ is the participant index ($J=18$), and $k=(1,\dots ,K)$ is the frequency band index ($K=6$), pointing to the six frequency bands defined above, and $n$ is the time index.

2.4. Data analysis

An estimation of the VFS was produced using the methods described in subsection 2.3, as VFS is assumed to be directly related to the NMA of the areas responsible for laryngeal control during phonation. Under the hypothesis stated in Sec. 1, this signal has been decomposed in frequency bands ${\xi }_{ij}^k(n)$ corresponding to EEG ($\delta ,\vartheta ,\alpha ,\beta ,\gamma$) activity. The amplitude distributions of each EEG-related frequency were estimated from their histograms. Distributions from post-stimulus recordings were compared with the corresponding ones from pre-stimulus conditions to produce explainable interpretations of potential behavioral changes in the phonation function.

The comparison methods proposed were based on LLRs and hypothesis tests. The LLR between two given probability density functions (pdfs) $p_i(\xi )$ and $p_j(\xi )$ may be defined as

where $\mathrm{\Omega }$ is the estimation interval. The LLR besides estimating the similarity between two distributions gives also a sense of comparison. Assuming the distributions under study $p_i(\xi )$ corresponding to the post-stimulus phonation to be compared to the pre-stimulus phonation $p_0(\xi )$, two possible situations may be considered, either $p_i(\xi )$ is more dispersed than $p_0(\xi )$, i.e. $\mathrm{\text{Var}}\{p_i(\xi )\}>\mathrm{\text{Var}}\{p_o(\xi )\}$, or $\mathrm{\text{Var}}\{p_i(\xi )\}<\mathrm{\text{Var}}\{p_o(\xi )\}$, where $\mathrm{\text{Var}}\{\cdot \}$ is the variance of the distribution. In the first case, let us assume $\xi \in {\mathrm{\Omega }}_b$ to be the interval where $p_o(\xi )>p_i(\xi )$, and $\xi \in {\mathrm{\Omega }}_a$ to be the interval where $p_i(\xi )>p_o(\xi )$, where ${\mathrm{\Omega }}_a\subset \mathrm{\Omega }$ and ${\mathrm{\Omega }}_b\subset \mathrm{\Omega }$ provided that ${\mathrm{\Omega }}_a\cup {\mathrm{\Omega }}_b=\mathrm{\Omega }$, and ${\mathrm{\Omega }}_a\cap {\mathrm{\Omega }}_b=Ø$. Therefore, from Eq. (1) it will follow that

$$\lambda (p_i|p_o)={\int }_{\xi ϵ{\mathrm{\Omega }}_a}log\{p_i(\xi )∕p_o(\xi )\}d\xi +{\int }_{\xi ϵ{\mathrm{\Omega }}_b}log\{p_i(\xi )∕p_o(\xi )\}d\xi ;$$(2)

where

It may be seen that the sign of one integral is the opposite of the other, therefore, it will be expected that in the case of unimodal distributions when there is a strong difference in variance $\mathrm{\text{Var}}\{p_i(\xi )\}\gg \mathrm{\text{Var}}\{p_o(\xi )\}$, $p_o(\xi )$ will be much narrower than $p_i(\xi )$, and consequently ${\mathrm{\Omega }}_b$ will be much smaller in size than ${\mathrm{\Omega }}_a$, and $\lambda (p_i|p_o)=I_a+I_b<0$. Conversely, when $\mathrm{\text{Var}}\{p_i(\xi )\}\ll \mathrm{\text{Var}}\{p_o(\xi )\}$ the opposite condition will prevail ($\lambda (p_i|p_o)>0$). This is especially evident in normal and near-normal distributions. To put it otherwise, given the properties of probability densities, if $p_i(\xi )>p_0(\xi )$ on the interval ${\xi \in \mathrm{\Omega }}_a$, most likely it will be expected $p_i(\xi )$ to be narrower (lower variance) than $p_0(\xi )$, or in other words, that the generating process of $p_i(\xi )$ would produce less dispersed outcomes (more stable) than that of $p_0(\xi )$. Therefore, it is expected that functional improvements in phonation will produce lower variance post-stimulus frequency band distributions, and therefore, positive LLRs, and on the contrary, worsening phonation conditions will produce negative LLRs.

The significance of the comparisons among pre- and post-stimulus distributions was assessed by Mann–Whitney (MW) $U$-tests, on the null hypothesis of equal medians, because typically $p_i(\xi )$ distribution patterns might differ from normality. The tests were carried out on EEG-related band-frequency feature samples ${\xi }_{ij}^k(n)$, therefore, hypothesis tests for a given participant $j$ and a given feature$k$ would be conducted on each post-stimulus sample ($i=1,\dots ,I$) concerning the pre-stimulus one ($i=0$), would be given by

where $T_{\mathrm{\text{MW}}}\{{\xi }_{ij}^k(n),{\xi }_{0j}^k(n)\}$ is the MW $U$-test between pre-stimulus and post-stimulus samples ${\xi }_{0j}^k(n)$, ${\xi }_{ij}^k(n)$, $h_{jk}$ being the $p$-value estimated by the test, assuming that ergodic conditions apply.³²

The above-described methods allow us to determine the behavior of each frequency band based on LLR estimations (1) to explore whether these features improve significantly as a result of the intervention by rTMS, with special attention to which frequency bands would be more sensitive to changes in functional behavior. A global score to describe the general behavior of potential improvements regarding a specific participant may be defined on the averages over time of each feature sample $〈{\xi }_{ij}^k(n)〉$. As ${\xi }_{ij}(n)$ is an unbiased estimation of VFS, it may be associated with the tremor (oscillating instability) of the vocal fold. Its amplitude is expected to be larger the more acute the functional disorder affection of each participant’s phonation. Therefore, lower values of $〈{\xi }_{ij}^k(n)〉$ will be associated with a less unstable phonation, and therefore, with larger stimulation beneficial effects. Having this consideration in mind, the following definition was used as a normalized weighted score associated with each frequency band feature per evaluation session :

where the weight $w_i=(d_i-d_0)∕(d_I-d_0)$ is a normalizing factor to take into account the time interval between each post-stimulus date $d_i$ and the corresponding pre-stimulus date $d_0$ normalized to the longest interval $(d_I-d_0)$ in days. In this way, long-lasting beneficial effects were given larger importance than short-duration effects. As it may be easily checked, improving phonation features produce negative scores $c_{ij}^k<0.$

Another relevant score was defined on the progression trend of improvements, based on the first difference of average estimations between successive recordings, as

therefore, progressive improvements in a given feature would produce also negative scores $d_{ij}^k<0.$ Both scores, for relative and progressive improvements, were fused in a single one as $g_{ij}^k=c_{ij}^k+d_{ij}^k$ The total score per participant was defined as the average of all post-stimulus evaluation sessions over all frequency bands : where $I-1$ is the number of post-stimulus evaluation sessions (four in this study), and $K$ is the number of frequency bands considered (from $\delta$ to $\gamma$).

3. Results

An example from a sustained emission of the vowel [a:] during 4 s by one of the participants actively stimulated with rTMS is shown in Fig. 1. This specific example is included as a prototype to describe the speech processing protocol, not to be taken as a mark of generalized behavior, but as a particular phenomenological description of the examination procedures conducted on each phonation analyzed, worth of being examined in detail as a singular case of study. The VFS was estimated from the GS using the methods described in subsection 2.3.

The $\vartheta$-band and $\gamma$-band frequency distributions, corresponding to the best-behaving active case (1400) are presented in Figs. 2(a) and 2(b). A similar set of $\vartheta$-band and $\gamma$-band distributions, corresponding to the worst-behaving sham case (2200) is shown in Figs. 2(c) and 2(d).

The results of evaluating the LLRs between the pre-stimulus (T0) and the four post-stimulus on the frequency bands following expression (1) are given in Table 2.

The results of the corresponding MW tests following expression (4), availing the relevance of the LLRs are given in Table 3, andthe global scores (7) are given in Table 4.

A very relevant question has to see with the effects of rTMS on perceptual assessment of dysarthric behavior performance. This meta-information might serve to complement other studies, therefore, it was evaluated by a speech therapist and a neurologist for each participant at each evaluation session. For such, a three-function battery of tests (3F) was used to evaluate faciokinesis, phonorespiration, and phonetics (articulation), each one of them scoring between 0 (lowest score, poorest performance) to 30 (highest score, best performance). A total dysarthric score (DS) was estimated by summing the three mentioned scores, ranging from 0 (anarthria), 1–15 (profound dysarthria), 16–31 (severe dysarthria), 32–47 (moderate dysarthria), 48–63 (mild dysarthria), 64–79 (very mild dysarthria), to 80–90 (no dysarthria). The following table presents the evaluation of DS for each participant in Table 1, and the relative percentage of evolution after T4 versus T0 to be estimated as $I_{\mathrm{\text{DS}}}=({\mathrm{\text{DS}}}_f-{\mathrm{\text{DS}}}_0)∕{\mathrm{\text{DS}}}_0$, where $f=3$ for case 0100, and $f=4$ for all other cases.

4. Discussion

As already stated in the introduction, this work is intended to characterize the stability of phonation in pre-stimulus and post-stimulus vocal emissions following rTMS by assessing the validity of features estimated on glottal neuromechanics from a limited set of PD participants. The example depicted in Fig. 1(a) was given to illustrate the protocol followed in the assessment. This example is especially meaningful, as it shows two events of phonation activity breaks typical in PD dysphonia, at intervals 4.0–4.2$$s and 5.1–5.2$$s, resulting in a muscle stiffness drop at 4$$s appreciated in Fig. 1(b), which is not that evident at 5.5$$s, although the stiffness intensification from proprioceptive correction is quite well appreciated in both cases on the $\delta$-band (Fig. 1(c)). A relaxation towards stable muscle tone is observed between 4.2$$s and 5.0$$s as a follow-up after stiffness restoration. Intuitively, it could be hypothesized that the $\delta$-band signal would be associated with the group activity of sets of motor units reacting to proprioceptive feedback (coarse tuning), whereas the $\gamma$-band could be related to auditory feedback-induced motor actions to maintain muscle tone activity as stable as possible according to tonal settings (fine-neuromotor control of the thyroarytenoid muscle tension). Although the findings presented in Fig. 1 are rather specific to that case, similar $f$0 blocking events are not infrequent in PD dysarthria. These events and their associated EEG-band activity reveal interesting information, which could help in clarifying the phenomena behind muscular blocking and neuromotor failure.^33,34 The first important observation is that an apparent correction takes place immediately to neuromotor blocking, probably from the auditory system, as revealed in the $\delta$-band, which shows clear corrective actions to recover muscle tone. The $\vartheta$-band (Figs. 1(f) and 1(g)) shows intensified activity before and after the first blocking, and a moderate degree of tremor all over the whole interval studied. But interestingly, the $\alpha$-band shows a low level of activity before each blocking, which is reactivated immediately after (Figs. 1(h) and 1(i)). It is also seen as a ripple in the VFS (Fig. 1(b)). Some similar phenomenon is observed in the $\beta$-band (Figs. 1(j) and 1(k)) where a strong burst of activity is seen after the first blocking following an interval of very low activity. On the contrary, the activity on the $\gamma$-band (Figs. 1(l) and 1(m)) is very intense before the blocking, to become more diffuse and less intense after the incident. The activity in the $\mu$-band (Figs. 1(n) and 1(o)) seems to intensify after each blocking. These observations open the possibility of conducting studies on $f$0 blocking, depending on the dual neuromotor mechanism influencing VFS, on the one hand, the cricothyroid muscle innervated by the inferior laryngeal nerve, and on the other hand, the thyroarytenoid muscle innervated by the recurrent laryngeal nerve,¹² pointing to misadjustment in the agonist–antagonist balanceexerted by these muscles when failing to ensure fine $f$0 tuning. Of course, this assumption would demand further studies to provide generalization insights on the VFS neuromotor driving functionality.

The information provided by the activity in the $\vartheta$ and $\gamma$ bands as presented in Figs. 2(a) and 2(b) for the same case (active stimulation case 1400) reveals an important decay of phonation instability after the stimulation, where the tremor distributions reduce notably their variance and average values, the effect being more strongly observed in the $\vartheta$-band (Fig. 2(a), T1). This drastic improvement, in value and dispersion, is slightly worsened in the next observations, although post-stimulus tremor amplitudes keep well under 1%. This observation is especially relevant, as the $\vartheta$-band is classically associated with the so-called “low-frequency tremor” by many studies in voice quality, as it is easily perceived by listening, and it is typically associated with many cases of PD phonation.^35,36

As it was mentioned before, bands $\vartheta$ and $\gamma$ are considered especially relevant in the study according to previous knowledge on EEG-related NMA.³⁷ In this sense, it will be of special interest to review the behavior of the best and worst functional behavior cases comparing pre- and post-stimulus feature distributions in both bands, corresponding to cases 1400 and 2200 shown in Figs. 2(a)–2(d). The boxplots in Fig. 2 show the distributions of the tremor amplitude in the $\vartheta$ (a) and $\gamma$ (b) frequency bands from case 1400, corresponding to the same participant (pre-stimulus code T0, post-stimulus codes T1, T2, T3, and T4). Assuming that tremor is associated with phonation instability, a decrement in the tremor amplitude should be considered an improvement. Therefore, as boxplots, T1–T4 (post-stimulus) show a smaller median and a much smaller interquartile dispersion than boxplot T0 (pre-stimulus) in both bands ($\vartheta$ and $\gamma$), it might be concluded that the post-stimulus tremor is less unstable than pre-stimulus one, and consequently, that an improvement is observed. The situation concerning Figs. 2(c) and 2(d) is just the opposite, as it may be seen that post-stimulus distributions (T1–T4) show larger medians and much larger interquartile dispersion than the pre-stimulus one (T0) in both bands ($\vartheta$ and $\gamma$), therefore, it might be concluded that the situation, contrary to being improved, has deteriorated substantially.

These two cases showed the largest improvements and the largest deterioration in phonation stability. Whereas improvements seen in 1400 could be attributed to the beneficial effects of rTMS, the increment in instability shown in case 2200 could not be attributed to the effects of nonstimulation, therefore, a possible explanation would be a worsening in the phonation conditions within a short time interval between T0 and T4. Some other circumstances could have influenced the worse behavior observed. The same criticism could apply to case 1400, which experienced a strong improvement because recording T0 showed a large instability, which was not shown in T1–T4. Of course, the possible benefits of rTMS might not be the only explanation for such evolution or not its single cause. Obviously, this question is fully open to discussion.

Table 2 presents interesting results when analyzing the phonation instability of each participant before and after stimulation. It may be appreciated that a subset of the participants who underwent active stimulation (0800, 1100, 1200, 1400) shows mainly improvements in most of the evaluation sessions on the two bands of reference. Another subset of sham cases (1600 and 1900) experienced unexpected improvements, whereas two other sham cases (1700 and 1800) showed worsening behavior (quite strong in 1800, possibly pointing out to an extremely stable condition in T0). Another subset ofactive cases (2000, 2700, 2800, 2900) experienced slight or no improvements at all. A fourth subset of sham cases showed worsening behavior (2200 and 2400), whereas the remnant sham cases expressed mixed behavior (2300, 2500, and 2600).

As it is expected, a reduction in variance following a phonation stability improvement is observed when comparing each single pre-stimulus frequency band of the VFS estimate to the corresponding frequency band estimate from each post-stimulus recording in a longitudinal sequence (intra-participant) as given in Table 2 in terms of LLRs, validated by corresponding nonparametric MW $U$-tests shown in Table 3, whereit may be seen that most of the tests reject the null hypothesis, pointing to the fact that the majority of tremor bands examined showed significant differences between each first exploration in time (T0) and the subsequent ones (T1–T4). The tests failing to reject the null hypothesis correspond to small negative LLR values ($-0.1<\lambda <0$) in Table 2 (italicized) associated with potentially slight deterioration, but its significance is not sustained by the tests given in Table 3, except for case 1600 T1$\gamma$-band, producing a residual LLR of 0.012, its significance not being availed by the MW $U$-test ($p$-value of $0.172>0.050$). Besides, a transversal (inter-participant) $t$-test has also been conducted to assess how general improvement scores given in Table 4 compare in the context of all participants.

Table 4 presents the summarized overall results indicating which cases showed better phonation stability conditions confronting pre- and post-stimulus examinations considering all participants, evaluation sessions, and frequency bands. The global score $g_j$ per participant $j$, given in the fourth column from left, is estimated including all post-stimulus evaluation sessions (T1–T4) and five out of the six frequency bands ($\delta ,\vartheta ,\alpha ,\beta ,\gamma$), as per expression (7). The fifth column (Agreement) assigns a “1” to the cases where active stimulation produced a negative global score (True positives, meaning improvement of tremor instability, in agreement with active stimulation expectations), and those cases of sham (non) stimulation which presented a deterioration of tremor instability (true negatives, corresponding to positive global scores, meaning no improvement). The true positives and true negatives are listed in columns TP and TN. Correspondingly, those sham cases (nonstimulus), showing improvement (false positives), and active stimulus cases showing deterioration (false negatives), are listed in columns FP and FN, respectively. The total amounts of TP, TN, FP, and FN are given at the bottom of the respective columns (sixth to ninth). Fixing the detection threshold at $g_t=-0.1$ (assuming a slight improvement), the sensitivity, specificity, and overall accuracy of the stimulation methodology would be around 56%, 78%, and 67%, respectively, corresponding to an F1 score of 0.63. The agreement rate including TP and TN concerning the whole set size would be around 71%. Twelve cases presented results aligned with expectations: five active cases produced scores under the threshold $g_t=-0.1$ (improvement in phonation stability), and seven sham cases produced scores over the threshold (worsening in phonation stability). The other six cases produced results contrary to expectations: four active cases produced scores over the threshold (worsening in phonation stability), and two sham cases produced scores under the threshold (improvement in phonation stability). As a general comment, two out of three cases included in the study behaved accordingly to expectations after having been submitted to the stimulation protocol, including active and sham cases.

A possible explanation for the relatively large amount of FP and FN cases could be attributed to the small sample size, as well as to the effects of possible confounding factors affecting the phonation of participants during the tests, such as possible different types of HD involved, and the comorbid and emotional uncontrolled conditions affecting vocal emissions during the evaluation sessions (aging, depression and anxiety, medication, among them³⁸). These effects could play a rather important role in the evaluation results, considering the amount of time between the pre-stimulus and the post-stimulus examinations, extending between 93 (minimum) and 119 days (maximum, see Appendix A Table A.1). Long time intervals might make the evaluation protocol susceptible to the appearance of disruptive events, in the participants’ daily life conditions, such as episodes of emotional disturbance, changes in mood, dissatisfaction, discouragement, helplessness, depression, or despondency, as well as transitory respiratory or pharyngolaryngeal disorders, which could affect vocal production in one or another way, not being considered. Given the apparent high sensitivity of frequency bands to changing conditions in vocal fold unexplained alterations, minor changes in these confounding phonation conditions could produce substantial changes in amplitude distributions such as altering LLRs or null hypothesis rejection results. A further extended study should concentrate on comparing each examination from each participant to observe the extent to which improvement or worsening phonation stability conditions could be observed among them.

An important aspect to be discussed is the robustness of the approach, and if it can establish a causal relationship between rTMS and improvements in PD HD. Given the premises of the experimental setup, the comparison between pre-stimulus and post-stimulus results was undertaken using LLRs, because these ratios could be interpreted in terms of similarity and improvement sense. The fact that the evaluations conducted include a full longitudinal analysis with four post-stimulus evaluations taking into account the time interval having passed after stimulation, added important robustness to the study because improvement assessment was not based on a single post-stimulus evaluation, therefore, a tendency may be inferred from expressions (5) and (6) summarized in a single global improving score given by Eq. (7), because each single summarized score is supported by four longitudinal estimations comparing $〈{\xi }_{ij}^k(n)〉$ vs versus $〈{\xi }_{0j}^k(n)〉$, plus three others comparing $〈{\xi }_{ij}^k(n)〉$ vs versus $〈{\xi }_{i-1j}^k(n)〉$, having also the time intervals T_i–T₀ into account. The LLRs given in Table 2 have been estimated using pdfs generated from normalized histograms of 100 bins summarizing tremor amplitude time series on a cycle-synchronous basis, which convey around 400 samples per phonation segment in a male voice, and about double size in a female voice. These same distributions have been used in computing MW $U$-tests given in Table 3 to avail LLRs’ significance. Because the checking for false-discovery cases has been conducted already using MW $U$-tests (Table 3) validating the results given in Table 2 of not being produced by chance when their respective $p$-values are below 0.05 (confidence level), global scores have to be taken as reliable enough in the best and worst behaving cases, with the exceptions marked in bold in Table 3. Therefore, it could be asserted that rTMS seems to have some influence on phonation stability in the cases commented with a low margin of error. The reasons why these beneficial effects are not seen in other active cases can also be asserted with little margin for error. Why phonation stability may benefit from rTMS in some cases, and not in others, is a matter which cannot be determined within this study framework, therefore, the debate on causality is to be left for further study.

The general performance reliability has been also evaluated by a transversal (inter-participant) $t$-test, been conducted to assess how general improvement scores given in Table 4 may be compared in the context of all participants. This study follows up on our previous work published by Brabenec et al.⁵ In this previous work, we evaluated the effect of rTMS perceptually utilizing the 3F test dysarthric profile.³⁹ The test consists of three subtests: I. Faciokinesis, II. Phonorespiration, and III. Phonetics. The phonorespiration subtest could be further divided into the assessment of: (a) respiration, (b) respiration during phonation, and (c) phonation. We observed that the rTMS had a significant inter-speaker effect on articulation, prosody, and intelligibility (as assessed by the III. Phonetics subtest). However, we did not identify any significant impact on the II. Phonorespiration sub-score, which is in line with our findings, because in this study, the general performance reliability has been also evaluated by a transversal (inter-participant) $t$-test. It has been conducted to assess how general improvement scores given in Table 4 may be compared in the context of all participants. The hypothesis $t$-test between the global scores from the set of actively stimulated cases versus the set of sham cases in Table 4 will produce a $p$-value of 0.2, which does not reject the equal means between both sets. Does it mean that the functional assessment methodology proposed fails in its objectives? Evidently, not. The most plausible interpretation is that because actively stimulated cases 2700, 2800, and 2900 show a large deterioration after stimulation, the average global score of the active set (0.36) is not far enough from the average global score of the sham set (0.87) for the two sets overlapping variances (1.18 and 1.61, respectively). Nevertheless, this observation does not invalidate the individual results, because these are evaluated longitudinally on their timeline statistics, independently of how any other participants arebehaving, as this transversal comparison points out. It must be added that longitudinal tests would be in full alignment with the old medical lemma “Treat the patient, not the disease”, which is the grounding floor emphasizing the balance ofintra-participantover inter-participant studies.

Regarding the results from the perceptual score evaluation given in Table 5, it may be seen that some improvements are appreciated in all the cases, although only those over 5% could be relevant in alignment with active stimulation, corresponding to cases 0100, 1100, 1200, 1400, and 2700, and correspondingly, in sham cases 1600, 1700, 1800, 1900, 2200, 2300, 2400, 2500, and 2600, being irrelevant in the remnant cases. It might be interesting to compare the global scores after the phonation assessment (Table 4) to the dysarthria assessment results (Table 5) to verify if there is some relationship between phonation and perceptual improvements. Pearson’s correlation coefficient between both scores gives a result of –0.06, indicating that there is almost no correlation between phonation and articulation results, although a two-tail $t$-test on equal medians and different variances does not reject the null hypothesis under a 0.05 significance level ($p$-value equal to 0.11). This result should not be unexpected, because phonation and articulation are under the control of quite different neural pathways, and they need not be affected similarly in all the cases considered.⁴

In summary, the findings shown in Table 4 support the use of LLRs and conducting hypothesis tests focusing in particular on the $\vartheta$- and $\gamma$-bands, as these bands are believed to summarize well both neuromotor and cognitive activity.^21,40,41 In this sense, it must be stressed that besides its clear neuromotor character, phonation is also associated with a clear cognitive profile. The compelling argument is that speakers require mastering a real-time perception of their phonation performance to implement and conduct pitch and loudness fine adjustments, therefore, phonation has a marked undeniable cognitive character.⁴² Different brain areas, such as the primary motor phonation area (VmCSP and DsP), the cerebellum, the basal ganglia, the superior temporal and lateral prefrontal cortex, and the SM area, among others, are involved in the direct, indirect, and feedback neuromotor control circuitry,^4,43,44 in what seems to be a clear attention-driven cognitive function. The apparent synchronized activity around phonation motor blocking and disruption events could demand coordinated action of motor units spiking at different frequencies.^36,45,46,47 Of course, the clarification of all these observations would require further studies combining other cooperating methods, such as EEG, although at the cost of complicating signal acquisition and pre-filtering to remove facial muscle activity during speech production, possibly by the use of surface electromyography (sEMG).

Although the reach of this study is hampered by evident limitations, such as the nonexhaustive protocol and the relatively small sample size and gender unbalance, some insights on the efficacy of the stimulation methodology and the data analytics methodology used might still be drawn from the results presented, although they had to be taken more as speculative assumptions to open new research lines than proven facts based on exhaustive testing.

On the one hand, the effects of rTMS on phonation stability offer mixed results. Some of the active cases studied report notorious improvements in the frequency bands studied, which could not be explained by a by-chance effect, whereas some other actively stimulated cases do not seem to show phonation stability improvements, or even show more unstable phonation. Conversely, some sham cases show clear improvements, whereas some others show undeniable worsening phonation stability. Relating both behaviors to the absence of stimulation is not an easy task, opening an issue for further discussion. It seems that phonation stability behavior might be too much sensitive to confounding factors to serve as a unique marker by itself, and it should be combined with other speech-based traits, including articulation acoustic features as well. Besides, tremors might not be a bi-univocal feature of PD HD.³⁸ In summary, it could be asserted that rTMS seems to have some influence on phonation stability in the cases commented with a low margin of error, although these beneficial effects are not seen in other active cases, also with a little margin of error. Why stimulation may seem beneficial for phonation stability in some cases, and not in others, is a matter which cannot be determined within this study framework, and it is to be left for further study.

On the other hand, it seems that the study of phonation instability to monitor disruption events in VFS, and its association with EEG-related frequency bands would be fully justified as a powerful introspective methodology to disclose interesting NMA in cortical areas affecting larynx control. In this sense, this approach would be well aligned with relevant studies on the structural complexity of the brain from EEG frequency sub-bands,⁴⁸ the functional connectivity pattern assessment from EEG signals,^49,50 or the detection of movement intention in brain-computer interface (BCI) systems from EEG signals,⁵¹ among others.

The analysis being proposed uses instability distributions of the VFS estimated from speech utterances instead of EEG recordings to evaluate differences in functional behavior based purely on acoustic signal analysis, shedding light beyond what can be provided by classical acoustic analysis. The ultimate objective driving future extensions of this study is to reproduce EEG-aligned descriptions of phonation estimated from audio recordings only, which could be used in neurodegenerative speech characterization. The justification for characterizing the VFS biomechanical correlate by frequency bands aligned with those of EEG studies, as if it were an EEG channel signal, comes from preliminary work done on the NMA of the muscles involved in articulation,^52,53 using formants and glottal features. This relationship seems natural because the activity of the larynx, oropharyngeal, and facial muscles is the ultimate cause of source modulation and framing into acoustically perceived speech. Therefore, if muscles would respond to neuromotor stimulation, and eventually, to brain cortical activity, a way to study the intervening links of interest in modeling possible dysfunctions in the activation chain, would consider the important advancements in modeling and activity coding of CMC, to align EEG and acoustic signal descriptions on a common code.

This alignment offers the possibility of conducting a more ambitious, extensive, and exhaustive study including the combination of EEG recordings and other traits derived from articulatory movements to include EEG, sEMG, and audio recordings. In this sense, recent advancements in brain connectivity combining EEG, MEG, fMRI, and NIRS characterization by graph theory^54,55,56,57 and probabilistic neural networks⁵⁸ could offer new insights for future studies. The application of this methodology to synchronized mixed EEG-audio databases⁵⁹ should offer new insights into speech production comprehension and eventually might allow further disclosure of brain functionality and physiological responses in PD.

Other possible methodologies for decomposing and characterizing the biomechanical correlate VFS derived from the acoustic phonation signal would include Gabor transform,⁶⁰ Wavelet filter banks,⁶¹ Fourier-based synchrosqueezingtransform,⁵¹ transferentropy,⁶² or fuzzy synchronization likelihood,⁴⁸ among others. These studies would require a further extension of the present paper out of reasonable limits; therefore, they are considered for future research. Likewise, phonation improvements from rTMS based on comparisons using glottal source features, such as the harmonic-noise ratios, first-second harmonic ratio, cepstral peak prominence, parabolic spectrum coefficient, open and closed quotients, and other indices of voice quality analysis⁶³ could enrich the assessment protocol within a future study framework. Another important question left for future studies is the potential comparison of dysarthrophonia tests and limb bradykinesia based on a new assessment protocol.

5. Conclusions

In this study, the possibilities of predicting the interactions on the EEG-related $\vartheta$–$\gamma$ frequency bands of the NMA from the phonation acoustical signal have been explored. Although the size of the sample studied is a limit to the findings observed and these results are tentative given the limited number of participants and need to be verified in longer and larger trials, it may be asserted that the preliminary results are well aligned with ongoing studies in the field, especially in the use of LLRs to assess functional improvements in phonation after rTMS. This finding is based on the possibility of using the VFS to serve as a correlate to monitor disruption events in vocal emission attributable to PD consequences. Consequently, visualizing EEG-related frequency bands could help in understanding some of the phenomena underlying vocal emission disruptions. The positive effects of rTMS are evident in the results, although observations on phonation instability behavior from active and sham cases in pre-stimulus and post-stimulus recordings offered mixed results, concluding that the reasons behind findings need to be further explored and explained. There is a clear promise in these tentative findings grounded on previous work that can demonstrate differences in functional behavior based purely on acoustic signal analysis, which delves deeper and provides insights beyond what can be assessed by classical acoustic analysis and allow one to get closer to brain functionality and physiological responses in PD. However, we appreciate that there may be limitations and confounding factors when examined in detail such as emotional or comorbid conditions which might alter the sensitivity of functional assessment, and their effects should be taken into account in future studies.

Acknowledgments

This research received funding from European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement No. 734718(CoBeN), a grant from the Czech Ministry of Health, 16-30805A, a grant from EU– Next Generation EU (Project No. LX22NPO5107 (MEYS)), and grants TEC2016-77791-C4-4-R (Ministry of Economic Affairs and Competitiveness of Spain), and Teca-Park-MonParLoc FGCSIC-CENIE 0348-CIE-6-E (InterRegProgramme). Andrés Gómez-Rodellar holds a scholarship from the Medical Research Council Doctoral Training Programme in the Usher’s Institute (University of Edinburgh Medical School).