Recording activity in proximal muscle networks with surface EMG in assessing infant motor development

OBJECTIVE
To develop methods for recording and analysing infant's proximal muscle activations.


METHODS
Surface electromyography (sEMG) of truncal muscles was recorded in three months old infants (N = 18) during spontaneous movement and controlled postural changes. The infants were also divided into two groups according to motor performance. We developed an efficient method for removing dynamic cardiac artefacts to allow i) accurate estimation of individual muscle activations, as well as ii) quantitative characterization of muscle networks.


RESULTS
The automated removal of cardiac artefacts allowed quantitation of truncal muscle activity, which showed predictable effects during postural changes, and there were differences between high and low performing infants.The muscle networks showed consistent change in network density during spontaneous movements between supine and prone position. Moreover, activity correlations in individual pairs of back muscles linked to infant́s motor performance.


CONCLUSIONS
The hereby developed sEMG analysis methodology is feasible and may disclose differences between high and low performing infants. Analysis of the muscle networks may provide novel insight to central control of motility.


SIGNIFICANCE
Quantitative analysis of infant's muscle activity and muscle networks holds promise for an objective neurodevelopmental assessment of motor system.


Introduction
With the survival rate of preterm or low birth weight infants improving during the past few decades, there has been an increasing number of children at a neurodevelopmental risk. Currently more than every tenth infant is considered to be at a risk of developing motor, sensory, cognitive, and many other neurodevelopmental problems later in life. Motor problems can range from developmental coordination disorders to cerebral palsy (Hirvonen et al., 2014). Consequently, it would be important to identify specifically the infants with motor dysfunction who would benefit from an early intervention (Novak et al., 2017). Abbreviations: CC, correlation coefficient; CC P95 , the 95th percentile of surrogate correlation coefficient values; ECG, electrocardiogram; EEG, electroencephalography; EMG, electromyography; E q., Equation; GM, general movements; ICA, independent component analysis; LAbd, Left abdominal muscle (musculus rectus abdominis); LDel, Left deltoid muscle; LoBack, lower back muscles (musculus erector spinae at L4-L5 paraspinal level); LPect, Left pectoral muscle; MAI, muscle activation index; Neck, neck extensor muscles (musculus semispinalis capitis); Nedges, the number of significant edges; RAbd, Right abdominal muscle (musculus rectus abdominis); RDel, right deltoid muscle; RPect, Right pectoral muscle; SCC, surrogate correlation coefficient; SD, standard deviation; sEMG, surface electromyography; SST, surrogate signal test; UpBack, upper back muscles (musculus erector spinae at Th3-4 paraspinal level); Z-sum, Z-score sum.
Despite the significant progress in recording technology, there is a notable mismatch between current methodology and the clinical perception of the problems in infant neuromotor development. While the current literature is mainly focused on measuring muscle activation during specific tasks (van Balen et al., 2012a;Boxum et al., 2014;de Graaf-Peters et al., 2007;Hadders-Algra et al., 1992;Hamer et al., 2016;Van Der Heide et al., 2004;Hirschfeld, 1997;Ritterband-Rosenbaum et al., 2017a;Sansom et al., 2013;Sundermier et al., 2001;Washington et al., 2004;Xiong et al., 2018;Zhvansky et al., 2015), the clinical assessment typically observes coordination of muscle ensembles during spontaneous behaviour (Hadders-Algra, 2018a). To the best of our knowledge, only few studies have assessed the spontaneous motor activity of infants as part of basic research projects of motor development (Ritterband-Rosenbaum et al., 2017a;Hadders-Algra et al., 1992;Hadders-Algra et al., 1997;Sylos-Labini et al., 2020) but no clinical applications were derived. In addition, as the primary postural control during the first few months of infant development aims at the maintenance of head and trunk posture against the forces of gravity (Hadders-Algra, 2018a; van der Heide and Hadders-Algra, 2005; Pavão et al., 2013), efficient monitoring of truncal muscles in infants seems appealing.
The mismatch between the current sEMG literature and clinical practice (Campanini et al., 2020) may result from several technical challenges related to signal analysis: First, sEMG recordings from truncal and other proximal muscles include very large cardiac artefacts, which are challenging to remove, especially when the cardiac artefact changes during infant's movement activity (Willigenburg et al., 2012). Recent methodological developments suggest that this could be a solvable problem, though it requires customizing for infant recordings. Second, analysis of spontaneous movements in muscle ensembles cannot be approached by quantifying taskbased muscle energies, but it requires the design of novel analytic approaches to account for uncontrolled and variable motility patterns. Third, all the wireless electrode solutions are currently too extensive to implement and might significantly hinder infants spontaneous and voluntary movements. The endorsed normalization methods, e.g. the ones based on maximal voluntary effort (MVC), cannot be implemented among infant population (Burden, 2010). Therefore, most of the infant surface EMG studies have focused on analysing quantitative EMG-parameters (e.g. duration and co-activation pattern) during very specific and precisely defined voluntary activities or perturbations (van Balen et al., 2012a;Boxum et al., 2014;Chen et al., 2018;de Graaf-Peters et al., 2007;Hamer et al., 2016;Hirschfeld, 1997;Sundermier et al., 2001;Washington et al., 2004;Van Der Heide et al., 2004;Sansom et al., 2013;Xiong et al., 2018;Zhvansky et al., 2015).
A key characteristic of properly operating motor system is the body-wide integration of muscle control by central nervous system. This muscle coordination gives rise to apparent ''muscle net-works" (Kerkman et al., 2018) where activation in groups of muscles show significant temporal correlation. Such muscle behaviour renders statistical analysis of muscle coordination comparable to approaches used recently in studies of electroencephalography (EEG)-based functional brain networks (Bassett and Bullmore, 2009). Analysis of muscle networks has shown promise in adult studies evaluating neural control of motor coordination (Boonstra et al., 2015;Kerkman et al., 2018Kerkman et al., , 2017Naro et al., 2019).
This study was set out to develop methods for studying proximal muscle activity with surface EMG in infants during motor tasks and spontaneous motor activity. To this end, we developed methods for removing cardiac artefacts in recordings during movement activity, and we developed statistical methods to assess infant's muscle networks. Finally, we tested their technical suitability for measuring infant motility in early developmental screening.

