1 Introduction

Neurodegenerative diseases are progressive, currently incurable conditions affecting millions of people worldwide, with an increasing prevalence owing in part to the extensions in lifespan (Gitler et al., 2017). As multisystem disorders, neurogenerative diseases progressively affect multiple domains of function, including communication (speech and language), cognition (memory and executive function), swallowing, mobility, amongst others. The clinical manifestations and progression of neurodegenerative diseases are highly heterogenous across individuals, and the relation of such clinical heterogeneity to the underlying pathophysiology remains poorly understood. This knowledge gap posits a major barrier to the development of effective therapeutic treatments. Currently, the mainstay of the management of neurodegenerative diseases focuses on controlling symptoms, slowing functional declines, preserving quality of life, and prolonging survival, through both pharmacological and nonpharmacological interventions (Church, 2021; Norris et al., 2020). To optimize the outcomes of these interventions, a patient-centered approach has been advocated, emphasizing delivering the right intervention to the right person at the right time (Miller et al., 2018; Olszewska and Lang, 2023). This approach, also known as precision or personalized medicine in a broader context (Pokorska-Bocci et al., 2014), requires a better understanding of the pathophysiological processes and phenotypes that underlie the clinical heterogeneity to improve intervention tailoring based on individual characteristics.

To better understand the pathophysiological mechanisms of neurodegenerative diseases, a network-based perspective has been posited to link three aspects of neurodegeneration: (1) the variable spreading patterns of pathogenic processes, (2) the degeneration of distributed neural networks, and (3) their integrative influence on clinical manifestations (Vogel et al., 2023). The major contribution of this perspective is that it identifies an intermediate phenotype (i.e., neural networks) to account for the clinical heterogeneity of neurodegenerative diseases. Moreover, it also paves the way for investigating the highly heterogeneous clinical manifestations and progression of neurodegenerative diseases using a well-developed, powerful mathematical tool—network analysis (Newman, 2003). This analysis maps distributed networks as nodes (e.g., neural structures) and edges (e.g., anatomical or functional connections between neural structures), and characterize nodal and edgewise properties of the network using graph-based descriptors (Rubinov and Sporns, 2010). Such descriptors provide a window into the rules governing the behaviors (e.g., regulatory mechanisms of nodal activities) and structural organization (e.g., connectivity between nodes) of the network. By establishing individualized and dynamic profiles of nodal and edgewise features, this network approach may offer an effective means to characterize and quantify the neurodegenerative processes at an individual level to improve personalized diagnosis and management.

So far network analysis has been primarily applied to study the degeneration of brain networks related to dementia syndromes, such as in Alzheimer's disease (Oxtoby et al., 2017; Vogel et al., 2023). The other functional domains remain relatively underexplored. Notably, bulbar motor function (i.e., speech and swallowing), which is vital to daily functioning, psychosocial, emotional, and physical health, has been rarely studied from the network perspective. Speech and swallowing disorders constitute key components of neurodegeneration in various conditions such as amyotrophic lateral sclerosis (ALS). Compared with other functional domains (e.g., mobility), there is a paucity of validated objective quantifiable markers of bulbar involvement (Pattee et al., 2019; Yunusova et al., 2019). This gap significantly hampers the early detection, progress monitoring, and tailored care of speech and swallowing disorders in neurodegenerative diseases and further undermines the efforts for therapeutic treatment discovery due to the lack of an efficacious means to evaluate treatment effects on bulbar motor function in clinical trials. Therefore, the need for objective quantifiable bulbar markers is well recognized (Chiaramonte et al., 2019; Green et al., 2018; Tondo and De Marchi, 2022). To fill this need, network analysis presents a huge potential in providing an efficient computational tool to quantify the performance of the bulbar motor system.

The purpose of this study is to exploit the potential of network analysis for objective bulbar marker development, by constructing a multiplex functional mandibular muscle network during speech to characterize bulbar involvement in ALS across the prodromal and symptomatic stages. Bulbar involvement is a hallmark feature of ALS, affecting the majority of patients during the course of ALS progression regardless of disease onset (e.g., bulbar/spinal) (Green et al., 2013; Yorkston et al., 1993). According to a patient-based subjective experience report, loss of useful speech secondary to bulbar involvement has been identified as the most devastating consequence of ALS, significantly impacting social participation, emotional health, and quality of life (Hecht et al., 2002). Although the current clinical standards for bulbar assessment focus on symptoms and functional outcomes, it is well established in the research literature that subclinical bulbar involvement starts long before the onset of clinical symptoms and functional declines in ALS (Rong et al., 2015a, 2016). Developing objective markers to capture such subclinical bulbar involvement across the prodromal and symptomatic stages can therefore provide important mechanistic insights for phenotyping the heterogeneous bulbar presentations and progression trajectories in ALS to improve personalized care.

