The application of the NIF across the post-stroke cohort and both pre- and post-sessions identified five redundant (R1–R5 (Figure 1—figure supplement 1)) and seven synergistic (S1–S7 (Figure 1—figure supplement 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 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 Figure 1—figure supplements 1 and 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 Figure 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 Figure 1—figure supplement 1); however, noticeable differences relating to the underlying impairment of the affected-side across the cohort were also represented (e.g. R3 (Figure 1—figure supplement 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 Figure 1—figure supplement 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.

In assessing the statistical relationships between the activation of specific redundant and synergistic network components and gross motor functional impairment, we found that the degree of activation of S2 and S6 (both comprised of complementary muscular contributions to elbow flexion with shoulder girdle stabilisation (Figure 2A) and Figure 1—figure supplement 2) at baseline was associated with lower FMA-UE scores both at baseline ($\beta$=–174.5 ± 62.5 ($p=0.008$), $\beta$ = –159.7 ± 45.6 ($p=0.001$) respectively) and at follow-up ($\beta$=–172.5 ± 65.4 ($p=0.012$), $\beta$ = –155.1 ± 48.2 ($p=0.003$), respectively; Figure 2A). This suggests that specific networks of synergistic muscle interactions are strongly linearly related to 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 to participant and modality provided.

Figure 2

Download asset Open asset

Examination of the extracted activation 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 (Figure 2B). Specifically, we found that among non-responders S3 activation changed significantly after treatment than responders (t=3.1 ($p=0.004$)), while the VR group demonstrated a significant change in the activation 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 (Figure 2B, Figure 1—figure supplement 2), suggesting that non-responders express maladaptive activation adjustments of specific flexor-driven networks, while VR treatment promotes changes in other flexor-driven networks.

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 rehabilitation 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 Figure 3A and B.1-2) and non-responders (i.e. <5 point change in FMA-UE scores Figure 3C and D.1-2).

Figure 3 with 1 supplement see all

Download asset Open asset

Treatment responders experienced a significant decrease in redundancy on the affected-side (t=2.3; p < 0.05; Figure 3A.1) and significant increase in synergy on the affected-side also (t=–2.2; p < 0.05; Figure 3B.1). No noticeable changes were experienced for average redundancy or synergy on the unaffected-side across responders were observed (Figure 3A-B.1). Specific muscle interactions also reflected this trend for redundancy and synergy with rehabilitation among responders (see network edges in Figure 3A-B.2). Whereas for the corresponding muscle interactions on the unaffected-side (Figure 3A.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 (Figure 3B.2), except for BRACH-SU which was lower at follow-up. This was also the case for the unaffected-side of responders (Figure 3B.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 (Figure 3C-D.1–2). Several redundant and synergistic muscle interactions on the unaffected-side were found to increase from pre-to-post sessions (Figure 3C-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 re-differentiation of muscle interactions provided by effective rehabilitation.

Implementation of our clustering approach (see ‘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 activation patterns) identified two redundant and synergistic patient subgroups ($R_{\text{Pre-Post}}$, $S_{\text{Pre-Post}}$). Figure 3—figure supplement 1 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_{\text{Pre-Post}}$ and $S_{\text{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_{\text{Pre-Post}}$ (Cluster 1 (Figure 3E.1–2), Cluster 2 (Figure 3G.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; Figure 3E.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 strength with rehabilitation (Figure 3E.2). Meanwhile, cluster 2 participants demonstrated a non-significant increase in the average redundancy of the affected-side and no perceived changes on the unaffected-side (Figure 3G.1). $S_{\text{Pre-Post}}$ (Cluster 1 (Figure 3F.1–2), Cluster 2 (Figure 3H.1–2) demonstrated a significant increase in synergistic muscle interaction strengths among cluster 2 participants (t=–3.78; p < 0.001; Figure 3H.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 (Figure 3H.2). On the other hand, for cluster 1 participants (Figure 3F), we found a non-significant decrease in synergy on the affected-side and no perceivable changes on the unaffected-side (Figure 3F.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 (Figure 3F.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.

Separate cluster analysis of pre- and post-session modular activation patterns across the post-stroke cohort (see ‘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_{\text{Pre}}$ and $S_{\text{Pre}}$ (Figure 4A–B)) and at follow-up (i.e. $R_{\text{Post}}$ and $S_{\text{Post}}$ (Figure 4—figure supplement 1A-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.

Figure 4 with 1 supplement see all

Download asset Open asset

A significant difference (*) between clusters was found for $R_{\text{Pre}}$ for both pre-session FMA-UE (t=2.27 ($p=0.029$) (Figure 4A.2)) and post-session FMA-UE (t=2.37 ($p=0.023$) (Figure 4A.3)), but not for the change in FMA-UE score from pre-to-post sessions (ΔFMA-UE) (Figure 4A.4). The $S_{\text{Pre}}$ partition also significantly differentiated participants based on motor impairment at baseline (t=–2.9, p < 0.01(**) (Figure 4B.2)) and follow-up sessions (t=–2.5, $p=0.017$(*) (Figure 4B.3)), but similarly to $R_{\text{Pre}}$, could not distinguish participants based on ΔFMA-UE (Figure 4B.4). For both $S_{\text{Post}}$ and $R_{\text{Post}}$, presented in Figure 4—figure supplement 1A-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 Figure 4A-B, Figure 4—figure supplement 1A-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_{\text{Pre}}$ and $S_{\text{Pre}}$ patient clusters are illustrated in Figure 4C-D while the corresponding outputs for $R_{\text{Post}}$ and $S_{\text{Post}}$ are presented in Figure 4—figure supplement 1C-D. For interpretation, alongside 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_{\text{Pre}}$ (Figure 4C), 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_{\text{Pre}}$ (Figure 4D.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 (Figure 4D.1)) was again an important contributor to the observed patient groupings, however in contrast to the redundant modules (Figure 4C), 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_{\text{Pre}}$, demonstrating how the maladaptive merging of elbow flexor related interactions (i.e. BL-BS (S1); Figure 4D.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_{\text{Post}}$ (β = -2.32±0.99 ($p=0.019$) (Figure 4—figure supplement 1C.1)). The underlying contributions to $R_{\text{Post}}$ clusters were more complex than $R_{\text{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$); Figure 4—figure supplement 1C.2). Together, these contributions illustrate that cluster 1 participants of $R_{\text{Post}}$ had significantly greater fractionation of R4 and preservation of R1 compared to cluster 2 participants at follow-up.

Finally, patterns of fractionation (S1 [Figure 4—figure supplement 1D.1]) and of merging (S2 [Figure 4—figure supplement 1D.2]) were found to significantly contribute to the $S_{\text{Post}}$ clusterings, both associated with cluster 2 affiliation ($\beta$ = 26.7±12.1 ($p=0.027$) and $\beta$ =27.3±11.8 ($p=0.021$), respectively [Classification accuracy: 76.2%]). S1 here, like S3 in $S_{\text{Pre}}$ (Figure 4D.1), comprised of PM playing a central role in interaction with all other muscles (Figure 4—figure supplement 1D.1). There wasn’t however a clear correspondence between the contributing merging patterns for $S_{\text{Pre}}$ and $S_{\text{Post}}$ (Figure 4D.2 and Figure 4—figure supplement 1D.2 respectively). Instead, for $S_{\text{Post}}$ the merging pattern significantly differentiating participants at post-treatment predominantly involved the SU.

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 and 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 Figure 1—figure supplements 1 and 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 networks closely reflected motor impairment consistently across the intervention and both treatment-type and -responsiveness (Figure 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 and O’Dwyer, 2018; Cheung and 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 Figure 1—figure supplements 1 and 2; Ó’ Reilly and Delis, 2022; O’Reilly and Delis, 2024). 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 Figure 1—figure supplements 1 and 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).

The fractionation of motor patterns into separate modules in response to neurological insult emerged as a central difference between the patient clusters identified here (see Figure 4, Figure 4—figure supplement 1). We found new evidence for the adaptive role of fractionation among redundant muscle interactions but maladaptive role of fractionation among synergistic muscle networks (Figure 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., 2026). Although fractionation was found to also differentiate patient clusters at follow-up (Figure 4—figure supplement 1C-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 (Figure 4—figure supplement 1A&B.3). Finally, the observations of both the merging and preservation of motor patterns contributing to gross motor function 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 training experience (Cheung et al., 2020; García-Cossio et al., 2014).

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 (Figure 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_{\text{Pre-Post}}$, $S_{\text{Pre-Post}}$) where it was in fact further accentuated. Previous research has shown that functional redundancy increases post-stroke (Cheung et al., 2012; Clark et al., 2010), reflecting the characteristic loss of functional specificity (i.e. functional de-differentiation) of muscle interactions post-stroke. Enhanced synergy with treatment here thus reflects the functional re-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 mechanisms 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 affecting 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 and Baker, 2022; Sangari and 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 to execute the same function (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 and Delis, 2024; O’Reilly and Delis, 2025a).

Through the novel clustering of stroke survivors presented here, we consistently identified dichotomous subgroupings that reflected distinct patterns of motor impairment and treatment responsiveness (Figures 3 and 4, Figure 3—figure supplement 1, Figure 4—figure supplement 1). These clusterings are formed through the greater affinity of 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 and Delis, 2025b; O’Reilly and Delis, 2024), 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 and 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 five further subclusters, aligning with previous research (O’Reilly and Delis, 2025b; Scano et al., 2017; van der Vliet et al., 2020). Indeed, in future work we aim to apply manifold learning approaches to the co-membership matrix derived from this clustering algorithm, providing a continuous representation of the population structure. Nevertheless, 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).

Although the FMA-UE is a gold standard measure of post-stroke treatment outcomes (Fugl-Meyer et al., 1975; Page et al., 2012), it does not capture the impact of rehabilitation on patients’ ability to perform activities of daily living or to participate in daily life. Hence, interpretations of the identified biomarkers are currently limited to gross motor function impairment and recovery. Future research should employ this framework to quantify biomarkers that correspond to other important aspects of patients’ recovery (e.g. functional independence, subjective experiences), thus offering a more complete evidence base for its clinical utility.

We conducted a secondary analysis on previously published data (Maistrello et al., 2021; Pregnolato et al., 2026; 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 (1 hr/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 PT (no. of participants enrolled = 8) or VR (no. of 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., 2026; Figure 1B)). The Fugl-Meyer Assessment Upper Extremity section (FMA-UE), a gold standard clinical measure for assessing gross motor functions was used to quantify motor impairment in the affected upper-limb (Fugl-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 recording 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; Figure 1C]) 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, bi-directional Butterworth filter with zero-phase distortion.

To quantify the task-relevant information shared by interacting muscles 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 and Delis, 2024). More specifically, leveraging a higher-order information-theoretic measure known as co-information ($II$) (Figure 1A, Ince et al., 2017, McGill, 1954), 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 ($\tau$), formally disentangling it from the task-irrelevant space (yellow intersection) (O’Reilly and Delis, 2024). $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 $\tau$) 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;\tau )+I(m_y;\tau )$, and shared together, $I(m_x,m_y;\tau )$ (Figure 5). 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 and O’Dwyer, 2018), formally distinguishing between two alternative ways in which muscles may ‘work together’ towards task performance.

Figure 5 with 1 supplement see all

Download asset Open asset

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 and Delis, 2022) by concatenating the 10 repetitions of each task executed on a particular side (i.e. variables of length no. of timepoints x 10 trials) and quantifying $II$ with respect to this discrete task parameter codified to describe the motor task performed at each timepoint for each trial included. This computation was performed on all unique $m_x$ and $m_y$ pairings, generating symmetric matrices ($A$) (i.e. $A=A^T$) composed separately of non-negative redundant and synergistic values (Figure 5). Having quantified these distinct types of task-relevant muscle networks for each participant on both affected and unaffected sides, the following steps in the NIF pipeline were then performed:

• Step 1: Network sparsification. We empirically sparsified each $A$ separately 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.

• Step 2: Model-rank specification. We then determined the functional organisation of each $A$ using a link-based community detection method (implemented with single linkage) (Ahn et al., 2010; Figure 2C.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 and Delis, 2024; O’Reilly et al., 2025c). A consensus partition across participants was then subsequently determined by aggregating across all co-membership matrices and then applying the Louvain algorithm to the resultant matrix (Blondel et al., 2008; Rubinov and Sporns, 2010).

• Step 3: Dimensionality reduction. The number of modules identified in this final consensus partition of step 2 was used as the input parameter for dimensionality reduction, namely projective non-negative matrix factorisation (PNMF; Figure 1D, Yang and Oja, 2010). The choice of PNMF here which has been employed in earlier work on the NIF (Ó’ Reilly and Delis, 2022), in contrast to the space-time tensor decomposition employed in the parent study (O’Reilly and Delis, 2024), was chosen simply to maintain brevity by focusing the analysis on the spatial domain. Continuing, the input matrix for PNMF consisted of the sparsified $A$ on both affected and unaffected sides from all participants at both pre- and post-sessions concatenated together in their vectorised form. More specifically, the input matrix composed of redundant or synergistic values was configured such that the set of unique muscle pairings ($1\dots K$) on affected and unaffected sides ($x_{\text{aff}}$ and $x_{\text{unaff}}$ respectively) corresponding to an individual participant ($p$) from the group of participants ($1\dots 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 ($v$) across affected and unaffected sides and their underlying session- and participant-specific activation coefficients ($a$).

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 and 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 and 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 andcopy archived at O’Reilly, 2026b. To facilitate comparisons with existing approaches, we performed a conventional muscle synergy analysis on the post-stroke cohort (see Figure 5—figure supplement 1 and associated text). Further comparisons with conventional approaches can be found in our previous work developing this framework (Ó’ Reilly and Delis, 2022; O’Reilly and Delis, 2024).

To understand how the stroke survivors compared in terms of gross motor impairment severity and treatment-responsiveness using NIF, we conducted cluster analyses using the extracted activation coefficients ($a$). We focussed on $a$ here as the extraction of population-level functional modules enabled the buffering of individual differences into the space of modular activations, making them an ideal target for identifying population structure. 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 (see Figure 5—figure supplement 1 and associated text for a direct comparison with conventional approaches), perhaps resulting from the complex interdependencies between the modular activations. Therefore, we sought to improve the linear separability of the participants by mapping $a$ 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 (Figure 6; O’Reilly and Delis, 2025b). This mapping was performed separately at pre- and post-sessions for motor impairment and from baseline to follow-up for treatment responsiveness (Figure 6A.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.

Figure 6

Download asset Open asset

We found that for distinguishing participants’ motor impairment level using redundant modular activations ($R_{\text{Pre}},R_{\text{Post}}$), the linear kernel (i.e. $a_i^{(p^{1\dots N})}\times a_j^{(p^{1\dots N})}$) was most effective, while for synergistic activations ($S_{\text{Pre}},S_{\text{Post}}$), the Gaussian kernel (i.e. $\exp (-\Vert a_i^{(p^{1\dots N})}-a_j^{(p^{1\dots N})}\Vert /(2{\sigma }^2))$, $\sigma =2$) was best suited for distinguishing participants’ motor impairment. Finally, for treatment responsiveness ($R_{\text{Pre-Post}},S_{\text{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 (Figure 6):

• Step 1: Multiplex feature mapping. Local kernel matrices were generated from each pairwise combination of extracted module activation coefficients across all participants (i.e. $K(a_i^{(p^{1\dots N})},a_j^{(p^{1\dots N})})$) (see Figure 6A.1 for clustering based on motor impairment and Figure 6A.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 activation (Figure 6A.1). For impairment-based clustering (Figure 6A.1), the diagonal of $c_i$ (representing the self-similarity of each participant’s modular activation) 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 activation patterns across the intervention (Figure 6A.2).

• Step 2: Sparsification and patient clustering. The densely connected $c_i$ (Figure 6B.1) were then empirically sparsified using a modified percolation analysis before their modular structures were separately quantified using the Louvain algorithm (node colours in Figure 6B.2 represent different network communities) (Blondel et al., 2008; Gallos et al., 2012; Rubinov and 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 and Fortunato, 2012). Following this, co-membership matrices were yielded from all $c_i$ partitions, where entries indicate whether a participant pairing belongs to the same cluster or not (Figure 6B.3). These co-membership matrices were then subsequently aggregated into a single representative matrix (Figure 6).

• Step 3: Patient cluster identification. Finally, the Louvain algorithm is re-applied to this representative matrix to uncover patient clusters (Figure 6B.4).

Matlab code for this clustering algorithm has been published here: https://github.com/DelisLab/NICF and copy archived at O’Reilly, 2026a.

To investigate whether the distinct patterns of preservation, merging and fractionation underlie the identified patient clusters in terms of motor impairment at pre- and post-sessions separately (i.e. $[R_{\text{Pre}},R_{\text{Post}}],[S_{\text{Pre}},S_{\text{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 (${\mathbf{v}}^u$) to the reconstruction of the Pth participant’s affected-side network (${\mathbf{w}}_p^a$) 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 (${\mathbf{v}}^a$) from individual participants’ unaffected-side muscle network (${\mathbf{w}}_p^u$) (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 ${\mathbf{v}}^u$ and ${\mathbf{w}}_p^a$ (Equation 2.3). Here, $s_p^k$ 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., 2026), we did not threshold the resulting values or define discrete events. Instead, we took the raw $m_k^p$, $f_k^p$ and $s_p^k$ 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_{\text{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_{\text{Pre}},S_{\text{Pre}},R_{\text{Post}},S_{\text{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).

To determine the effects of rehabilitation on functional muscle network structure among responders and non-responders, we firstly took the average functional 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. $R_{\text{Pre}}^{\text{FMA}},R_{\text{Post}}^{\text{FMA}}$ and $S_{\text{Pre}}^{\text{FMA}},S_{\text{Post}}^{\text{FMA}}$), and on the patient response clustering established using our framework (i.e. $R_{(\text{Pre-Post})},S_{(\text{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 and 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 95th percentile value (i.e. p < 0.05) acted as the statistical significance threshold to compare against actual test values.

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 activation patterns that were independently associated with functional impairment. Additionally, we sought to determine if changes to the activation 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. $\frac{|a_i^{\text{Pre}}-a_i^{\text{Post}}|}{\left|a_i^{\text{\% Pre}}\right|+\left|a_i^{\text{Post}}\right|}$) between corresponding pre- and post-session activation coefficients ($a_i^{\text{Pre}},a_i^{\text{Post}}$) as a metric sensitive to small but potentially important proportional differences in modular activation 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.

To more directly illustrate the advantages of the proposed framework, we carried out a standardised pre-processing of the EMG data in line with conventional muscle synergy analysis. This included rectification, low-pass filtration (cut-off: 20Hz) and smooth resampling of EMG waveforms to 50 timepoints. All data for each participant at each session was separately normalised by channel-wise variance, concatenated together and input into non-negative matrix factorisation (NMF; ‘nnmf’ Matlab function, 10 replications) to extract 11 muscle synergies ($W_{1-11}$ of Figure 5—figure supplement 1(Left)) and their time-varying activations. The number of components to extract was determined in a conventional way as the number of components required to explain >75% of the data variance. The extracted muscle synergies included distinct shoulder- (e.g. $W_2$), elbow (e.g. $W_8$), and forearm-level (e.g. $W_1$) muscle covariation patterns along with more isolated muscle contributions (e.g. $UT$ in $W_3$, $TL$ in $W_{10}$).

Regarding the clustering results of our framework and how they compare to conventional approaches, to facilitate this comparison we applied agglomerative clustering to the time-varying activation coefficients of all participants, trials, tasks separately for pre- and post-sessions and employed the ‘evalclusters’ Matlab function (ward linkage clustering, Calinski Harabasz criterion, Klist search = 2:21) for each session. We identified two clusters both at pre-session (criterion = 1.69) and post-session (criterion = 1.81) as optimal fits to the population data (see Figure 5—figure supplement 1 -2 [right]). We found no associations between pre- or post-session cluster partitions and participants FMA-UE scores. Nevertheless, we did identify significant associations between the pre-session clusterings and $S_{pre}$ (${\chi }^2$ = 7.08, $p$ = 0.008) and between post-session clusterings and conventionally-defined treatment responders (${\chi }^2$ = 4.2, $p$ = 0.04). These findings, along with the similar two-way clustering structure found using the NIF, highlight important commonalities between these approaches. To summarise the main advantages of our framework over this conventional approach:

• Improved detection of clinically relevant population structure. The clustering component of our framework revealed distinct post-stroke subgroups with important clinical relevance, distinguishing moderately and severely impaired cohorts and treatment responders and non-responders from pre-treatment data.

• Enhanced interpretability of extracted components and clusters. As our framework maps muscle interactions to a specific task parameter, we yield population-level motor components that correspond more consistently to meaningful biomechanical and physiological functions that can be interpreted across the dimensions of the specified task parameter. The proposed clustering approach also offers enhanced interpretability, addressing key limitations in the application of clustering approaches to the activation space of conventional muscle synergy analysis (e.g. different activation timings; Scano et al., 2017).

• Integration of pairwise muscle relationships. By incorporating muscle-pair level analysis, our framework captures coordinated interactions between primary and stabilising muscles—relationships that conventional NMF approaches overlook.

• Separation of task-relevant and task-irrelevant activity. The NIF isolates task-relevant coordination patterns, distinguishing them from task-irrelevant interactions driven by biomechanical or task constraints. On the other hand, task-relevant and -irrelevant muscle contributions are intermixed in conventional muscle synergy analysis.

• Ability to identify complementary functional roles. The NIF characterises whether muscle pairs act in similar or complementary ways, providing richer insight into motor control strategies.

• Reduced dependence on variance-based optimisation. Unlike conventional methods that rely on maximising variance explained, our framework allows detection of subtle but functionally significant interactions that contribute less to total variance.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

David O’Reilly and Ioannis Delis were funded by the Biotechnology and Biological Sciences Research Council (BB/Y513799/1). This publication emanated from research supported by the European Union’s Horizon 2020 research and innovation programme (Marie Skłodowska-Curie grant agreement No 101034252).

Informed consent and consent to publication was obtained from all participants in the study. Ethical approval was provided by the Ethical Committee of the IRCSS San Camillo Hospital, s.r.l., Venice, Italy.

1. Sent for peer review: August 17, 2025
2. Preprint posted: August 18, 2025
3. Reviewed Preprint version 1: October 24, 2025
4. Reviewed Preprint version 2: February 18, 2026
5. Reviewed Preprint version 3: June 18, 2026
6. Version of Record published: July 23, 2026

You can cite all versions using the DOI https://doi.org/10.7554/eLife.108509. This DOI represents all versions, and will always resolve to the latest one.

© 2025, O'Reilly, Pregnolato et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.