Materials and methods
The overview of the data acquisition process and analyses of this study is represented in Fig. 1. After data collection the overall workflow consisted of (i) the design of methods to analyse infant muscle activation efficiently and reliably, (ii) the development of techniques for muscle network analysis in infants.

Participants
Eighteen infants born at term age were met at 3 months of age (mean age 3.3 months, 0.357 (SD), range 2.8-4.0; 40 % male) without any prior concerns of development or history of significant medical issues. The subject characteristics are presented in detail in the Supplementary Material (Table S1). This cohort was recruited from a larger ongoing research project. The age period between 2 and 4 months was chosen since it is considered to represent a clinically interesting transition in motor development. At this age period, infants still express endogenously generated spontaneous movements, which are gradually changing to voluntary movement repertoire with environmental adaptation (Hadders-Algra, 2018a). The study complies with the Declaration of Helsinki and the research was approved by the Ethics Committee of Children's Hospital, Helsinki University Hospital. A written informed consent was obtained from the parents.

Surface EMG recordings
The recordings were done using wireless amplifier Bittium BrainStatus TM (model 9305030A01, Bittium Biosignals Ltd, Oulu, Finland). The amplifier was chosen due to its small size, easy mobility and good signal quality. This amplifier design was also easy to position close proximity of the infant, usually near the head with electromyography (EMG) leads bundled together out of the body overall at the neck opening. The technical specifications for this amplifier were as follows: analog-to-digital converter precision 24 bits, input range ± 375 mV, sampling frequency 250 Hz, and common mode rejection > 100 dB. No online filtering was used. Prior studies have shown that EMG sampling rate of 250 Hz may be suboptimal for the full spectral content of the sEMG signal, however such undersampling does not significantly affect estimates of average sEMG amplitude or total sEMG area (Ives and Wigglesworth, 2003).
EMG was recorded by regular surface electrodes (AmbuBlueSensor NF-50-K). The skin areas were swiped with alcohol-containing solution (Septidin Solution, Takeda Ltd, www.takeda.com) the sensors were affixed to the skin and the electrodes were further secured by using hypoallergenic tape. We used inter-electrode distance of 30 mm for the dorsal extensor muscles and 15 mm for the other muscles. Left and right dorsal muscles were registered with the same bipolar signal pair, because we were only interested in differentiating the rostro-caudal differences in the dorsal extensor muscles. The recording reference was placed over sternum.
We also manufactured a full body overall with integrated textile electrodes to see whether single electrode placement could be replaced with a wearable solution for an improved wearing comfort and easing the setting up the system, when the suit provides the fixed location -platform for electrodes (Airaksinen et al., 2020;Colyer and McGuigan, 2018;Michelsen et al., 2020). Due to issues with electrode-skin interface, however, this solution did not provide high enough signal quality for further signal analyses (see the Supplementary Material for details).