While most existing network models are built upon brain imaging data, this study innovatively applied network analysis to electrophysiological data acquired by a clinically readily available, noninvasive instrumental technique—surface electromyography (sEMG). This methodological choice was contingent on three considerations. First, different from other neurodegenerative diseases such as Alzheimer's disease and Parkinson's disease, which primarily involve brain structures and pathways, ALS attacks various types of neurons in the cerebrum, brainstem, and spinal cord, as well as in the pathways connecting these central nervous system (CNS) structures with each other and with the muscles in the peripheral nervous system (PNS). A brain network model can only capture pathological changes in the CNS, but not in the PNS or in the interface between CNS and PNS. In contrast, the performance of a functional muscle network reflects the integrative functioning of all upstream CNS and PNS structures and pathways. Such a muscle network model can thus provide a more complete picture of the neuropathological mechanisms of ALS compared to a brain network model. Second, with immediate anatomical attachment to the effectors (e.g., speech articulators), the performance of a functional muscle network is directly associated with the behaviors of these effectors, which constitute an important part of oral mechanism exam in standard clinical bulbar assessment (Duffy, 2013; Yunusova et al., 2019). In this sense, the muscle network can offer explanatory insights into the complex, poorly understood clinical-neuropathological relationship in ALS, by providing an interface to link the clinically relevant behavioral patterns with the neuropathological underpinnings. Third, from a practical standpoint, the instrumental technique for constructing a functional muscle network (i.e., sEMG) is more accessible and less expensive compared with the neuroimaging techniques required for constructing a brain network. Lastly, the selected mandibular muscle groups for network construction have been demonstrated by previous studies to show measurable changes as early as during the prodromal stage of bulbar involvement in ALS (Rong and Jawdat, 2021; Rong and Pattee, 2021, 2022; Rong and Rasmussen, 2024). Taken together, these practical and empirical considerations support the potential of a sEMG-based functional mandibular muscle network for a clinically scalable objective bulbar assessment tool.

To construct the proposed multiplex functional mandibular muscle network, we developed a fully automated data processing and analytic approach to characterize the network performance using graph-based measures. We hypothesized that these measures would (1) effectively detect and profile bulbar involvement across the prodromal and symptomatic stages in ALS and (2) associate with clinically relevant behavioral patterns of the effector (i.e., jaw).

2 Materials and methods

This study was part of an ongoing larger-scope project. The protocol of this project was approved by the Institutional Review Board of the university medical center. Written informed consent was obtained from all participants. All study procedures were non-invasive and minimal risk. No adverse events were reported.

2.1 Participants

Fifteen individuals with ALS (nine men/six women; age: 38–77 years), including eight with overt clinical bulbar symptoms and seven without, and 10 neurologically healthy controls (three men/seven women; age: 38–81 years) participated in this study. The inclusionary criteria included: (1) being diagnosed with definite or probable ALS as per the revised El Escorial Criteria (Brooks et al., 2000) for participants with ALS or reporting no known neurological diseases or injury for healthy controls; (2) speaking American English as the first and primary language; (3) passing hearing screening at 1,000, 2,000, and 4,000 Hz at 30 dB in the better ear; (4) possessing adequate cognitive function to follow instructions and perform experimental tasks as per by standard cognitive screening procedures.

Participants with ALS were recruited from the multidisciplinary ALS clinic of the university medical center. As per standard clinical examination procedures (e.g., oral mechanism exam, clinician-based screening/evaluation, patient-reported outcomes), eight of the participants presented with overt clinical bulbar symptoms, providing samples representative of the symptomatic stage of bulbar involvement. The other seven participants exhibited no overt clinical bulbar symptoms. According to a prior longitudinal study (Rong et al., 2015a), the likelihood that these participants, who were on average 453 days post-diagnosis, had experienced subclinical bulbar involvement at the time of enrolling into this study was high despite the absence of clinical symptoms. Thus, these participants were regarded as samples representing the prodromal stage of bulbar involvement. Together, the inclusion of both prodromal and symptomatic cases offered a comprehensive sample set representative of the whole spectrum of bulbar involvement in ALS to address the research aims.

