Introduction

Stroke is a leading cause of death and disability worldwide, the incidence of which is predicted to rise globally in coming years (Saini et al., 2021). Stroke induces motor impairments among survivors, characterised by hemiparesis (i.e. unilateral loss of function especially of the upper limb), muscle weakness, dyscoordination and fatigue (Sathian et al., 2011; Teasell & Hussein, 2016). Clinical practice primarily focuses on early and intensive intervention to maximise movement restoration (Shiromoto et al., 2017; Teasell et al., 2012), however recent evidence advocates for the importance of long-term rehabilitation (Ballester et al., 2019; Kwakkel et al., 2023; Teasell et al., 2012). Nevertheless, current clinical assessment scales don’t effectively convey comprehensive information about stroke survivors at all recovery stages, partly due to the reliance of clinically important change metrics on acute-stage recovery rates, their inability to distinguish genuine recovery from compensation or to capture participation-level outcomes (Ballester et al., 2019; Bushnell et al., 2015; Page et al., 2012). Hence, a current research gap lies in effective assessment tools that effectively quantify post-stroke motor impairment and treatment responsiveness at all recovery stages.

Muscle synergy analysis is a data-driven and theoretically supported approach from the motor neurosciences that quantifies coordinated patterns of muscle activity (‘muscle synergies’) from surface electromyographic (sEMG) signals (Berret et al., 2019; Bruton & O’Dwyer, 2018; Cheung & Seki, 2021), interpreting these patterns as the output of functionally modular neural networks in the central nervous system. Through this approach, people affected by stroke have been shown to exhibit different neurophysiological responses depending on impairment severity and adaptations during recovery, including the preservation, merging and fractionation of unimpaired motor patterns along with the generation of novel motor patterns (Cheung et al., 2012; Clark et al., 2010; García-Cossio et al., 2014; Hashiguchi et al., 2016; Hong et al., 2021; Pregnolato et al., 2025). Indeed, useful biomarkers have been produced from this approach that convey linear relationships with clinically validated measures of motor impairment (Clark et al., 2010; Schwartz et al., 2016). However, restrictive model assumptions (e.g. linearity) and the lack of a direct mapping of muscle interactions to task performance has hampered the progression of this approach in practical applications (Alessandro et al., 2013; Berret et al., 2019; Cheung & Seki, 2021; de Rugy et al., 2013). In our recent work (O’Reilly & Delis, 2024; O’Reilly & Delis, 2022), we addressed these important limitations by incorporating network- and information-theoretic approaches into a novel, generalised framework (Network-Information Framework (NIF)) that maps muscle interactions directly to their functional consequences (Fig.1(A)). In doing so, we can more precisely quantify functional muscle co-variations, distinguishing between task-relevant (i.e. muscle activation patterns associated with task-specific movements (pink-orange chequerboard intersection (Fig.1(A))) from task-irrelevant (i.e. muscle activation patterns present across motor tasks (yellow intersection exclusively shared by m_x and m_y (Fig.1(A)) contributions intermixed in conventional approaches and going further in characterising their functional roles as functionally-similar (i.e. redundant (pink in pink-orange chequerboard intersection (Fig.1(A)) and -complementary (i.e. synergistic (orange in pink-orange chequerboard intersection (Fig.1(A))). Hence, the NIF offers novel research opportunities and practical use cases for muscle synergy analysis that may be beneficial in addressing these current limitations in the clinical assessment of the stroke population.

Figure 1.

Therefore, in the current study we applied the NIF to a cohort of 42 stroke survivors performing 7 standardised upper-limb movements pre- and post-20 sessions of intensive upper-limb motor training through either a virtual reality (VR) or conventional physical therapy (PT) intervention (Fig.1(B)). In doing so, we aimed to provide novel characterisations and biomarkers of post-stroke impairment and therapeutic responsiveness that address these current limitations in post-stroke clinical assessments. To do so, we firstly extracted redundant and synergistic modules that reflected the diverse biomechanical contributions of collective muscle interactions towards the 7 standardised motor tasks and performed an in-depth analysis on their underlying activation patterns (Fig.1(C-D)). Through our novel clustering algorithm (Fig.1(E)), we identified patient clusters aligned with and accentuating established clinical measures of motor impairment (moderately and severely-impaired subgroups) and treatment responsiveness (responders and non-responders subgroups). Finally, by probing the underlying contributions of merging, fractionation, and preservation patterns to the observed patient clusters, we shed mechanistic light on their role in the manifestation of movement pathology. Taken together, this study provides an in-depth functional characterisation of post-stroke impairment and treatment responsiveness whilst progressing the clinical use case for muscle synergy analysis as an independent assessment tool for post-stroke rehabilitation.

Results

The cohort of stroke survivors overall experienced a statistically significant increase at FMA-UE (Pre-treatment: 43.1±13.2, Post-treatment: 49.1±13.6 (t= −7.84, p < 0.001)), representing a clinically important effect from rehabilitation (Page et al., 2012). Although the VR group displayed a greater improvement in FMA-UE scores compared to PT, no significant differences between treatment conditions were found at baseline (PT: 42.75±16.3, VR: 43.2±12.6 (t= −0.69, p > 0.9)) or at follow-up (PT: 45.62±15.2, VR: 49.9±13.3 (t= −0.723, p > 0.45)). 21 responders and 21 non-responders were iden-tified using the conventional MCID of >5 points on the FMA-UE (responders (PT=3/8, VR=18/34), non-responders (PT=5/8, VR=16/34)) (Page et al., 2012). No significant differences were found between VR and PT responders and non-responders at baseline or follow-up (p > 0.05).

Functional muscle networks reflect the biomechanical contributions of collective muscle interactions

The application of the NIF across the post-stroke cohort and both pre- and post-sessions identified five redundant (R1-R5 (Supplementary Materials Fig.1)) and seven synergistic (S1-S7 (Supplementary Materials Fig.2)) networks of co-occurring muscle interactions across affected and unaffected sides. These functional muscle networks were comprised of muscle interactions that contributed in functionally-similar and -complementary ways respectively towards the seven standardised motor tasks performed (see Table.2 of Materials and Methods section for details). Subsequent examination of the muscle interaction strengths (relative network edge width), principal muscles (relative node size) and submodular structure (node colour) enabled a biomechanical function to be attributed to each extracted component (see the supporting information for Supplementary Materials Fig.1-2 for detailed descriptions and ‘Extraction of Redundant and Synergistic Muscle Networks’ section of the Materials and Methods section for technical details). Muscle acronyms are described in the caption of Fig.1 here.

To briefly summarise, the post-stroke cohort were characterised by multiple extensor- and flexor-driven muscle networks comprised of predominantly redundant and synergistic muscle interactions respectively. The redundant muscle networks specifically were also consistently composed of functional contributions of shoulder girdle and/or elbow stabilisers (e.g. TEM, Infra, PT and SU) supporting prime-mover musculature (e.g. AD, TL and TM). There was a functional correspondence between the affected- and unaffected-side redundant networks (e.g. R2 and R5 (Supplementary Materials Fig.1)), however, noticeable differences relating to the underlying impairment of the affected-side across the cohort were also represented (e.g. R3 (Supplementary Materials Fig.1), involving interactions between SU and the shoulder musculature on both sides, displayed a specific reliance on PM on the affected-side whereas for the unaffected-side, the interactions were spread across PM, AD and LD, representing the characteristic loss of muscle selectivity at the shoulder level post-stroke). Meanwhile the synergistic muscle networks were generally characterised by widespread dependencies between distal and proximal musculature (e.g. S6 (Affected-side) Supplementary Materials Fig.2) that were less comparable across affected- and unaffected-sides. Contributions from joint stabilising musculature such as TEM and Infra were also found among synergistic networks, suggesting they contribute in both functionally-similar and -complementary ways to task performance with respect to prime-mover muscles.

Network recruitment patterns are associated with motor impairment, treatment type and responsiveness

In assessing the statistical relationships between the recruitment of specific redundant and synergistic network components and functional impairment, we found that the degree of recruitment of S2 and S6 (both comprised of complementary muscular contributions to elbow flexion with shoulder girdle stabilisation (Fig.2(A)) and Supplementary Materials Fig.2) at baseline was associated with lower FMA-UE scores both at baseline (β = −174.5±62.5 (p = 0.008), β = −159.7±45.6 (p = 0.001) respectively) and at follow-up (β = −172.5±65.4 (p = 0.012), β = −155.1±48.2 (p = 0.003) respectively) (Fig.2(A)). This suggests that specific networks of synergistic muscle interactions are strongly linearly related to functional impairment at the group-level in such a way that persists after clinical intervention. Meanwhile, no significant relationships were found between individual activation coefficients and FMA-UE scores after treatment, suggesting that perhaps the activation coefficients after treatment reflected features specific of participant and modality provided.

Figure 2.

Examination of the extracted recruitment patterns relationship with treatment responsiveness among responders and non-responders and among PT and VR participants revealed further evidence for the crucial role of synergistic muscle interactions (Fig.2(B)). Specifically, we found that among non-responders S3 recruitment changed significantly after treatment than responders (t = 3.1 (p = 0.004)), while the VR group demonstrated a significant change in the recruitment of S4 compared to participants in the PT group (t = −2.2 (p = 0.044)). Interestingly, both S3 and S4 were functionally related to elbow flexion with shoulder girdle stabilisation (Fig.2(B) and Supplementary Materials Fig.2), suggesting that non-responders express maladaptive recruitment adjustments of specific flexor-driven networks, while VR treatment promotes changes in other flexor-driven networks.

A redundancy-to-synergy transformation characterises effective rehabilitation

To characterise treatment responsiveness from the perspective of the NIF, we took the redundant and synergistic muscle networks of each participant on both affected- and unaffected-sides and statistically compared muscle network-level (i.e. average network interaction strengths) and interaction-level differences (i.e. individual muscle connection strengths (see ‘Characterising reha-bilitation effects on functional muscle interactions among responders and non-responders’ section of Materials and Methods)) between sessions among conventionally-defined responders (i.e. >5 point change in FMA-UE scores) (Fig.3(A-B.1-2)) and non-responders (i.e. <5 point change in FMA-UE scores) (Fig.3(C-D.1-2)).

Figure 3.

Treatment responders experienced a significant decrease in redundancy on the affected-side (t= 2.3 (p < 0.05) (Fig.3(A.1)) and significant increase in synergy on the affected-side also (t= −2.2 (p < 0.05) (Fig.3(B.1)). No noticeable changes were experienced for average redundancy or synergy on the unaffected-side across responders were observed (Fig.3(A-B.1)). Specific muscle interactions also reflected this trend for redundancy and synergy with rehabilitation among responders (see network edges in Fig.3(A-B.2)). Whereas for the corresponding muscle interactions on the unaffected-side (Fig.3(A.2)), a mixture of significantly greater and lower magnitudes from baseline were found. The most significant synergistic muscle interactions of responders’ affected-side were also predominantly greater at follow-up (Fig.3(B.2)), except for BRACH-SU which was lower at follow-up. This was also the case for the unaffected-side of responders (Fig.3(B.2)), where TEM-PT was the only muscle interaction of those depicted to be lower at follow-up.

Meanwhile non-responders demonstrated no significant changes in average redundancy or synergy from pre-to-post sessions on either affected- or unaffected-sides (Fig.3(C-D.1-2)). Several redundant and synergistic muscle interactions on the unaffected-side were found to increase from pre-to-post sessions (Fig.3(C-D.2)), suggesting rehabilitation may have influenced control of the unaffected-side among non-responders rather than the affected-side.

Altogether, these findings characterise post-stroke motor recovery as a redundancy-to-synergy information conversion, demonstrating quantitatively the enhanced functional differentiation of muscle interactions provided by effective rehabilitation.

Data-driven patient clusters more precisely quantify treatment response patterns

Implementation of our clustering approach (see Fig.5 and ‘Clustering stroke survivors based on impairment severity and therapeutic (non-) responsiveness’ section of the Materials and methods for details) to quantify treatment response patterns (i.e. pre-to-post session changes in modular recruitment patterns) identified two redundant and synergistic patient subgroups (R_Pre-Post, S_Pre-Post). Supplementary Materials Fig.3 illustrates these clusters which did not demonstrate significant differences (p > 0.05) for pre- and post-session FMA-UE scores and the change in FMA-UE scores across the intervention (ΔFMA-UE). However, direct comparison of R_Pre-Post and S_Pre-Post clusters to the FMA-based responsiveness classification in terms of network- and interaction-level differences across the clinical intervention (see ‘Characterising rehabilitation effects on functional muscle interactions among (non-) responders’ section of Materials and methods for further details) demonstrated a more precise quantification of treatment responsiveness.

More specifically, for R_Pre-Post (Cluster 1 (Fig.3(E.1-2)), Cluster 2 (Fig.3(G.1-2)), we found that cluster 1 participants significantly reduced redundancy on the affected-side from pre-to-post sessions (t= 3.5 (p < 0.001) (Fig.3(E.1)), while the unaffected-side experienced no significant changes. Several of the affected-side muscle interactions were found to decrease in strength at post-treatment, while the unaffected-side demonstrated mostly increased with rehabilitation (Fig.3(E.2)). Mean-while, cluster 2 participants demonstrated a non-significant increase in the average redundancy of the affected-side and no perceived changes on the unaffected-side (Fig.3(G.1)). S_Pre-Post (Cluster 1 (Fig.3(F.1-2)), Cluster 2 (Fig.3(H.1-2)) demonstrated a significant increase in synergistic muscle interaction strengths among cluster 2 participants (t= −3.78 (p<0.001) (Fig.3(H.1))), while no significant changes were found on the unaffected-side. The most significantly different synergistic muscle interactions on both affected- and unaffected-sides were all found to be greater at post-treatment (Fig.3(H.2)). On the other hand, for cluster 1 participants (Fig.3(F)), we found a non-significant decrease in synergy on the affected-side and no perceivable changes on the unaffected-side (Fig.3(F.1)). The most significant synergistic muscle interactions were mixed in both directions for cluster 1 participants affected-side and all greater at follow-up for the unaffected-side (Fig.3(F.2)).

Taken together, we found that the NIF-derived treatment response clusters not only closely aligned with but further accentuated the redundancy-to-synergy conversion underpinning conventional metrics (see ‘Effective post-stroke rehabilitation as a redundancy-to-synergy information conversion’ section of the Results), suggesting a more precise quantification of treatment responsiveness.

Stroke survivors cluster into severely and non-severely impaired subgroups with distinct physiological responses to neurological insult

Separate cluster analysis of pre- and post session modular recruitment patterns across the post-stroke cohort (see Fig.5 and ‘Clustering stroke survivors based on impairment severity and therapeutic (non-) responsiveness’ section of the Materials and methods for details) consistently revealed two redundant and synergistic patient subgroups at baseline (i.e. R_Pre and S_Pre (Fig.4(A-B))) and at follow-up (i.e. R_Post and S_Post (Supplementary Materials Fig.4(A-B))). Table.1 below provides an overview of the patients assigned to each cluster identified, showing participants across treatment groups (VR vs PT) and both clinical assessment- and NIF-defined responders (R) and non-responders (NR) were present in each cluster.

Table 1.

Figure 4.

Figure 5.

A significant difference (*) between clusters was found for R_Pre for both pre-session FMA-UE (t= 2.27 (p = 0.029) (Fig.4(A.2))) and post-session FMA-UE (t= 2.37 (p = 0.023) (Fig.4(A.3))), but not for the change in FMA-UE score from pre-to-post sessions (ΔFMA-UE) (Fig.4(A.4)). The S_Pre partition also significantly differentiated participants based on motor impairment at baseline (t= −2.9, p < 0.01(**) (Fig.4(B.2))) and follow-up sessions (t= −2.5, p = 0.017(*) (Fig.4(B.3)), but similarly to R_Pre, could not distinguish participants based on ΔFMA-UE (Fig.4(B.4)). For both S_Post and R_Post, presented in Supplementary Materials Fig.4(A-B) respectively, the extracted partitions could not significantly differentiate participants in terms of pre- or post-session FMA-UE scores or ΔFMA-UE.

To shed mechanistic light on the underlying contributions to the patient clustering’s observed in Fig.4(A-B) and Supplementary Materials Fig.4(A-B) here, we sought to determine if they captured distinct pattern of merging, fractionation and preservation among the post-stroke cohort (see ‘Quantifying the contributions of preservation, merging and fractionation to impairment-based patient clusters’ section of Materials and methods for details on their quantification). The network components and their corresponding physiological responses found to be significantly different between R_Pre and S_Pre patient clusters are illustrated in Fig.4(C-D) while the corresponding outputs for R_Post and S_Post are presented in Supplementary Materials Fig.4(C-D). For interpretation, along-side each significantly contributing module, we have illustrated a corresponding network from a representative participant involved in each calculation.

To summarise, we found that for R_Pre (Fig.4(C)), the prevalence of fractionation in R3 (comprised of interactions among the shoulder prime-movers (i.e. PD, AD and LD)) was associated with the less impaired cluster 1 participants (β = −5.84±2.23 (p < 0.01), 78.6% classification accuracy). This evidence suggests that the fractionation of specific redundant modules is an adaptive response to neurological insult among non-severely impaired stroke survivors.

Fractionation of S3 along with merging of S1 at baseline contributed to the patient clustering’s observed in the more severely impaired subgroup of S_Pre (Fig.4(D.1-2)), with both associated with cluster 1 patients (β = −34.82±10.9 (p = 0.001), β = −16.4±8.6 (p = 0.056) respectively (Classification accuracy= 83.3%)). The fractionation of muscle groups involving the shoulder girdle (i.e. PM (Fig.4(D.1))) was again an important contributor to the observed patient groupings, however in contrast to the redundant modules (Fig.4(C)), here fractionation of synergistic modules was associated with the severely impaired patient cluster. The first instance of merging patterns contributing to patient subgroupings here was found for S_Pre, demonstrating how the maladaptive merging of elbow flexor related interactions (i.e. BL-BS (S1) (Fig.4(D.2)) underpins the more severely impaired stroke survivors.

Fractionation of R4 (centralised around redundant interactions between UT and both scapula stabilisers and wrist extensor muscles) at follow-up was also associated with the cluster 1 participants of R_Post (β = −2.32±0.99 (p = 0.019) (Supplementary Materials Fig.4(C.1))). The underlying contributions to R_Post clusters were more complex than R_Pre however, including an additional contribution by the preservation patterns in R1 (an elbow extensor related functional module) (β = −15.2±6.5 (p = 0.018)) (Supplementary Materials Fig.4(C.2))). Together, these contributions illustrate that cluster 1 participants of R_Post had significantly greater fractionation of R4 and preservation of R1 compared to cluster 2 participants at follow-up.

Finally, patterns of fractionation (S1 (Supplementary Materials Fig.4(D.1)) and of merging (S2 (Supplementary Materials Fig.2(D.2))) were found to significantly contribute to the S_Post clusterings, both associated with cluster 2 affiliation (β = 26.7±12.1 (p = 0.027) and β =27.3±11.8 (p = 0.021) respectively (Classification accuracy: 76.2%)). S1 here, like S3 in S_Pre (Fig.4(D.1)), comprised of PM playing a central role in interaction with all other muscles (Supplementary Materials Fig.4(D.1)). There wasn’t however a clear correspondence between the contributing merging patterns for S_Pre and S_Post (Fig.4(D.2) and Supplementary Materials Fig.4(D.2) respectively). Instead, for S_Post the merging pattern significantly differentiating participants at post-treatment predominantly involved the SU.

Discussion

This study leveraged a novel muscle network analysis framework to provide a data-driven account of post-stroke motor impairment and therapeutic responsiveness and, in doing so, address current challenges to the standardised assessment of stroke survivors. Multiple flexor- and extensor-driven modules were uncovered from a post-stroke cohort undergoing intensive upper-limb motor training, collectively representing the diverse biomechanical contributions of functional muscle networks to task performance. Their recruitment patterns were closely associated with motor impairment and both treatment type and responsiveness, promoting their use as clinical biomarkers. Deployment of a novel clustering algorithm revealed consistent dichotomous patient subgroupings reflecting clinical assessment-based classifications of motor impairment (i.e. severely vs. non-severely impaired patients) and treatment responsiveness (i.e. responders and non-responders) underpinned by distinct physiological markers of impairment (i.e. muscle network merging, fractionation, and preservation). Crucially, we provided a novel characterisation of effective post-stroke rehabilitation as a redundancy-to-synergy information conversion that highlights the enhanced functional differentiation of muscle interactions with effective treatment, a response pattern closely reflected in conventional clinical assessment-based classifications and further accentuated in the derived patient clusters. Our work here provides novel neurobiological insights into and biomarkers of movement pathology whilst offering a principled approach for precision motor assessments in stroke survivors.

Altered composition and recruitment of functional muscle networks post-stroke

The functional characterisation of muscle interactions as redundant (i.e. functionally-similar) and -synergistic (i.e. functionally-complementary) is a unique capability of the NIF (O’Reilly & Delis, 2024). By extracting these distinct types of functional interaction, we uncovered multiple task-relevant extensor- and flexor-driven synergies that, remarkably, were found predominantly among redundant and synergistic muscle networks respectively (see Supplementary Materials Fig.1-2). These findings align with recent work demonstrating multiple modules underpin the post-stroke flexor synergy (Kim et al., 2024a), going further in showing that multiple modules also underpin the post-stroke extensor synergy and that both synergies are comprised predominantly of distinct types of functional interaction. Continuing, the activation of several flexor-driven synergistic muscle net-works closely reflected motor impairment consistently across the intervention and both treatment-type and -responsiveness (Fig.2), exemplifying the central role of the flexor synergy in manifestations of functional impairment post-stroke (Alt Murphy et al., 2022; Ellis et al., 2017). We present the NIF here as a principled methodology for the separate quantification of abnormal flexor- and extensor-synergy expression along with their distinct contributions to motor function and treatment responsiveness.

Impaired control of not only prime-mover muscle groups but also secondary, stabilising musculature are prevalent post-stroke (Paci et al., 2007; Ryerson et al., 2008; Silva et al., 2018). In contrast to traditional muscle synergy analysis (Bruton & O’Dwyer, 2018; Cheung & Seki, 2021; Turpin et al., 2021), the NIF does not rely on optimising data variance explained and therefore the functional contributions of more subtle interactions often marginalised in traditional muscle synergy approaches were equivalently emphasised alongside higher amplitude muscle covariations (see Supplementary Materials Fig.1-2) (O’Reilly & Delis, 2024; Ó’Reilly & Delis, 2022). This additional capability revealed widespread redundant and synergistic interactions among secondary stabilising musculature (e.g. teres major, infraspinatus and rhomboids) and with prime-mover muscles (e.g. triceps and biceps brachii, deltoids) that frequently demonstrated high network prominence and interaction strengths (see Tables in Supplementary Materials Fig.1-2). Hence, the NIF offers novel clinical utility to the assessment of these subtle muscle interactions underlying motor functions commonly impaired post-stroke such as anticipatory postural adjustments and segmental/joint stabilisation (Paci et al., 2007; Silva et al., 2018).

Table 2.

The fractionation of unaffected motor patterns into separate modules in response to neurological insult emerged as a central difference between the patient clusters identified here (see Fig.4 and Supplementary Materials Fig.4). We found new evidence for the adaptive role of fractionation among redundant muscle interactions but maladaptive role of fractionation among synergistic muscle networks (Fig.4). This provides evidence towards crucially unresolved questions on the role of muscle synergy fractionation from the original work presenting these physiological markers (Cheung et al., 2012), demonstrating that it indeed depends on the type of functional muscle interaction involved. Further, the NIF demonstrated here an enhanced capability over traditional approaches to identify these crucial patterns, as earlier work on related versions of this dataset could not identify any differentiable fractionation events across the cohort (Pregnolato et al., 2025). Although fractionation was found to also differentiate patient clusters at follow-up (Supplementary Materials Fig.4(C-D)), the exact role fractionation played in treatment responsiveness was unclear here as no significant differences between post-session patient clusters in terms of FMA-UE were observed (Supplementary Materials Fig.4(A&B.3). Finally, the observations of both the merging and preservation of motor patterns contributing to motor impairment and perhaps treatment responsiveness aligns with the existing literature (Cheung et al., 2012; Clark et al., 2010; García-Cossio et al., 2014; Hashiguchi et al., 2016), and offers the potential for additional insight into their frequently contrasting effects on motor function depending on the individual’s condition—whether it’s impairment severity, developmental stage, or prior training (Cheung et al., 2020; García-Cossio et al., 2014).

Post-stroke motor recovery as a redundancy-to-synergy information conversion

In this study, we provided a novel characterisation of treatment-induced motor recovery among stroke survivors as a redundancy-to-synergy conversion of task information (Fig.3). Specifically, a reduction in redundancy and increase in synergy on the affected-side of responders was found while no such changes were observed among non-responders or on the unaffected-side. This pattern was not only observed in the functional muscle networks of patients grouped in the conventional way using the MCID of the FMA-UE measure, but also in the NIF-derived patient clusters (i.e. R_Pre-Post, S_Pre-Post) where it was in fact further accentuated. Enhanced synergy with treatment here reflects the functional differentiation of predominantly flexor-driven muscle networks towards different, complementary task-objectives across the seven upper-limb motor tasks performed (Kim et al., 2024b), leading to improved motor function among responders. Although the exact neural substrates underpinning redundant and synergistic muscle interactions remain a future research matter, abnormal muscle interactions post-stroke are known to originate from damage to the corticospinal tract (CST), particularly effecting flexor muscle groups and which are compensated for by reticulospinal tract (RST) upregulation (García-Cossio et al., 2014; Kim et al., 2024a; McPherson et al., 2018). Both neural pathways are known to exhibit specific flexor-extensor biases (Baker, 2011), but both indeed provide parallel input to motoneurons to coordinate groups of muscles towards complementary functional objectives (e.g. fine vs. gross motor functions) (Glover & Baker, 2022; Sangari & Perez, 2020). Hence, treatment-induced enhancement of functional complementarity among predominantly flexor-driven muscle networks here could reflect CST functional recovery and therefore, a reduced reliance on alternative neural pathways (Kim et al., 2024a, Kim et al., 2024b; McPherson et al., 2018). These observations promote the clinical use case of NIF in the precise assessment and monitoring of neural impairment over traditional approaches (O’Reilly & Delis, 2024; O’Reilly & Delis, 2025; O’Reilly & Delis, 2024).

Dichotomous clustering of motor impairment and treatment responsiveness in the post-stroke population

Through the novel clustering of stroke survivors presented here, we consistently identified dichotomous subgroupings that reflected distinct patterns of motor impairment and treatment responsiveness (Fig.3-4 and Supplementary Materials Fig.3-4). These clustering’s are formed the greater affinity among some stroke survivors compared to others in terms of their modular activations and changes thereof with treatment. Based on our own and others’ previous work (García-Cossio et al., 2014; O’Reilly & Delis, 2024; O’Reilly & Delis, 2025), this functional grouping results from the preservation of unimpaired modules among certain cohorts and the generation of new motor patterns in others, a distinction seemingly dependent on residual CST function (García-Cossio et al., 2014). These findings link well also with the increasing awareness in the rehabilitation literature for the distinct aetiology of severe and non-severe impairments in stroke survivors that is not merely an artefact of measurement or unfounded simplification (Bonkhoff et al., 2022; Byblow et al., 2015; García-Cossio et al., 2014; Kim et al., 2024a; O’Reilly & Delis, 2024). This awareness comes from the distinct rates of proportional recovery observed among stroke survivors with and without an intact CST (Bonkhoff et al., 2022; Byblow et al., 2015). These observations do not rule out the possibility of more fine-grained patient subgroups, as preliminary applications without sparsification and related work employing the clustering algorithm employed here revealed 5 further subclusters, aligning with previous research (O’Reilly & Delis, 2025; Scano et al., 2017; van der Vliet et al., 2020). Nonetheless, the principled quantification of more generally applicable yet crucially discernible patient clusters offers new opportunities in the rehabilitation of severely impaired stroke survivors, a cohort that experiences significantly greater challenges during rehabilitation (Angerhöfer et al., 2021; McGlinchey et al., 2020).

Methods and Materials

Experimental setup and Data collection

We conducted a secondary analysis on previously published data (Maistrello et al., 2021; Pregnolato et al., 2025) (ethical approval provided by the local Ethical Committee of the IRCCS San Camillo Hospital s.r.l, Venice, Italy), including 42 stroke survivors (Age (Years): 62.4 ±12.2, Gender (Male/Female): 28/14, Diagnosis (ischaemic/haemorrhagic): 36/6, Time from lesion onset (Months): 7.5±14.1). All participants had previously been hospitalised with their first stroke. Exclusion criteria for this study included: 1) Mini Mental State Examination (MMSE) score lower than 20 points indicating moderate cognitive impairment; 2) severe verbal comprehension deficits, defined >13 errors on the Token Test; 3) evidence of apraxia and visuospatial neglect, evaluated by neurological examination; and 4) history of behavioural disorders (e.g., depression, aggressiveness, apathy) or neurological/vascular co-morbidity. This pre-post longitudinal study involved 4 weeks of intensive treatment consisting of 20 sessions of upper-limb exercises (one hour/session, 5 sessions/week). Participants that performed less than 80% of the planned sessions (< 16/20 sessions) were considered as study drop-outs and not included in further analyses. Two treatment groups were defined for this study which participants were randomly assigned into either conventional physical therapy (PT, No. participants enrolled= 8) or virtual reality (VR, No. participants enrolled=34) intervention (VRRS®, Khymeia Group Ltd., Noventa Padovana, Italy) groups (see the parent papers of this experimental protocol for further details (Maistrello et al., 2021; Pregnolato et al., 2024) (Fig.1(A-B)). The Fugl-Meyer Assessment Upper Extremity section (FMA-UE) was used to quantify motor impairment in the affected limb (Meyer et al., 1975), both prior to and immediately following the 4-week intervention. A minimal clinically important difference (MCID) of 5 points change at FMA-UE scores after treatment was used to classify responders and non-responders in both groups (Page et al., 2012). Moreover, before and after treatment an instrumental assessment was conducted recordings sEMG (EMG-USB2+, OT Bioelettronica, sampling rate: 2000Hz) from 16 muscles (triceps -medialis (TM) and lateralis (TL); biceps short head (BS) and biceps long head (BL); anterior deltoid (AD); lateral deltoid (LD); posterior deltoid (PD); upper trapezius (UT); rhomboid major (RM); brachioradialis (BR); supinator (SU); brachialis (BRACH); pronator teres (PT); pectoralis major (clavicular head) (PM); infraspinatus (Infra); teres major (TEM) (Fig.1(C))) of both affected and unaffected sides, while participants performed 10 repetitions of seven standardised motor tasks (Table.2). To prepare the sEMG signals for subsequent analyses, they were full-wave rectified and low-pass filtered at 20Hz using a 4th-order Butterworth filter with zero-phase distortion.

Extraction of Redundant and Synergistic Muscle Networks

To quantify the task-relevance of muscle interactions and then characterise them as having either functionally-similar (i.e. redundant) or -complementary (i.e. synergistic) roles, we employed the NIF to map networks of muscle interactions to task performance and extract the underlying low-dimensional components (O’Reilly & Delis, 2024). More specifically, leveraging a higher-order information-theoretic measure known as co-information (II) (Ince et al., 2017; McGill, 1954) (Fig.1(A)), we firstly quantified the task-relevant information (pink-orange intersection) shared between all possible muscle pairings (e.g. [m_x, m_y]) with respect to a given task parameter (τ), formally disentangling it from the task-irrelevant space (yellow intersection) (Alessandro et al., 2013).

II quantifies the balance between redundancy (i.e. the task-specific variations that are equivalently shared by m_x and m_y and so provide less information when combined) and synergy (i.e. the task-specific variations of m_x and m_y that when observed together result in additional predictive information about τ) by quantifying the difference between the task information shared (i.e. the mutual information (Ince et al., 2017)) provided by m_x and m_y separately, I(m_x; τ) + I(m_y; τ), and shared together, I(m_x, m_y; τ). Consequently, the functional role of each muscle interaction can be principally quantified as net redundant (negative II; pink shading of pink-orange intersection) or synergistic (positive II; orange shading of pink-orange intersection). These distinct quantities generalise the conceptual underpinning of muscle synergies (Bruton & O’Dwyer, 2018), formally distinguishing between two alternative ways in which muscles may ‘work together’ towards task performance.

For this study’s analysis, we generated redundant and synergistic networks from the affected and unaffected side muscles separately with respect to a discrete task parameter describing the seven motor tasks each participant performed. We carried out this application in the spatial domain (i.e. interactions between muscles across time (Ó’Reilly & Delis, 2022)) by concatenating the 10 repetitions of each task executed on a particular side and quantifying II with respect to this discrete task parameter codified to describe the motor task performed at each timepoint within each trial. Where a muscle interaction was determined to be net redundant/synergistic, their corresponding network edge in the other muscle network was set to zero. Having quantified these distinct types of task-relevant muscle networks, the following steps in the NIF pipeline were then run:

1. We empirically sparsified each network by identifying robust muscle network structures not found in equivalent random networks using a modified percolation analysis (Gallos et al., 2012), setting non-significant muscle interactions to zero.

2. We then determined the functional organisation of each participant’s muscle network using a link-based community detection method (implemented with single linkage) (Ahn et al., 2010) (Fig. 2(C.1)), essentially unravelling the nested structure of redundant or synergistic interactions into separate co-membership matrices for each module identified (i.e. entries indicate where a pair of nodes are functionally affiliated or not) (O’Reilly et al., 2025; O’Reilly & Delis, 2024). A consensus partition across participants was then subsequently determined by aggregating across all of these adjacency matrices and applying the Louvain algorithm to the resultant adjacency matrix (Blondel et al., 2008; Rubinov & Sporns, 2010).

3. The number of muscle network modules identified in this consensus partition was used as the input parameter for dimensionality reduction, namely projective non-negative matrix factorisation (PNMF) (Fig. 1(D)) (Yang & Oja, 2010). The input matrix for PNMF consisted of both affected and unaffected sides redundant/synergistic muscle networks from all participants at both pre- and post-sessions concatenated together in their vectorised forms. More specifically, the input matrix was configured such that the set of unique pairings (1 … K) of affected and unaffected muscles (m_aff and m_unaff respectively) corresponding to an individual participant (p) and session from the group of participants (1… N) are concatenated row-wise and then the set of pre-session (pre) and post-session (post) networks are concatenated column-wise (Equation 1.1). Factorisation of this matrix results in the extraction of, for the jth module, a representative set of co-occurring muscle interaction weightings (υ) across affected and unaffected sides and their underlying session- and participant-specific activation coefficients (𝒶).

The extracted components can be reshaped into fully connected adjacency matrices (i.e. muscle networks) where the magnitude of the edges represents the interaction strength between muscle pairs. For the purpose of illustrations, the most prominent muscle interactions (i.e. proportional thresholding (p<0.05)) within each extracted module are presented as a network over a human body model (Rubinov & Sporns, 2010). To visualise the extracted networks’ submodular structure, we applied the Louvain algorithm and depicted the partitions as distinct node colours on a human body model (Blondel et al., 2008; Rubinov & Sporns, 2010). Moreover, to illustrate the principal muscles involved in each component, the interaction strengths and relative importance of individual muscles were used to indicate edge-widths and node sizes respectively. Relative importance of muscles was quantified via the eigenvector centrality (i.e. the eigenvector associated with the largest eigenvalue in the network) (Newman, 2008), where muscles were considered prominent when they were well-connected with other well-connected muscles. Matlab code for this analytical pipeline has been published here: https://github.com/DelisLab/EMG2Task.

Clustering stroke survivors based on impairment severity and therapeutic responsiveness

To understand how the stroke survivors compared in terms of motor impairment severity and treatment-responsiveness using NIF, we conducted cluster analyses using the extracted activation coefficients (𝒶). Towards finding an effective approach to clustering participants in this data based on differences in impairment severity and therapeutic (non-)responsiveness, we found that conventional clustering algorithms (e.g. agglomerative, k-means etc.) could not provide substantive outputs, perhaps resulting from the complex interdependencies between the modular activations. Therefore, we sought to improve the linear separability of the participants by mapping 𝒶 into a higher-dimensional feature space, and then applying network community detection protocols already incorporated into the NIF here to discern the underlying population structure (Fig.5) (O’Reilly & Delis, 2025). This mapping was performed separately at pre- and post-sessions for motor impairment and from baseline to follow-up for treatment responsiveness (Fig.5(A.1-2)). To carry out this feature mapping, we conducted a grid-search over different kernel functions (K) (i.e. linear, Gaussian, and polynomial functions) and identified parsimonious clusterings.

We found that for distinguishing participants’ motor impairment level using redundant modular activations (R_Pre, R_Post), the linear kernel (i.e.) was most effective, while for synergistic activations (S_Pre, S_Post), the Gaussian kernel (i.e. exp , σ = 2) was best suited for distinguishing participants’ motor impairment. Finally, for treatment responsiveness (R_Pre-Post, S_Pre-Post), we found that the linear kernel provided a meaningful participant clustering. In the following, we briefly outline the steps involved in this clustering algorithm:

1. Kernel matrices were generated from each pairwise combination of extracted module activation coefficients across all participants (i.e.) (see Fig. 5(A.1) for clustering based on motor impairment and Fig. 5(A.2) for clustering based on treatment responsiveness). This results in a multiplex network (C) comprised of layers (c_i) for each module combination that described the similarity of each pair of participants’ modular recruitment (Fig. 5(A.1)). For impairment-based clustering (Fig. 5(A.1)), the diagonal of c_i (representing the self-similarity of each participant’s modular recruitment) was set to zero (Murphy et al., 2018), while for treatment responsiveness clustering the diagonal was maintained as it provided important self-referential information regarding the changes in participants’ modular recruitment patterns across the intervention (Fig. 5(A.2)).

2. The densely connected c_i (Fig. 5(B.1)) were then empirically sparsified using a modified percolation analysis before their modular structures were separately quantified using the Louvain algorithm (node colours in Fig. 5(B.2) represent different network communities) (Blondel et al., 2008; Gallos et al., 2012; Rubinov & Sporns, 2010). To offset any stochasticity in the solutions found in the final application of the Louvain algorithm, we re-applied it n = 5000 times and determined a consensus partition using an established method (Lancichinetti & Fortunato, 2012). Following this, co-membership matrices were yielded from all c_i partitions, where entries indicate whether a participant pairing belong to the same cluster or not (Fig. 5(B.3)). These co-membership matrices were then subsequently aggregated into a single representative matrix (Fig. 5(B.3)).

3. Finally, the Louvain algorithm is re-applied to this representative matrix to uncover patient clusters (Fig. 5(B.4)).

Matlab code for this clustering algorithm has been published here: https://github.com/DelisLab/NICF.

Quantifying the contributions of preservation, merging and fractionation to impairment-based patient clusters

To investigate whether the distinct patterns of preservation, merging and fractionation underlay the identified patient clusters in terms of motor impairment at pre- and post-sessions separately (i.e. [R_Pre, R_Post], [S_Pre, S_Post]), we implemented the following analysis.

To identify patterns of merging and fractionation among the functionally-similar and -complementary muscle networks of this post-stroke cohort, we aimed to identify linear combinations of representative affected- and unaffected-side muscle networks extracted across participants that could reconstruct individual participants’ muscle networks. More specifically, we concatenated all participants’ affected and unaffected-side redundant/synergistic muscle networks in a similar way as described in Equation 1.1 (see Extraction of Redundant and Synergistic Muscle Networks section of Materials and Methods) but separately for pre- and post-session data. We then determined the optimal number of components to extract for each session across participants using PNMF (see steps 2–3 of ‘Extraction of Redundant and Synergistic Muscle Networks’ from the Materials and Methods section).

To quantify merging events for each participant at baseline and follow-up sessions, we determined the weighted contributions of the extracted unaffected-side network components (v^u) to the reconstruction of the pth participant’s affected-side network using a non-negative least squares algorithm (‘lsqnonneg’ Matlab function) (Equation 2.1) (Cheung et al., 2012). Here, the non-negative, participant-specific merging coefficient (m) describes the contribution of the kth unaffected-side component from the set of extracted session-specific network components (N) to the reconstruction of each individual participant’s network:

In a similar vein, fractionation was quantified by essentially reversing this computation to determine participant-specific fractionation coefficients (f) for the kth affected-side network component extracted (v^𝒶) from individual participants’ unaffected-side muscle network (Equation 2.2):

To determine patterns of preservation for the pth participant at baseline and follow-up sessions, we quantified the inverse of the Euclidean distance as a measure of similarity (s) between the kth v^u and (Equation 2.3). Here, ranges from 0 to 1 with values close to 1 indicating perfect similarity between the networks:

Unlike previous research investigating merging and fractionation in stroke survivors (Cheung et al., 2012; Hashiguchi et al., 2016; Pregnolato et al., 2025), we did not threshold the resulting values or define discrete events. Instead, we took the raw and values as participant-specific measures of each physiological response at baseline and after treatment. Using these raw values, we applied a feature selection method (‘f_classif’ function in the sklearn Python package) to identify statistically significant (i.e. p < 0.05) independent contributors to the corresponding patient clusters related to motor impairment (i.e. R_Pre, S_Pre, R_Post, S_Post). Finally, to provide a parsimonious model of these significant contributors to each partition, we included them as explanatory variables in a binary logistic regression model with forward selection (Wald’s criterion: inclusion < 0.05, exclusion > 0.1) (SPSS Statistics 28).

Characterising rehabilitation effects on functional muscle interactions among responders and non-responders

To determine the effects of rehabilitation on functional muscle network structure among responders and non-responders, we firstly took an average interaction magnitude from the sparsified networks of participants affected and unaffected sides at each session. Using independent samples t-tests, we then determined statistical differences between patient groupings categorised based on the conventional MCID derived from the FMA-UE scale (i.e and ), and on the patient response clustering’s established using our framework (i.e. R_(Pre-Post), S_(Pre-Post)) (see ‘Clustering stroke survivors based on impairment severity and therapeutic responsiveness’ section of Materials and Methods). This enabled us to characterise the effects of rehabilitation on muscle network interdependencies and contrast between conventional and functional muscle network-based perspectives on therapeutic responsiveness.

Additionally, we aimed to quantify specific differences between responders and non-responders at each session at the level of individual muscle interactions. To do so, we employed a permutation-based approach to empirically derive interaction-specific significance thresholds and counteract family-wise error accumulation (Maris & Oostenveld, 2007). Briefly, for each grouping (i.e. clinical assessment-based and NIF) and muscle interaction across participants we randomly shuffled the group labels and determined statistical differences between the permuted groups using independent samples t-tests. This procedure was repeated n = 1000 times, generating a ground-truth distribution of test values from which the 95^th percentile value (i.e. p < 0.05) acted as the statistical significance threshold to compare against actual test values.

Statistical analyses

The association between the extracted redundant- and synergistic pre- and post-session activation coefficients and both pre- and post-session FMA-UE scores respectively were assessed using univariate linear regression analyses (“f_regression” function in the sklearn python package). In this way, we could identify modular recruitment patterns that were independently associated with functional impairment. Additionally, we sought to determine if changes to the recruitment of specific redundant and synergistic components reflected differences between participants based on treatment group (i.e. PT vs. VR) and treatment responsiveness (i.e. MCID on the FMA-UE (> 5 points = responder, ≤ 5 points = non-responder) (Page et al., 2012). To do so, we firstly computed the Canberra distance (i.e.) between corresponding pre- and post-session activation coefficients as a metric sensitive to small but potentially important differences in modular recruitment across sessions. Using these distance values, we then determined differences between treatment responders and non-responders using independent t-tests. To control for the imbalanced groups, differences between experimental groups were examined instead using the Welch’s t-test. Significance was set for all statistical analyses here a priori to p < 0.05.

Data Availability

The data is available upon reasonable request.

Acknowledgements

David O’Reilly and Ioannis Delis were funded by the Biotechnology and Biological Sciences Research Council (BB/Y513799/1). Giacomo Severini was funded by Research Ireland (Grant No. 21/FFP-P/10228). Giorgia Pregnolato was funded by the European Union’s Horizon 2020 research and innovation programme (Marie Skłodowska-Curie grant agreement No 101034252).

Supplementary Materials

Supplementary Information

The application of the NIF to the post-stroke cohort across both pre- and post-sessions identified five redundant (R1-R5 (Supplementary Materials Fig.1)) and seven synergistic (S1-S7 (Supplementary Materials Fig.2)) networks of co-occurring muscle interactions across affected and unaffected sides. We found a wide range of muscle interaction strengths (illustrated as the relative edgewidths on the human body models), principally contributing muscles (relative node sizes on human body model networks) and submodular affiliations (i.e. modules-within-modules) (node colours on the human body model networks) (see tables in Supplementary Materials Fig.1-2 for summary statistics including the range for the entire network (i.e. ‘Full Network’) and the ‘Extraction of Redun-dant and Synergistic Muscle Networks’ section of the Materials and Methods for further details). Interpreting these interaction strengths and principal muscles together with the submodular structure (see tables in Supplementary Materials Fig.1-2), we were able to attribute a biomechanical function to each extracted component.

Beginning with the redundant muscle networks (Supplementary Materials Fig.1), we found the biomechanical function of most components included the stabilisation of the shoulder girdle and/or elbow joint. These task-specific contributions to musculoskeletal stability were in some cases the predominant function of the muscle network (e.g. R3 (Affected-side), R5 (both sides)) and for others played a supportive role to elbow extension (e.g. R2 (both sides), R4 (Unaffected-side)) or forearm/shoulder rotation (e.g. R1 and R3 (both sides), R4 (Affected-side)). Examination of the sub-modular structure of the redundant muscle networks, which typically included 2-3 submodules, supported these functional interpretations. For example, for the affected-side R1 which related primarily to elbow extension and shoulder flexion (Fig.2), we found TM and TL were affiliated with separate submodules (green and red nodes respectively) including AD and LD (green nodes) and PD and BRACH (red nodes). This functional grouping reflects the distinct prime-mover and supportive contributions of the triceps brachii heads to elbow extension/shoulder flexion that are shared with specific elbow/shoulder muscles (i.e. TL contributes more to voluntary elbow extension than TM and is affiliated with shoulder/elbow stabilisers; TM contributes to shoulder stabilisation more so than TL and is coupled with shoulder flexors). Hence, knowing the activation state of TL or TM alone was sufficient to predict the motor task performed in an equivalent way as to observing the activation state of PD/BRACH or AD/LD respectively here. Differences were observed between the affected- and unaffected-sides of the redundant muscle networks indicative of motor impairment. R3 for instance (Supplementary Materials Fig.1), which involved interactions between SU and the shoulder musculature on both sides, displayed a specific reliance on PM on the affected-side whereas for the unaffected-side, the interactions were spread across PM, AD and LD, perhaps representing a loss of muscle selectivity at the shoulder level across the post-stroke cohort.

The extracted synergistic motor components (Supplementary Materials Fig.2) presented with widespread functional connections between proximal and distal musculature (e.g. S6 (Affectedside)). The extracted components were also less comparable across affected and unaffected sides compared to the redundant motor components (Supplementary Materials Fig.1), suggesting that the synergistic interactions more frequently reflected functional deficits induced by hemiparesis in the post-stroke cohort. Moreover, the submodular structure of the synergistic networks was more complex on the affected-side, ranging from 2-5 submodules compared to the 2-3 submodules typically found on the unaffected-side. The biomechanical function of several synergistic components was related to elbow flexion (i.e. S1, S3-S4 (both sides) and S2 (affected-side) (Supplementary Materials Fig.2)) in contrast to the elbow extensors common to the redundant components (Supplementary Materials Fig.1), while all extracted components on both sides were comprised of some form of task-specific joint stabilisation (i.e. shoulder girdle or elbow joint). This is intuitive as prime-mover and stabiliser muscles may also have synergistic roles (i.e. they predict complementary features of the motor task). The unaffected-side S7 exemplifies this intuition (Supplementary Materials Fig.2), comprising of the prime-mover AD during shoulder flexion coupled with several shoulder girdle stabilisers that offer complementary information regarding the specific orientation of the arm when executing the motor tasks. Similarly to the separate submodular affiliation of TM and TL in the affected-side R1 described previously (Supplementary Materials Fig.1), the affected-side S7 here comprised of separate affiliations for the biceps brachii heads (BS (orange nodes) and BL (black nodes) (Supplementary Materials Fig.2)). The complementary task information provided by BS here was more closely aligned with AD while BL was affiliated with BRACH, RM and Infra. These distinct functional affiliations were present despite a prominent functional interaction between BS and BL, suggesting that their complementary functions towards elbow flexion and shoulder girdle stabilisation to some extent overlapped.

Supplementary Figures

Figure S1.

Figure S2.

Figure S3.

Figure S4.