The choice of the recording constellation
The major focus of this study was on the truncal musculature because of the clinical and research-based notion that problems of adaptive postural muscle control are one of the earliest clinical signs of impaired motor development 3 . We were also motivated by the concept that truncal muscles would theoretically be easier to incorporate into possible wearable solutions in the future. Therefore, the following muscles were selected (Fig. 1A): neck extensors (over m. semispinalis capitis, ''Neck"), thoracic extensors (m. erector spinae at Th3-4 paraspinal level, ''UpBack"), lumbar extensors (m. erector spinae at L4-L5 paraspinal level, ''LoBack"), right and left m. deltoideus (''RDel" and ''LDel"), right and left m. rectus abdominis (''RAbd" and ''LAbd"), right and left m. pectoralis major (''RPect" and ''LPect"). SENIAM guidelines (Hermens et al., 1999) were used for electrode locations when available, or the maximal muscle body was identified by palpation.
Due to unforeseen technical issues, one recording did not include neck muscles (subject 17), one recording did not have pectoral muscles (subject 18), and some infants had dorsal muscles recorded monopolarly instead of the intended bipolar derivation. These were taken into account in cohort analyses, and our further testing indicated no significant effects by these deviations (see the Supplementary Material for further information).

The course of the recording session
The recordings were carried out in a clinic like setting in 30-60 minutes long sessions that consisted of two phases, specific events, and spontaneous movements. The specific events included lying in supine, prone, pre-reaching/reaching, traction, suspension, and lateral tilting, where a pediatric physiotherapist assisted in the given postures and actions. A detailed verbal explanation of each event type is offered in the Supplementary Material. These events were chosen because they are typical components of standardized neurological assessments (Haataja et al., 1999;Romeo et al., 2016) and thus they would provide systematic and intuitive feedback for the analytic development. Lateral tilting was performed 1-3 times per side depending on infants tolerance. Spontaneous movement was assessed during prone and supine positions that were recorded for up to 3-5 minutes both, or as long as tolerated by the infant. In addition, an experienced pediatric physiotherapist assessed each participant with a structured neurological examination (Romeo et al., 2016).
All sessions were video recorded to allow an accurate offline epoch segmentation, as well as to evaluate the general movements (GM) as per Prechtl's method (Einspieler et al., 2004). The GMs and overall motor performance of each participant were later graded from the video recordings by an experienced child neurologist and a pediatric physiotherapist, allowing their posthoc grouping Overview of the data acquisition and analysis. a) Surface EMG (sEMG) was recorded from 18 infants at the age of 3 months. Nine truncal muscles were monitored during predetermined events of interest (prone, supine, tilting left, tilting right). b) Electrocardiogram (ECG) artifact was removed from EMG signals using dynamic template matching approach. Raw EMG (grey) was filtered within an infant-typical frequency band 15-70 Hz (black). A QRS-complex template (average across whole recording) was fitted for each cardiac cycle (red), and the ECG artifact was removed by subtracting the template from the pre-filtered EMG. Finally, for the cleaned EMG (black) amplitude envelope (green) was computed and smoothed (bold black) for further analyses. c) Magnitude of muscle activation during specific events was assessed using muscle activation index (MAI), which was compared to the inactive baseline. Functional interactions between muscles, or muscle networks, were assessed by computing correlations between (smoothed) amplitude envelopes, yielding connectivity matrices for each infant and event. The network correlations were assessed for statistical significance by using surrogates, created by time shuffling. The muscle activation index analyses were visualized by showing both the baseline activation (black circles) and the event activation (filled green circles), where the circle diameter is indicating the MAI magnitude. The muscle networks were visualized topographically as shown in the figure, where the thickness of the links/edges (red) reflects the correlation strength of significant interactions. into high vs low-performing groups (N = 8 and N = 7, respectively). These different motor performance groups highlight the wide spectrum of normal developmental variation. Those infants who were classified in high-performing group had already reached mature control and coordination in supine/prone and tilting postures. Those who were in low-performing group had still immature features in control or coordination of supine/prone or tilting postures (Haataja et al., 1999). Group details are provided in the Supplementary Material Table S1.

Data analysis
After visual review, annotation and epoch selection, the EMG data was analysed computationally, including (i) band-pass filtering for EMG frequency range, ii) removal of ECG artefacts, iii) assessment of muscle activation indices, and iv) assessment of muscle networks (Fig. 1b-d).
The visual EMG review was done using an open-source viewer EDFBrowser (www.teuniz.net/edfbrowser/) while all the subsequent computational analysis was performed using our custom scripted algorithm in MATLAB environment (MATLAB Ò R2017b, The MathWorks, Inc., Natick, Massachusetts, United States).
The analysis was carried out at three levels. First, we estimated the level of muscle activation during specific tasks relative to baseline, in order to assess the feasibility of the analytic pipeline in fully controlled postural changes. Second, co-activation of intuitively reasoned pairs of muscles were examined by estimating correlations between instantaneous muscle activities, i.e. EMG envelopes. Third, pairwise correlations in the global muscle network were estimated by pairwise comparison of instantaneous activity in all muscle pairs. In addition, the effect of ECG removal was carefully estimated to ascertain that, for instance, the tiny ECG residuals do not give rise to spurious muscle network measures (see the Supplementary Material for details).
3.1. Pre-processing of EMG data