To evaluate the severity of the overall and bulbar functional disabilities, all participants with ALS completed the ALS Functional Rating Scale-Revised (ALSFRS-R) (Cedarbaum et al., 1999). ALSFRS-R is the most widely accepted means of monitoring disabilities in ALS clinical practice and clinical trials. It consists of 12 questions that assess four domains of function, including gross motor, fine motor, bulbar, and respiratory functions, with each question being rated on a scale of 0 to 4 points (4 being normal; 0 being the highest level of disability). Based on the participant's responses, the total score for all 12 questions and the bulbar subscore for the three questions related to bulbar items (i.e., speech, salivation, swallowing) were calculated to index the overall and bulbar functional disabilities, respectively. The clinical, functional, and demographic characteristics of the participants are provided in Supplementary Table S1, and a statistical summary of these characteristics is displayed in Table 1.

Table 1. Summary of demographic, clinical, and functional characteristics of participants.

2.2 Experimental task, setup, and data collection

To construct the proposed functional muscle network, the myoelectric activities of three bilateral pairs of jaw muscles—anterior temporalis (TEMP), masseter (MAS), anterior belly of digastric (ABD)—were recorded by a wireless sEMG system (BIOPAC) during a speech task. The selection of target muscles was primarily based on accessibility consideration. The morphological characteristics of anterior temporalis, masseter, and anterior belly of digastric muscles make them easily accessible and reliably measurable by surface electrodes, as demonstrated by a body of experimental work by our research team and others (Koole et al., 1991; Moore et al., 1988; Rong and Pattee, 2021).

Regarding task selection, speech production, as a fine oromotor behavior, requires rapid contractions and complex coordination of a variety of bulbar muscles. Speech tasks such as oral reading and narrative have been consistently demonstrated by prior studies to be sensitive for detecting subtle bulbar involvement (e.g., during the prodromal stage) in ALS (Rong and Heidrick, 2022, 2024; Rong and Pattee, 2022). In this study, we selected an oral reading task, where participants read a standard phonetically balanced reading passage, namely Rainbow Passage, at their habitual rate and loudness. This passage consists of 19 sentences with a total of 329 words, which cover the entire phonetic inventory with the frequency of each phoneme matching their distribution in English. Oral reading of the Rainbow Passage has proven by prior work to elicit robust activations of mandibular muscles in both neurologically healthy and impaired speakers (Rong and Jawdat, 2021; Rong and Pattee, 2022). The theoretical and empirical evidence base provides the rationale for our task selection to enhance both the robustness and generalizability of the results of this study.

To maximize reproducibility and minimize artifacts (e.g., crosstalk), we followed the best practice guidelines for surface electrode selection, setup, and recording (Castroflorio et al., 2008; Merletti and Muceli, 2019; Stepp Cara, 2012). First, the target areas of skin were prepared using an alcohol swab to increase skin conductance. Following skin preparation, self-adhesive bipolar Ag/AgCl electrodes with 11 mm diameter size and the closest inter-electrode distance possible (e.g., about 20 mm) were attached to the skin over the belly of each target muscle, parallel to fiber orientation. To enhance the reproducibility of this procedure, craniofacial anatomical landmarks were used to guide the placement of electrodes along the cantho-gonial line for masseter, vertically at the coronal suture for anterior temporalis, and submentally along the posterior-inferior direction for anterior belly of digastric. Ground electrodes were attached to the participant's shoulder. Electrode placement was verified by calibration gestures such as jaw oscillation and clenching to ensure data validity. The analog signals acquired by the surface electrodes were pre-amplified by 2,000 and band-passed filtered at 5–500 Hz by the wearable BioNormadix modules, digitized at 2,000 Hz by the MP160 module, and finally recorded by the Acknowledge software.

To examine the relation of the proposed functional mandibular muscle network to behavioral patterns of the jaw, jaw motion was recorded in three dimensions by an electromagnetic tracking system (Wave, Northern Digital Inc.). During the recording, participants were seated next to a field generator to the right; a small wired sensor was attached to the center of lower chin, and a reference sensor was attached to the center of forehead for head movement correction. Jaw motion was captured by the sensor on the lower chin relative to the forehead and was recorded at 100 Hz by the WaveFront software. Additionally, the midsagittal contour of the hard palate was traced by a manufacturer-supplied probe, from the posterior edge at the intersection between the hard and soft palate to the anterior edge at the location of the lower central incisor. This palatal trace served as an anatomical reference to characterize the motion pattern of the jaw in the subsequent kinematic analysis.