Visual review
Multi-channel EMG signals were reviewed offline with the open-source viewer EDFBrowser to identify visually obvious technical and other artefacts for rejection from further analysis. Each video recording was independently annotated for infants posture and movement by one of the authors to allow independent epoch selection.

Removal of cardiac artefact and filtering
The electrocardiographic (ECG) artefacts in the EMG signals from truncal muscles create a significant challenge to all EMG analyses, which can only be solved by their adequate removal. However, their removal is also challenging due to their complicated and time-varying waveforms. After failing to remove the ECG artefact with the available methods, we designed, constructed, and implemented our own removal tool based on adaptive template matching (Xu et al., 2020). Illustration of the principle is shown in Fig. 1B, and more technical details are explained in the Supplementary Material. After ECG removal, the signals were filtered within 15-70 Hz frequency band and power line noise 50 Hz was removed with a notch filter.

Estimation of signal envelopes
Both muscle activation index (MAI) and muscle networks used amplitude time courses (envelopes) of EMG signals. In turn, envelopes were computed by taking absolute values of complex analytic signal (obtained by using Hilbert transform), and then smoothed with median filter (time window 0.4 s).

Muscle activation index (MAI)
Normalized muscle activation data was used to establish posture discrimination between different activity states: i) supine versus prone, and ii) tilting left versus tilting right. To achieve a more comparable analysis between the different muscles and subjects and to minimize the effect of noise, a specific EMG baseline was introduced for each infant with a total duration of approximately 10 sec. This baseline was visually determined by choosing representative periods of muscle inactivity based on both EMG signal and synchronized video. Muscle inactivity (''relaxed -looking state") was chosen for normalization because the standard normalization with maximal activity is not feasible in the non-cooperative infants. The muscle activation during the different activity states (events) of interest was then compared to the muscle activation during the inactive baseline. To achieve this, we established a muscle activation index (MAI) for each event and muscle separately.
Using baseline epochs from each separate EMG-channel, we computed normalization coefficients (C bline ) as the 75th percentiles in the amplitude distributions of their smoothed envelopes. Next, the whole EMG envelopes for each muscle signal were normalized as following: where EMG m sm t ð Þ and EMG m norm t ð Þ are initial smoothed and normalized envelope signals respectively, C m bline denotes the normalization coefficient for a considered muscle, m 2 M, M is the number of analysed muscle. Eventually, muscle activation index (MAI) for a given event and muscle was calculated as the mean activation as follows:MAI m event ¼ 1 where Dt refers to the length of time that child spends in the given position.
We assessed the magnitude of changes between each musclé s MAI during event relative to its baseline MAI, using the following ratio: The values from this equation were further depicted in the muscle activation visualization represented in the Fig. 1d (the radius of the green circles corresponds to the R m event value). In addition, to estimate the overall signal quality, an average amplitude for EMG envelope was defined for every channel both during the baseline and during the different events.

Muscle network analysis
We employed amplitude envelope correlations for infant muscle network analysis. Only prone and supine events were accepted for the network analysis because of the longest continuous duration of these events. In addition, only the subjects with over 30 seconds of continuous data for each event were included. The exact event durations for each subject can be found in Table S2 in the Supplementary Material. Eventually, muscle networks were reconstructed from EMG signals acquired during these two different states separately.The synchrony between muscle activations was estimated by computing Pearson correlation between corresponding (smoothed) EMG amplitude envelopes. Pearson correlation coefficient was chosen as the coupling metric due to its wide use in comparable studies (Omidvarnia et al., 2014). We therefore compared the correlation analysis estimates between Pearson and Spearman methods (see Supplementary Figure S1) which proved highly comparable results in datasets of this kind. From the point of view of network theory, individual muscles are network nodes and the strength of pairwise interactions between them (measured with Pearson correlation or corresponding zscores) are functional edges.

Statistical testing with surrogates
To assess the significance of the amplitude envelope correlations for a given muscle pair, we used surrogates. To test each connection for significance, we generated N = 100 surrogate correlation coefficients. To do so, we split one of the envelope signals into 2 second blocks and randomly shuffled them. At every iteration correlation between shuffled envelope and another original envelope was computed. Further, we computed the 95th percentile (CC p95 ) of surrogate correlation coefficient values. The correlations between original signals were considered significant if their correlation was greater than the CC p95 . Further rationale of surrogate testing is described in the Supplementary Material.
The network correlation coefficients (or network edges) were rendered comparable across infants by normalizing them via Zscores, which were obtained using individual's surrogates as: Where R is the correlation coefficient between a pair of muscles l surr ; and r surr are correspondingly mean and standard deviation of the surrogate correlation coefficients distribution. These Z-scores allowed computing more extensive intra-and interindividual statistics.

Statistical analyses
The time series data were analysed with MATLAB R2017b, while some statistical analyses were carried out in IBM SPSS Statistics version 25.0. The threshold of statistical significance was taken as p < 0.05, while significance level was adjusted with Bonferroni correction when multiple comparisons were performed. To assess the relationship between amplitude envelopes of different muscles, Pearsoń s correlation was calculated. The Mann-Whitney U test was used to assess differences between subject groups, while pairwise comparisons of different activations at a subject level were executed by using Wilcoxon signed-rank test.

Results
A total of 15 measurements were eventually accepted for further analysis. Three subjects were rejected because of either profuse amounts of signal artefacts (subjects 7 and 16) or poor cooperation resulting in discontinuation of the recording session (subject 13). One infant (subject 15) was later found to have significant neurological findings in the clinical assessment, and hence was excluded from the group statistics, unless otherwise mentioned. The average acceptance percent of the data was 81.0 % (range 74.5 -87.7 %, SD +/-16.2 %). See Table S2 in the Supplementary Material for additional details on data quality.
The absolute amplitudes during both baseline and event were highest in the dorsal muscles (see also the Supplementary Material), which might be due to simultaneous recording of a combination of bilateral extensor muscles and a longer interelectrode distance. Average amplitude during events was at lowest in the abdominal and deltoid channels, which is likely due to the abundant subcutaneous tissue in these regions.

Muscle activation during postural changes
Muscle activation index (MAI) was used to test how two paired activities can be distinguished: i) supine versus prone and ii) tilting left versus tilting right. These activations were chosen to assess the reliability and repeatability of our algorithm before progressing to network analysis. We expected to see a clear and intuitively interpretable change in MAI values between these event pairs in this sample of healthy infants.

Supine versus prone
A total of 13 subjects were included in the analysis (subject 4 was excluded due to poor cooperation). As expected, we found a significantly increased dorsal extensor muscle tone in prone position (neck z = 3.18, p < 0.01; upper back z = 3.30, p < 0.01 and lower back z = 3.30, p < 0.01). During change from supine to prone position, the relative average increase in MAI ((R Prone -R Supine )/R Supine -Â 100 %) was for neck muscles 289.5 % (SD +/-244), upper back muscles 247.8 % (SD +/-166) and lower back muscles 243.1 % (SD +/-142). In addition, there was a significant increase in median MAI in prone in most other muscles as well (see Fig. 2).
This finding appeared to be universal, as there were no significant differences (p > 0.05, Mann-Whitney U) in the median MAI levels between age groups (<3.4 months, N = 7 and ! 3.4 months, N = 6) or between infants with different motor performance (high, N = 7 vs. low performing, N = 6).

Tilting left versus tilting right
A total of 13 normal subjects were included in the analysis (subject 4 was excluded due to insufficient data). We found an asymmetric muscle activation during tilts in every subject (see Fig. 2b). The MAI in right abdominal muscle was significantly lower during tilting right than tilting left, z = -3.30, p = 0.01. Left abdominal muscle had significantly higher median MAI during tilting right than tilting left, z = 3.18, p = 0.01. The average relative change in MAI between tilting left and right was 86.0 % (SD +/-58) increase for right abdomen and 31.6 % (SD +/-19) decrease for left abdomen. In addition, right deltoideus also had a statistically significant higher median MAI when tilting right (z = -2.14, p = 0.01) with the average of 9.9 % (SD +/-18). There were no statistically significant differences in the activation energy of any other muscles between tilting left and right (p > 0.05).

Correlations in muscle pair activation during postural changes
We computed correlations between pairs of muscles that could be assumed to show synchronized activation patterns during prone or supine positions. Summary of muscle pair correlations are shown in Fig. 3. Comparison to surrogate data showed that all r > 0.15 values were statistically significant (p < 0.05) in all participants. See the Supplementary Material Table S3 and S4 for detailed results on subgroup analyses.

Group differences in prone position
The correlation of UpBack -LoBack activity was statistically significantly higher in the low performing group (median = 0.46) in comparison to the high performing group (median = 0.31) (U = 25, z = 2.08, p = 0.04, N = 11, Mann-Whitney U). Likewise, the correlation of UpBack -Neck activity was significantly higher in low performing group (median = 0.29) when compared to the high performing group (median = 0.07) (U = 26, z = 2.27, p = 0.02) (Fig. 3). No statistically significant difference was found in the correlations between the different age groups). These did not survive correction for multiple comparisons (three comparisons; corrected p threshold 0.017); however, the study was exploratory, and the findings are compatible with intuitive reasoning.