In addition to the sEMG and kinematic recordings, speech was audio-recorded by a head-mounted microphone (DPA dfine 4188) placed ~5 cm away from the left lip corner. The audio signal was preprocessed by the Behringer Xenyx 802 sound conditioner and recorded at 22,050 Hz by the WaveFront software, simultaneously with the kinematic signal. The audio recording served as a reference to assist with sEMG and kinematic data segmentation and interpretation. Figure 1 provides an overview of the experimental paradigm along with the procedures for signal processing and analysis, as elaborated below.

Figure 1. Overview of experimental paradigm.

2.3 sEMG data processing

sEMG is known to be susceptible to various sources of electrical and mechanical artifacts (e.g., motion at the skin-electrode interface, crosstalk, power line noise). To enhance the signal-to-noise ratio of the sEMG data, we used a multistep signal processing approach to minimize these artifacts, following the recommended guidelines in the sEMG literature (De Luca et al., 2010; Rong and Pattee, 2022; Stepp Cara, 2012). First, to minimize the effect of crosstalk, we adopted a blind source separation (BSS) algorithm from Kilner et al. (2002) to remove potential contamination of nearby muscles from each sEMG channel. BSS is a statistical signal processing technique to reconstruct the original unobserved signals (e.g., sources reflecting true myoelectric activities) based on the decomposition of the measured signals (i.e., surface-detected signals) (Talib et al., 2019). BSS has been successfully applied to reduce crosstalk from facial sEMG recordings in previous research (Sato and Kochiyama, 2023) and is therefore adopted by this study to serve the same purpose. After crosstalk removal, all sEMG channels were notch-filtered at 60 Hz and high-pass filtered at 20 Hz to further remove power line noise and low-frequency movement artifacts, respectively. Lastly, DC offsets were removed from each channel.

2.4 sEMG data analytics

Based on the processed sEMG signals, a weighted multiplex functional mandibular muscle network was constructed. This network consisted of six nodes, each representing a muscle; the weights of the edge between nodes reflected the functional connection between each pair of muscles in different frequency bands (see details below). To characterize the structural organization and behaviors of the network in relation to bulbar neuromuscular pathology in ALS, a fit-for-purpose data analytic program was developed and implemented in MATLAB (R2023a) to extract a variety of nodal and edgewise features from the sEMG signals for each sentence of the Rainbow Passage [total N = (25 participants × 19 sentences) − 2 errors = 473]. An overview of these network features is provided in Table 2, and the pipeline for feature extraction is outlined in Figure 2, with methodological details described below.

Table 2. Overview of network features.

Figure 2. Pipeline for constructing the multiplex functional mandibular muscle network. The network consists of six nodes corresponding to six mandibular muscles: right temporalis (RTEMP), left temporalis (LTEMP), right masseter (RMAS), left masseter (LMAS), right anterior belly of digastric (RABD), left anterior belly of digastric (LABD). Each pair of nodes are connected by an edge, with the weight of the edge (denoted by line width) reflecting the strength of functional connection between muscles (measured by intermuscular coherence) in three frequency bands (theta/alpha, beta, low gamma). Step ①: Nodal feature extraction. Using the RABD node as an example for demonstration, the surface electromyography signal for the node is first converted to a local standard deviation series, which is then transformed to a visibility graph. Based on this graph, two visibility descriptors—density and spectral radius ratio—are derived to measure the overall level and variability of nodal activity, respectively. Step ②: Edgewise feature extraction. Using the edge between the RABD and RMAS nodes as an example, the magnitude-squared coherence spectrum is calculated; after verifying significance, the mean coherence within the target frequency bands (highlighted by different colors) is computed. Using these coherence values as edge weights, a set of network measures, including mean nodal strength, assortativity coefficient, global efficiency, weighted clustering coefficient, and laterality, are derived to characterize the edgewise features of the muscle network.

2.4.1 Nodal features

Nodal features were extracted from each sEMG channel to characterize the regulatory mechanisms for modulating the myoelectric activity of each jaw muscle. To this end, each sEMG signal was first transformed into a local standard deviation series using Equation 1:

Vj=∑(j-1) * L+1jL(Ui-Ūj)2L-1    (1)

where [Uj] = [U1, U2, …, UN] is the original sEMG signal, [Vj] = [V1, V2, …, VM] is the transformed local standard deviation series, and L is the length of interval for standard deviation calculation, which was set to 50 ms (i.e., L = 100), in line with our prior study (Rong et al., 2024). These local standard deviation series were then transformed to visibility graphs. This transformation treated each data point in the series as a node and determined the connection between nodes based on Equation 2:

Vy-Vzy-z>Vy-Vxy-x    (2)

where Vz is an arbitrary node between two given nodes Vx and Vy. If all nodes between Vx and Vy meet the criterion in Equation 2, Vx and Vy are regarded as being connected. A node corresponding to a local peak in the series (i.e., a state of high-level activation) tends to have high connectivity with its neighboring nodes (i.e., high visibility). This notion relates the structure of the visibility graph with the dynamics of myoelectric activities, allowing the descriptors of the graph to inform the regulatory mechanisms for modulating myoelectric activities.

To characterize the structure of the visibility graph, two descriptors—density and spectral radius ratio—were derived. Density represents the overall node connectivity of the graph, which is defined as the ratio of the total number of the graph's edges to the largest possible number of edges, as per Equation 3:

density=2×mM(M-1)    (3)

where M and m are the total number of nodes and edges of the graph, respectively. Spectral radius ratio represents the variation of node connectivity, which is defined as the ratio of the principal eigenvalue of the adjacency matrix of the graph to the mean node degree (Meghanathan, 2014).

Together, density and spectral radius ratio provided quantitative means to assess the overall level and variability of myoelectric activities. Given the well documented changes in motor unit recruitment and firing patterns in ALS (de Carvalho et al., 2014, 2012), myoelectric activities tend to be reduced and become more irregular, as observed in prior research (Rong and Pattee, 2022). Density and spectral radius ratio were intended to target and quantify these changes related to the nodal properties of the muscle network.

2.4.2 Edgewise features

Edgewise features were calculated based on intermuscular coherence in three frequency bands—theta/alpha (4–12 Hz), beta (12–30 Hz), and low gamma (30–60 Hz). By definition, coherence is a measure of temporal correlation between two signals in the frequency domain (Laine and Valero-Cuevas, 2017; Reyes et al., 2017). It is well known that neural oscillations in different frequency bands are related to distinct neural drives for modulating muscle functions. Specifically, beta and low-gamma oscillations originate directly from the motor cortical network, which have been associated with submaximal tonic contractions and attentionally more demanding, stronger tonic and phasic contractions, respectively (Boonstra et al., 2015; Brown, 2000; Laine et al., 2015). Theta/alpha oscillations have been linked to diverse subcortical and cortical sources (e.g., brainstem, sensory cortex) outside of the motor cortical network, contributing to indirect motor control (e.g., sensorimotor integration and adaptation) (MacKay, 1997; Maezawa, 2017; Suppa et al., 2016). Therefore, measurements of intermuscular coherence in these bands can provide a window into the multiplex functional connection between muscles related to different neural sources. The feasibility of measuring intermuscular coherence in these bands from sEMG signals has been established in the literature in both neurologically healthy and impaired individuals (Carlsen et al., 2023; Fisher et al., 2012; Flood et al., 2019; Issa et al., 2017; Laine and Valero-Cuevas, 2017; Rong and Pattee, 2021).

To calculate intermuscular coherence, all sEMG signals were full wave rectified to maximize the timing information about muscle activation (Halliday et al., 1995). The rectified signals were then reconstructed by identifying and concatenating stationary 1-s epochs around the bursts for each channel. Magnitude-squared coherence spectrum was calculated based on the reconstructed signals for each pair of muscles, using Equation 4:

|Rxy|2=|Sxy(f)|2Sxx(f)Syy(f)    (4)

where |Rxy|2 is magnitude-squared coherence spectrum, Sxy(f) is the cross-spectrum between two muscles, and Sxx(f) and Syy(f) are the auto-spectra for each muscle. This analysis was implemented by a 4,096-point Fast Fourier Transform applied over a sliding 1,024-point Hamming window with 75% overlap, following the recommendations by Terry and Griffin (2008).

For each coherence spectrum, a significance level corresponding to the upper 95% confidence limit under the hypothesis of independence between muscles was calculated as: S=1-0.051/(L^-1), where L^ is the adjusted number of overlapped segments in coherence calculation (Halliday et al., 1995; Terry and Griffin, 2008). Weak, non-significant coherence usually implies a spurious connection between muscles; such a connection could obscure the topology of strong, significant connections in the muscle network (Rubinov and Sporns, 2010). Therefore, the significance level was applied as a threshold to rule out spurious connections. Within each target band (theta/alpha, beta, gamma), if there was a lack of peaks above the significance level, the coherence for the band was set to zero; otherwise, the mean coherence within the band was calculated. Finally, all coherence values were transformed to Fisher z-scores for variance stabilization (Halliday et al., 1995).