Group differences in supine position
The correlations between pectoralis muscle activations were significantly higher in the low performing group (median = 0.53) as compared to the high performing group (median = 0.38) (U = 35, x = 2.72, p = 0.01, Mann-Whitney U, N = 12). There were no significant differences between these groups in other muscle pairs, or between the age groups.

Outliers
There were two outliers (subjects 9 and 15) during supine and one outlier (subject 14) during prone position, as assessed as being greater than 3 box-lengths from the edge of the box in a box plot. Fig. 2. Results from muscle activation and average amplitude analysis. We estimated the muscle activation during specific events relative to baseline and motor performance. Muscle activation index (MAI) was used to establish posture discrimination between the following event pairs: i) supine versus prone (a) and ii) tilting left versus tilting right (b). a) All the dorsal extensor muscles showed a significantly higher increase in the median MAI during prone position in comparison to supine position when assessed by Wilcoxon signed-rank test. In the upper plot we show the muscle activation pattern of one of the test subjects (subject 5) during the different positions. In the lower plot we show the median muscle activation index values for each subject during supine (the small red dots) vs. prone (the small blue dots). The group median values are portrayed as filled circles and interquartile ranges (from 25th percentile to 75th) as lines. The significance of the difference in median value is visualized as a bigger filled bigger red circle over the muscle in question and the significance was further quantified as follows: * for p < 0.05, ** for p < 0.01, and *** for p < 0.001. b) In these plots the corresponding results for lateral tilting are portrayed. Left abdominal muscle had significantly higher median MAI during tilting right than tilting left as assessed by Wilcoxon signed-rank test. The opposite was true for the right abdominal muscle. In addition, the right deltoid muscle also had a statistically significantly higher MAI during tilting right. Results for the left deltoid muscle or other muscles were not statistically significant. In the upper plot the median values for one individual (subject 5) are shown as an example and in the lower plot all the individuals are mapped separately. c) Average amplitude of the EMG envelope was constructed over all the events (supine, prone and lateral tilting) combined and during baseline for each muscle. The average amplitude for all the muscles was 5.5 mv (SD 1.0) during baseline and 19.5 mV (SD 6.8) during the combined events. Muscle name abbreviations: Neck = neck extensors over musculus semispinalis capitis; UpBack = upper back muscles over m. erector spinae at Th3-4 paraspinal level; LoBack = lower back muscles over musculus erector spinae at L4-L5 paraspinal level; RDel = right deltoid muscle; LDel = Left deltoid muscle; RAbd = Right abdominal muscle (musculus rectus abdominis); LAbd = Left abdominal muscle (musculus rectus abdominis); RPect = Right pectoral muscle; LPect = Left pectoral muscle. Subjects 9 and 15 received the lowest motor grading of all participants. Subject 15 was hypotonic and had abnormal GMs, which justified exclusion from group statistics. Other outliers were included in the analysis since they did not have significant neurological findings. During supine position, subject 9 showed slightly atypical motor patterns and GMs. Subject 14 had external noise in the pectoral channels and an asymmetric prone posture observed offline in the video record.