Following the procedures above, a total of 45 coherence metrics (15 edges × 3 bands) were extracted for the multiplex functional muscle network. Using these metrics as the weights of the network's edges, graph network analysis was applied to derive a set of graph descriptors to characterize the edgewise properties of the network. These descriptors included mean nodal strength, assortativity coefficient, global efficiency, and weighted clustering coefficients.

Mean nodal strength is the average sum of edge weights across nodes, which reflects the overall functional connectivity of the network. There is a mounting body of evidence showing impaired functional connection between muscles in individuals with ALS, especially in the beta band (Fisher et al., 2012; Issa et al., 2017; Rong and Pattee, 2021). Mean nodal strength was intended to detect and quantify such impairments of functional connection in ALS.

Assortativity coefficient measures the correlation of the (weighted) degrees of all nodes between two opposite sides of an edge, providing an insight into selective functional connectivity of the network. In general, an assortatively mixed network (i.e., characterized by mutually interconnected high-degree hubs) tends to have a positive coefficient, whereas a disassortative network usually has a negative coefficient. Prior studies have reported selective vulnerability of functional connection between different muscle groups in individuals with ALS (Rong and Jawdat, 2021; Rong and Pattee, 2021). Such selective vulnerability may differentially affect the connectivity between different nodes and in turn influence the assortativity of the muscle network.

Global efficiency is the average inverse shortest path length between all pairs of nodes, which provides a global index of functional integration of the network. Weighted clustering coefficient represents the connectivity of each node to its neighboring nodes, which can be interpreted as a measure of functional segregation/specialization, that is, the ability to perform specialized functions by segregated muscle groups within the network. A complex, high-performing biological system should be both functionally integrated and specialized (Tononi et al., 1998). This notion, applied to the functional muscle network in this study, implies that the activity of functionally segregated muscle groups (e.g., agonists vs. antagonists) should be integrated, in order to generate complex, coherent, and adaptive behaviors during speech. Yet, prior work has reported overall less complex, less coherent, and more irregular speech behaviors in individuals with ALS (Rong, 2021). These changes could be attributed to impaired functional integration and/or segregation of the speech motor system due to neurodegeneration. Along this line, global efficiency and weighted clustering coefficient were employed to provide complementary insights into functional integration and segregation of the muscle network.

Besides the standard graph descriptors above, another set of edgewise features was extracted to characterize the laterality of the multiplex functional muscle network. Here laterality was indexed by the difference between the weight of a specified edge on the right side and its counterpart on the left side, normalized by the sum of all edge weights in the network. This index was calculated for three types of edges, which connected temporalis with masseter, temporalis with digastric, and masseter with digastric, respectively. These indices evaluated the functional lateralization of one agonist (temporalis-masseter) and two antagonist (temporalis-digastric, masseter-digastric) muscle groups within each target band. While the majority of bulbar muscles are bilaterally innervated, a contralateral dominance has been reported in beta-band corticomuscular coherence in previous studies, suggesting that a larger proportion of corticobulbar projection comes from the contralateral hemisphere (Maezawa, 2017; Maezawa et al., 2014). For individuals with ALS, a recent evolutionary perspective of neurodegeneration has associated evolutionary hemisphere specialization with preferential involvement (Henderson et al., 2019). It posits that ALS tends to preferentially involve cerebral structures and pathways that have evolved recently in human evolution, including those in the dominant left hemisphere serving later evolved functions such as speech (Henderson et al., 2019). Together, the interplay between the disease-related asymmetrical cerebral involvement and contralateral dominance of corticobulbar projection may result in a decrease in right-to-left functional lateralization of the muscles. The laterality index was intended to capture such a change in functional lateralization.

2.5 Kinematic data processing and analysis

The 3D jaw sensor data were low-pass filtered at 15 Hz by a second-order, zero-lag Butterworth filter to remove high-frequency movement artifacts. The Euclidian distance between the jaw sensor and the lower central incisor (i.e., anterior edge of the palatal trace) was calculated at each sampled time point, generating a data series representing the global movement pattern of the jaw.

This data series was submitted to a custom-developed analysis to extract three metrics: acceleration time, mean acceleration, and stiffness. Acceleration time and mean acceleration are temporal and spatial descriptors of the acceleration profile. Scientifically, according to Newton's second law, force is the product of mass and acceleration. As such, descriptors of the jaw acceleration profile can provide insights into the force exerted by the jaw muscles on the mandible. Clinically, the slowness of orofacial movement is one of the most prevalent signs of bulbar involvement in ALS (Green et al., 2013; Yorkston et al., 1993); such a sign can manifest in the acceleration profile as increased acceleration time and reduced mean acceleration, as demonstrated by prior research (Bandini et al., 2018; Rong and Heidrick, 2021). Linking these acceleration descriptors with the network measures can inform the relation of the proposed functional muscle network to the force generation capacity of the mandibular system.

From the perspective of dynamical systems, the generation of the desired force/acceleration for a specific task requires the neuromuscular system to modulate the dynamic properties of the effectors based on task demands. Along this line, several theoretical models have been established to characterize dynamic speech behaviors with a spring-mass second-order dynamical system. These models, including the original task dynamic model proposed by Saltzman and Munhall (1989) and various extended versions (Parrell and Lammert, 2019; Simko and Cummins, 2010), posit that force/acceleration is associated with two dynamic parameters—stiffness and damping. Given that these two parameters are not mutually exclusive (i.e., damping is associated with stiffness and mass), we focused on stiffness in this study. The functional significance of stiffness is well established in the motor speech literature. Increased stiffness is associated with various mechanical advantages for generating rapid movement and/or maintaining stability against perturbation to ensure movement precision (Humphrey and Reed, 1983; Moore, 1993; Moore et al., 1988). As such, linking stiffness with the network measures can inform the relation of the functional muscle network to the dynamic control mechanism of the mandibular system.

To obtain these kinematic/dynamic metrics, the first and second derivatives of the jaw Euclidian distance data series were calculated, resulting in two new data series representing velocity and acceleration traces, respectively. Acceleration time was calculated as the average duration from the onset (i.e., velocity = 0) to the peak velocity of movement across all jaw opening and closing cycles during the speech stimulus. Mean acceleration was calculated as the average acceleration across the accelerating phase of all jaw opening and closing cycles. Unlike the acceleration descriptors, stiffness was not directly measurable from the kinematic recording but can be estimated by empirical means. In this study, stiffness was estimated by the ratio of maximum speed to maximum distance of motion (Berry, 2011).

2.6 Statistical and computational analyses

All analyses were conducted in the R statistical computing program (R Core Team, 2022). For statistical analysis, the significance level was set to p < 0.05 for main effects and was adjusted using the Bonferroni method for post-hoc tests.

2.6.1 Data reduction

As shown in Table 2, a total of 48 network features were extracted, which combined characterized the structural organization and behaviors of the multiplex functional muscle network. To reduce the dimensionality of the feature set, a data reduction technique was applied. This procedure aimed to prevent overfitting of the subsequent machine learning (ML) analyses due to the dilemma between a large number of variables and a small set of training samples, as commonly encountered in health data applications (Berisha et al., 2021). To this end, maximum likelihood factor analysis (factanal) with oblimin rotation was applied to the feature set to cluster the features that shared common variance into a lower-dimensional set of latent factors, where the number of factors was determined by parallel analysis (fa.parallel; Revelle, 2020). The fit of the factorization model was verified by chi-square goodness of fit test. After the verification, the scores of all factors were calculated using the tenBerge method. These scores represented the neuromuscular performance of the participants in the new multidimensional factor space.

2.6.2 Utility of the multiplex functional muscle network for detecting and profiling bulbar involvement in ALS

To determine the utility of the multiplex functional muscle network for detecting and profiling bulbar involvement in ALS, we first evaluated the disease effect on all network features using linear mixed effects (LME) models. For nodal features, the LME models were constructed with group (i.e., ALS vs. healthy control), node (i.e., six muscles), and group-by-node interaction as the fixed effects and a subject-dependent intercept as the random effect. Post-hoc between-group comparisons were conducted by node based on estimated marginal means (emmeans; Lenth, 2020). For edgewise features, the LME models were constructed with group, band (i.e., theta/alpha, beta, gamma), and group-by-band interaction as the fixed effects and a subject-dependent intercept as the random effect. Post-hoc between-group comparisons were conducted by band based on estimated marginal means.