Amplitude envelope correlations in global networks
The main results from the global muscle network analysis are represented in the Fig. 4. Visualization of individual muscle networks showed wide network patterns with significant correlations between most muscles. However, clear differences were also found between body postures (Fig. 4a).
The significant correlations were found between almost all muscle pairs in almost all infants (Fig. 4b), hence the main difference between infants or postures was to be found in the correlation strength, which could be summarized as the number of significant edges (or functional pairwise interactions measured with Pearson coefficient; Nedges) or taking into account the normalized strengths as Z-score sum (Z-sum) over all edges. As shown in Fig. 4c, the absolute correlations did not show posture-difference, but Z-score sum of edges was statistically significantly higher during supine position compared to prone position (z = -2.073, p = 0.038, Wilcoxon sign-rank test). No significant difference was observed in network summary metrics between groups with different motor performance (p > 0.05, Mann-Whitney U test).
While supine and prone were found to differ with respect to Zsum at the global level, there was only limited difference in the consistent edges. As shown in Fig. 4d, only 22.2% of all edges were Fig. 3. Amplitude envelope correlations between intuitive muscle pairs. Co-activation of intuitively reasoned muscle pairs were examined by estimating their mutual amplitude envelope correlation. Box plots show the median, lower, and upper quartiles (boxes) and the whiskers show the interquartile ranges (from 25th percentile to 75th). Mann-Whitney U test was used to assess the difference between the motor performance groups and p < 0.05 was considered statistically significant. The muscle pairs with statistically significant difference in the median amplitude envelope correlation values between the motor performance groups are marked with asterisk (*). a) In this plot we show the pairwise correlation coefficient values during supine position plotted separately for the different motor performance groups. The correlations between right and left pectoralis muscle activations were significantly higher in the low performing group (N = 7 for the high performing and N = 6 for the low performing group). b) The pairwise muscle correlation coefficient values during prone position in the different motor performance groups. The correlations between upper back (UpBack) -lower back (LoBack) and UpBack -Neck activities were statistically significantly higher in the low performing group (N = 7 for high performing and N = 4 for the low performing group). significantly different between the postures (Wilcoxon sign-rank test, N = 9).
Finally, we wanted to assess possible relationships between event durations and network measures when the minimum epoch length was set to 30sec. To achieve this a Spearmań s rank order correlation was run. Neither Z-sum nor Nedges were found to correlate with the duration of prone recordings (Nedges: r s (9) = 0.29, p = 0.39 ; Z-sum r s (9) = 0.48p = 0.13). In the supine posture, however, both Z-sum nor Nedges were increased with epoch duration (Nedges: r s (10) = 0.60, p = 0.04; Z-sum r s (10) = 0.75p = 0.01). The average supine event duration was 136.1 (SD +/-46.2) seconds in the high performing group and 172.9 (SD +/-57.5) seconds in the low performing skills group. The supine event duration for the subject with neurological findings and abnormal GMs was 143 seconds.
The average Z-sum among all normal infants with 9 monitored muscles was 182.1 (SD +/-52.0). The only infant with neurological findings and abnormal GMs (subject 15), also had the lowest Zsum (92.0) during supine position and the Nedges was also the lowest of all participants (24 vs average 30.7 (SD +/-4.8). This infant did not tolerate prone position long enough to be analysed. There were no significant differences in the medians of these network parameters between the different age or motor skills groups (Mann-Whitney U test, p > 0.05).

Discussion
Our study shows that a quantitative assessment of infant's event-related and spontaneous muscle activity is possible using multichannel surface EMG recordings of truncal and proximal muscles. Our work developed several novel analysis algorithms to first remove the large amplitude cardiac artefacts, and then to allow statistically rigorous assessment of spontaneous muscle activations, as well as quantitation of muscle networks. We validate Fig. 4. Muscle networks analysis. a) Examples of subject wise muscle networks during supine (left) and prone (right) position. The edges are shown after thresholding with surrogate testing, and the correlation strength is shown with the edge width.b) This network shows how many (N) infants out of all N = 13 infants did not show a significant correlation in the given edge in supine position. c) The muscle network summary metrics in the different positions and according to the motor performance groups (N = 13 for supine and N = 12 for prone). Significant differences were observed only in the Z-score sum between supine and prone positions (marked with asterisk). Box plots show the median, lower, and upper quartiles (boxes) and the whiskers show the interquartile ranges (from 25th percentile to 75th). The median Z-score sum was statistically significantly higher during supine position in comparison to prone position (p = 0.038, Wilcoxon sign-rank test). d) The difference between postures (supine vs prone) in the Z-score values at individual network edges, i.e., muscle pairs. The matrix shows p-values with different shades (Wilcoxon sign rank test, N = 9), and a thicker red line for those 22.2 % of all edges that were deemed to be significantly different between supine and prone posture. However, the emergent muscle pairs did not form any seemingly clear anatomical distribution. Muscle name abbreviations: Neck = neck extensors over musculus semispinalis capitis; UpBack = upper back muscles over m. erector spinae at Th3-4 paraspinal level; LoBack = lower back muscles over musculus erector spinae at L4-L5 paraspinal level; RDel = right deltoid muscle; LDel = Left deltoid muscle; RAbd = Right abdominal muscle (musculus rectus abdominis); LAbd = Left abdominal muscle (musculus rectus abdominis); RPect = Right pectoral muscle; LPect = Left pectoral muscle. our novel analyses pipelines by showing that they can disclose quantitative and qualitative differences in muscle networks between different postures and movements, as well they can distinguish groups of motorically high and low performing infants.
The relatively small sample size of this study limits clinical generalizability of the observations; however, we found the cohort large enough to provide a technical proof of concept for the practical and clinical feasibility of the full end-to-end solution from recording settings to visualization of analysis outputs. For instance, we developed the first of its kind ECG removal that successfully cleans the changing ECG artifacts in a moving infant, we introduced muscle activation index (MAI), and we developed a novel method for global muscle network analyses based on statistical testing with data-driven surrogate approach.
The present work extends prior to the previous literature by demonstrating the feasibility of muscle network analysis in infant population for the first time. In addition, the previous research on muscle networks in adult population has mainly focused on intermuscular coherence (Boonstra et al., 2015;Kerkman et al., 2018Kerkman et al., , 2017Naro et al., 2019). However, it has recently been argued that coherence might describe rather poorly neural drive (Dideriksen Jakob Lund AND Farina, 2019). Amplitude envelope correlation requires temporally less precise coupling of signals than coherence because the envelope of bandpass filtered signals does not change as rapidly as the signals themselves. Therefore, we employed amplitude correlations for infant muscle network analysis. Muscle synergies have also been utilized in several occasions in the assessment of central neural strategies of children (Tang et al., 2015a;Xiong et al., 2018). However, with the muscle synergies method the focus has been on voluntary movement control (Aoi and Funato, 2016) and the method itself has shown to be overly sensitive to many methodological errors (Banks et al., 2017). It has even been argued that the muscle synergies per se might not actually present true central motor control but are simply a by-product of simultaneous muscle activity due to biomechanical and task dependent constraints (Ranganathan et al., 2016).
In this study, all the infants showed significantly denser muscle networks during supine and sparser networks during prone position. All the typically developing infants in this study also portrayed general movements, which in turn might correlate to the denser muscle networks seen during supine position. Interestingly, the only infant with incidental neurological findings and absent fidgety raising suspicion for a likely developmental delay, also had the sparsest muscle network overall during supine position (prone was not assessed for this participant). Indeed, it has been shown that there is also a significant EEG-EMG coherence in infants during the exact time window of fidgety general movements (Ritterband-Rosenbaum et al., 2017a). We also found that the duration of the supine session significantly affected the quality of muscle networks. The longer the infant maintained the supine position, the denser the networks were. This might be explained at least partly by the fact that the infant was likely to express both quantitatively and qualitatively more general movements when mastering supine for an extended time period. It is widely accepted that the quality of the general movements correlates well with individual neurodevelopmental outcome. (Einspieler et al., 2016;Hadders-Algra, 2018b;Kwong et al., 2018;Novak et al., 2017;Øberg et al., 2015). Instead, quantitative analyses of general movements have failed to demonstrate any significant differences between high-risk and low-risk preterm infants (Bos et al., 1997;Einspieler et al., 2016). However, our results also prove that if the intention is to bring out delicate group differences or even individual characteristics in motor performance with muscle networks analysis in the future, stricter controlling of the event duration and quality will be required.
Previous studies have also shown that children with cerebral palsy recruit fewer muscle synergies during crawling (Xiong et al., 2018) and walking (Tang et al., 2015a;Steele et al., 2015) than typically developing children. In addition, it is also known that movements of even typically developing infants depend much more on coactivation of synergists and antagonists than adult movements (Cahill-Rowley and Rose, 2014). Therefore, the sparser muscle networks seen during prone position in these typically developing infants might represent more mature and versatile voluntary movement repertoire and thus reduced obligatory coactivation of synergist muscles (Cahill-Rowley and Rose, 2014). In addition, we demonstrated that the upper back amplitude envelope correlation to both lower back and neck was lower in infants with more advanced motor skills. We suggest that this finding could be explained by the better control of prone posture and thus voluntary muscle control seen in the motorically higher performing group.

Conclusions
The hereby described methodology for infant's muscle networks analyses can benefit many further research directions. As a primary proof of concept, we showed that muscle networks analysis may differentiate between high and low performing infants by the strength of the correlation between key muscle pairs. There were also preliminary indications that the overall muscle network assessment might provide a new tool to recognize abnormal central control of posture and movement. A possible practical application of this method could be tracing atypical patterns in motor development. In addition, the methodology could serve in theoretical studies, ranging from the development of central motor control (Ritterband-Rosenbaum et al., 2017a) to searching mechanisms in global muscle networks in subjects of any age groups (Boonstra et al., 2015;Kerkman et al., 2018Kerkman et al., , 2017Naro et al., 2019).

CRediT authorship contribution statement
Sini Hautala, Sampsa Vanhatalo, Leena Haataja, Anton Tokariev, Oleksii Roienko, Elina Ilen and Taru Häyrinen conceived the study. Sampsa Vanhatalo and Leena Haataja supervised the work. Anton Tokariev and Oleksii Roienko made software development. Sini Hautala and Taru Häyrinen made data curation. Elina Ilen was responsible for wearable methodology. Sini Hautala, Sampsa Vanhatalo and Anton Tokariev wrote the original draft, while all authors revised, edited, and approved the later versions of the article.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.