In the second step, we further evaluated the efficacy of the factor scores as derived above for classification between the ALS and healthy control samples, using supervised ML algorithms. Three ML algorithms were employed, including random forest (RF), support vector machine (SVM) with the radial basis function kernel, and k-nearest neighbors (KNN), to allow for comparison of classification performance. RF is an ensemble ML algorithm that uses a combination of decision trees to solve classification or regression problems (Breiman, 2001). SVM is a flexible ML algorithm that allows raw data to be mapped into linear or nonlinear space using different kernels (Drucker et al., 1997). KNN relies on the similarity of data points (i.e., “nearest neighbors”) to assign labels (for classification) or values (for regression) (Taunk et al., 2019). All classification models were cross-validated through five-fold cross-validation repeated 10 times, and model performance was evaluated by precision (i.e., ratio of true positives to total predicted positives), recall (i.e., ratio of true positives to total actual positives), and F1 score (i.e., harmonic mean of precision and recall, defined as F1=2×Precision×RecallPrecision+Recall). Moreover, the Receiver Operating Characteristic (ROC) curve and the area under the curve (AUC) were calculated for each model.

2.6.3 Relation of the multiplex functional muscle network to clinically relevant behavioral patterns of the jaw

A stepwise approach was chosen to construct LME models for estimating the relation of the multiplex functional muscle network to the three kinematic/dynamic metrics of the jaw. For each metric, an initial LME model was constructed with group and the scores of all latent factors as the fixed effects, and a subject-dependent intercept as the random effect. Next, the fixed effects of the initial model were screened in a backward stepwise fashion to optimize the Akaike Information Criterion. The marginal and conditional R2 of the final models were calculated to evaluate model fit.

3 Results 3.1 Factorization

The factor analysis clustered the 48 network features into 10 latent factors. The rotated factor loadings are shown in Figure 3. This 10-factor model was determined to be sufficient based on chi-square goodness of fit test, χ2(693) = 4, 714.69, p = 0. The cumulative variance accounted for by the model was 68.2%.

Figure 3. Rotated factor loadings. The height of the colored bars reflects the loading of the network features on each factor. Features that load >0.3 (i.e., absolute loading >0.3) on each factor are by convention regarded as the primary component features of the factor and are marked in red; the remaining features are marked in blue. specradRatio, spectral radius ratio; nodstr, mean nodal strength; ac, assortativity coefficient; ge, global efficiency; wcc, weighted clustering coefficient; lat, laterality index. RTEMP, right temporalis; LTEMP, left temporalis; RMAS, right masseter; LMAS, left masseter; RABD, right anterior belly of digastric; LABD, left anterior belly of digastric.

Based on the rotated factor loadings, features that loaded >0.3 on each factor were conventionally identified as the primary component features of the factor. As such, the primary component features of the 10 factors were identified as: all nodal features for factor 1; gamma-band weighted clustering coefficients for factor 2; theta/alpha-band weighted clustering coefficients for factor 3; beta-band weighted clustering coefficients for factor 4; theta/alpha-band mean nodal strength and theta/alpha-band global efficiency for factor 5; laterality indices for the antagonist muscle groups (i.e., temporalis-digastric; masseter-digastric) across all three bands for factor 6; beta-band mean nodal strength and beta-band global efficiency for factor 7; laterality indices for the agonist muscle group (i.e., temporalis-masseter) across all three bands for factor 8; gamma-band mean nodal strength and gamma-band global efficiency for factor 9; assortativity coefficients across all three bands for factor 10.

3.2 Utility of the multiplex functional muscle network for detecting and profiling bulbar involvement in ALS 3.2.1 Disease effects on network features

Descriptive boxplots and statistical results of the LME models for the nodal features are provided in Figure 4 and Table 3, respectively. Note that the descriptive boxplots are displayed for each subgroup of the ALS cohort (i.e., corresponding to the prodromal and symptomatic stages of bulbar involvement) alongside the healthy controls, to allow for visual examination of stage-dependent changes in these features. For both density and spectral radius ratio, a significant main effect was found for group, but not for node, nor did the interaction between group and node show a significant effect. Post-hoc comparisons between the ALS and healthy control groups revealed a significant decrease in density and a significant increase in spectral radius ratio, for all nodes. These results implied reduced level and increased variability of myoelectric activities in individuals with ALS relative to healthy controls. Most of these disease effects were visually identifiable as early as in the prodromal stage and were incremental as bulbar involvement progressed from prodromal to symptomatic stages, as shown in Figure 4.

Figure 4. Boxplots for nodal features (density, spectral radius ratio) of the multiplex functional muscle network, displayed by node (RTEMP, right temporalis; LTEMP, left temporalis; RMAS, right masseter; LMAS, left masseter; RABD, right anterior belly of digastric; LABD, left anterior belly of digastric) and subgroup (ALSwB, individuals with amyotrophic lateral sclerosis, with overt clinical bulbar symptoms; ALSwoB, individuals with amyotrophic lateral sclerosis, absent of overt clinical bulbar symptoms; Control, healthy controls).