Kinetic Interjoint Coordination in Lower Limbs during Gait in Patients with Hemiparesis

: The coordination of joint moments in the same limb—otherwise known as kinetic interjoint coordination—during gait in patients with hemiparesis remains unclear. This study clariﬁes the characteristics of kinetic interjoint coordination in the lower limbs using a principal component analysis (PCA). Using a three-dimensional motion analysis system and force plates, the kinematic and kinetic data from 29 patients with hemiparesis and 12 healthy controls were measured when they walked along a 7 m walkway. The spatiotemporal principal components (PCs) of the hip, knee, and ankle joint moments were calculated using a PCA and the motor modules during gait were identiﬁed. We adopted a case–control study design to clarify the kinetic interjoint coordination characteristics during gait in patients with hemiplegia. As the results of comparisons between the patients and controls showed, the peak timing of the ﬁrst PC, which had high loadings of hip and ankle joint moments on the paretic side, was signiﬁcantly earlier than that on the other sides. The loading of the knee joint moment for the ﬁrst PC on the paretic side was signiﬁcantly lower than that on the non-paretic side ( p < 0.05), which was highly variable with negative and positive values. The results demonstrated that the ﬁrst motor module comprising hip and ankle joint moments on the paretic side during gait in patients with hemiparesis may be merged with knee joint ﬂexion or the extension moment, and may have an atypical temporal component. The index of kinetic interjoint coordination would be a useful tool for robotic-based systems for effective rehabilitation, which would signiﬁcantly contribute to the acceleration of collaborative research in the ﬁelds of engineering and rehabilitation medicine.


Introduction
For most patients with hemiparesis, community ambulation is an important and essential ability [1]. Community ambulation is generally associated with quality of life and depression [2]. However, community ambulation can be restricted by an impaired walking ability after a stroke [1,3]. Many studies have demonstrated that for patients with hemiparesis, decreased community ambulation is related to deficits in gait speed, 6 min walking distance, and the energy cost of walking [3][4][5][6]. In addition, previous studies have suggested that a reduced walking speed and 6 min walking distance are related to forward propulsion, step length asymmetry (i.e., longer step length on the paretic side (PS)), and temporal asymmetry (i.e., shorter single stance time on the PS) in patients with hemiparesis during gait [7][8][9]. Therefore, forward propulsion and spatiotemporal parameters are essential gait parameters for walking ability and primary targets of poststroke gait rehabilitation in patients with hemiparesis.
Muscle activities during gait in healthy controls are grouped into motor modules [10,11]. The motor modules can be identified using the non-negative matrix factorization (NNMF) Biomechanics 2022, 2 467 results from surface electromyography (EMG) signals [11]. There are five modules during gait: (1) hip and knee extensors in early stance; (2) ankle plantar flexors in late stance; (3) tibialis anterior and rectus femoris in swing; (4) hamstrings in late swing and early stance; (5) hip flexors in pre-and early swing [12]. An abnormal synergy pattern is one of the motor impairments that affects a patient following a stroke, and is coupled with joint movement due to difficulty in gaining independent control over certain muscle groups. Muscle activity has been demonstrated to result in impaired muscle coordination during gait for these patients, which indicates a less complex motor module due to an abnormal synergy pattern [11,13]. In patients with hemiparesis, the motor module related to propulsion in the late stance during gait is combined with the other module that is related to the body support during weight acceptance, the deceleration of the leg at the end of the swing, and propulsion in the early stance [11]. The motor module combined with the other modules represents the merged motor module. Consequently, the merged motor module causes more braking on the PS and more propulsive asymmetry than the healthy controls [11,13]. Moreover, patients with a merged motor module have different peak timing of the ankle plantar flexor module on the PS from those with healthy controls [14]. In the previous studies, the analyses of the muscle coordination by applying electromyography of the lower limb muscles did not include the muscle activity of the hip flexor (e.g., iliopsoas). Since the hip flexor moment impulse is negatively related to the propulsive impulse, an analysis of the muscle coordination may need to include the muscle activity of the hip flexor [15].
Several researchers have also used a principal component analysis (PCA) to analyze interjoint coordination based on the kinetic variables in the lower limbs during gait using healthy controls [16,17]. Sadeghi et al. revealed that in healthy controls, weight support and the generation of propulsion are important functions of the interjoint coordination between the ankle and hip moments during gait [16]. Therefore, interjoint coordination among the joint moments could be an index that considers the effect of the hip joint moment, unlike muscle coordination. The PCA results using joint moments in the lower limb during gait in patients with hemiparesis represent the kinetic intralimb coordination impairment, mainly due to abnormal synergy, as shown in a previous study using NNMF [11]. The high variance described by the PCs of the kinetic variables in the lower limbs indicates the motor module merging. Abnormal peak timing of kinetic PC variables in the lower limbs can explain impaired motor module timing. The kinetic characteristics in the lower limbs during gait in patients with hemiplegia have been analyzed regarding their quantity but not quality in many previous studies, although the interjoint coordination of the joint moments in the lower limbs during gait seems to be a useful clinical assessment marker for the selection and adjustment of orthotics and wearable powered exoskeletons using the data obtained from patients with hemiparesis [18][19][20][21][22][23][24][25].
This study aims to investigate the difference in variance and the PC timing of the kinetic variables in the lower limbs on the PS and non-PS in patients with hemiparesis, and on the left side (LS) of the healthy controls. This study's significance is due to characterizing the kinetics of the impaired motor module on the paretic side in patients with hemiparesis. The results can provide useful information for individuals using wearable powered exoskeletons and orthotics that require kinetic control at each joint in the lower limb in hemiplegia. We hypothesized that the variance explained by the PCs of the kinetic variables on the paretic side, which is involved in propulsive deficits during walking in patients with hemiplegia, would be higher than the controls. As shown, the motor module related to propulsion in the late stance phase during gait was combined with the other modules of muscle activity in the previous study in the patients with hemiplegia [11]. Secondly, the PC timing of the kinetic variables on the PS in hemiplegia patients was earlier than in controls, as high ankle joint stiffness and premature muscle activity in the early stance phase on the paretic side were observed in previous studies [25,26].

Participants
A total of 29 stroke patients (age range, 30-74 years) and 12 healthy controls of comparable age were recruited (age range, 31-78 years) in the study. Table 1 shows the demographic and clinical subjects' characteristics. No significant difference in age, weight, or height was found. consists of three different parts: the upper extremity, hand, and lower extremity. b The SIAS assesses neurologic impairments (upper and lower limb motor function, muscle tone, sensory function, range of motion, deep tendon reflexes, pain, trunk function, visuospatial function, and speech). There are 22 items and each item is rated from 0 (severely impaired) to 3 (normal) for muscle tone, sensory function, range of motion, pain, trunk, higher cortical function, and unaffected side function, or to 5 (normal) for motor function. The total score is 76. c Walking or wheelchair item of FIM-the ability to perform the scale items is rated on a 7-point ordinal scale ranging from 1 (total dependence/maximum assistance) to 7 (total independence).
The inclusion criteria for subjects included: (1) first-time stroke and mild to severe paresis; (2) a Brunnstrom recovery stage score in the lower limb on the paretic side of less than or equal to VI; (3) the stroke being caused by either an ischemic or hemorrhagic supratentorial lesion; and (4) the ability to walk at least 7 m without using an assistive device. The exclusion criterion for both patients with hemiparesis and controls was as follows: (1) brainstem or cerebellar lesion presence and orthopedic problems. Additionally, not all patients with higher brain dysfunction were excluded, only those with higher brain dysfunction who had difficulty walking 7 m on a straight path.
We assessed hemiparetic patients based on hemiparesis severity and the ability to perform movements outside the extensor and flexor synergy patterns using the Brunnstrom recovery stages [27]. We also additionally assessed neurological impairments including abnormal muscle tone using the Stroke Impairment Assessment Set (SIAS) in patients with hemiparesis [28] and functional impairments of ambulation using the walking or wheelchair item from the Functional Independence Measure (FIM) [29]. An experienced physical therapist (Y.S.) performed these tests while applying standardized protocols. This study was conducted in accordance with the Declaration of Helsinki's tenets. The study approval was granted by the institutional review board (no. 2016-1-354). All participants provided written informed consent before data collection.

Gait Analysis
The subjects were asked to walk 7 m over 2-10 trials until the five stride data were recorded. The patients were instructed to walk at a self-selected, comfortable pace without assistive devices. If the patients complained of fatigue, they were allowed to rest. We determined the time it took the healthy subjects to walk 7 m to match their already measured walking speed. The healthy participants were instructed to walk 7 m in the set time. The healthy controls practiced the 7 m walk several times before gait measurements were taken. An average of more than five strides during the successful trials was used for analysis. Moreover, 33 reflective markers were attached to 12 segments. An 8-camera motion analysis system (MAC 3D System; Motion Analysis Corporation, Santa Rosa, CA, USA) was used to record the three-dimensional coordinates for the 33 reflective markers attached to the various body segments (Table 2) during gait, operating at 120 Hz. The ground reaction force data were obtained at a 1200 Hz sampling rate using four 90 cm × 60 cm force plates (Anima Corporation, Chofu, Tokyo, Japan), which were embedded in the middle of a walkway. Table 2. Marker placement on the body.

Placement of Markers
Trunk Spinous process of the 7th cervical vertebrae, spinous process of the 10th thoracic vertebrae, jugular notch where the clavicles meet the sternum, xiphoid process of the sternum, and the position in the middle of the right scapula.

Upper arm Both acromions and both elbow lateral epicondyles Forearm
Both elbow lateral epicondyles and ulnar styloid processes and radius Pelvis Both anterior superior iliac spines and both posterior superior iliac spines Thigh Both greater trochanters and both knee lateral and medial epicondyles Shank Both knee lateral epicondyles and both lateral and medial malleoli Foot Both the 1st and 5th metatarsal heads, both lateral and medial malleoli, and both calcanei A bidirectional fourth-order Butterworth low-pass filter with cut-off frequencies of 20 Hz and 200 Hz was used to smooth the three-dimensional coordinates and ground reaction force data. A 12-segment model based on anthropometric data, in accordance with the Dumas' work [30], consisted of the feet, shanks, thighs, pelvis, thorax, upper arms, and forearms. A joint coordinate system was used for kinematic data calculation for each joint in the lower extremities [31]. The inverse dynamics was additionally used for kinetics estimation of the lower extremity joints [32]. All kinematic, kinetic, and PC data were time-normalized to 100% of one gait cycle. The representative parameters for the gait in the patient group were extracted from the kinematic and kinetic data in accordance with a previous study by Kinsella [33]. We normalized the kinetic data to the patient's body weight. The peak propulsion was identified as the maximum value of the ground reaction force curve in the forward direction during gait. A customized software program created with MATLAB (Mathworks Inc., Natick, MA, USA) was used for the calculation of the parameters.

PCA Using Singular Value Decomposition (SVD)
We performed spatiotemporal decomposition using a PCA based on the joint moments in the lower limb [34,35] to calculate the kinetic joint coordination of the lower limbs and the loading of each lower limb joint. The time-series joint moment data were represented as a column in the matrix R (τ,t) as follows: where τ hip t j , τ knee t j , τ ankle t j represents the normalized moments at each jth step (%) using the absolute maximum moments max τ i t j during one gait cycle [j: 0-100 (%)]. The PCA was conducted by determining the eigenvectors and eigenvalues of the covariance matrix of the object data, for example conv(R(τ, t)). Additionally, the singular value decomposition (SVD) analysis was employed because it can divide the time-series joint moment data into spatial and temporal PC components using MATLAB (Math Works Incorporated, Natick, MA, USA). The eigenvectors z i (τ) of conv(R(τ, t)) can be calculated as follows using SVD: where λ i v i (t) denotes time-series values with the corresponding eigenvectors z i (τ). Therefore, the time-series data R[τ, t] of the joint moments can be decomposed into the timeinvariant spatial components z i (τ) and temporal characteristic components λ i v i (t) of the corresponding spatial components z i (τ). Additionally, the variance percentage accounting (%) for the ith PC was calculated using The parameters of interest were the percentage of variance accounting for each PC, the peak timing of the first PC, and the factor loadings of the three joints in each PC. The variance percentage accounted for each PC and the peak timings represent the spatial and temporal components of the motor module. A high variance percentage for a PC indicates that the kinetic motor module variable is merged with the other motor module, as shown in the previous study on the motor modules of muscle activities during gait in patients with hemiparesis [11]. The first PC peak timings were extracted from the first PC time curve (Figure 1). The factor loadings of the three joint moments in the spatial PCs indicate the amount of variance in the joint moment that is captured by the PC (i.e., the weighted contribution of each joint moment to the motor module), referring to the previous studies on kinematic coordination and muscle synergy [11,36].

Statistical Analysis
The mean value of the five stride data for each subject was used for analysis. The gait speed, stride time, and stride length between the patient group and the healthy controls were compared using unpaired t-tests. Using a one-way ANOVA followed by a Tukey test, this study examined the differences in the spatiotemporal, representative kinetic, and PCA-related parameters along the LS of the healthy controls and on the PS and non-PS in the patient group. The eta-squared (η 2 ) value was calculated to estimate the effect size. The significance level for all the tests was set to p < 0.05. All statistical tests were performed using the Statistical Package for the Social Sciences (SPSS) (version 22.0; SPSS Inc., Chicago, IL, USA). Table 3 shows all the gait parameters in the patient group and healthy controls. No significant differences were found between the gait speeds of the two groups. However, significant differences were observed in the swing time (F 2,67 = 15.648, p < 0.001, η 2 = 0.32), single-leg stance time (F 2,67 = 6.659, p = 0.002, η 2 = 0.17), and step length (F 2,67 = 6.419, p = 0.003, η 2 = 0.16) among the PS and non-PS in the patients and the left side in the healthy controls. The swing time in the non-PS group was the shortest and the single-leg stance time in the healthy controls was the longest among the PS and non-PS in the patients and the left side in the healthy controls. The step length on the non-PS was shorter than that in the healthy control group (p = 0.002). Table 3. The mean and standard deviation results of spatial-temporal and kinetic data.

Control
Paretic  Figure 2 shows the representative time curve of the three joint moments in the lower limb during gait in a patient with hemiparesis and a healthy control. Significant differences were observed among the PS and non-PS in the patients and the left side in the healthy controls in the maximum hip extension moment in the early stance (F 2,67 = 5.079, p = 0.009, η 2 = 0.13; as shown in Table 2), maximum ankle plantar flexion moment (F 2,67 = 9.299, p < 0.001, η 2 = 0.22), and knee flexion moment (F 2,67 = 4.271, p = 0.018, η 2 = 0.11) in the stance phase and the maximum knee extension moment (F 2,67 = 5.078, p = 0.009, η 2 = 0.13) in the late stance. The maximum hip extension moment in the early stance for the non-PS group was higher than that in the PS group (p = 0.006). The maximum knee flexion moment in the stance phase for the PS was higher than that for the non-PS (p = 0.016). For the PS, the maximum ankle plantar flexion moment in the stance phase was the lowest among the PS and non-PS in the patients and the left side in the healthy controls. In addition, for the non-PS group, the maximum knee extension moment in the late stance was the highest among the PS and non-PS in the patients and the left side in the healthy controls.  Table 4 shows the PCA results. The PCA of the joint moments demonstrated that the first and second PCs accounted for total variances of 82% and 17% for the LS of the healthy controls; 75% and 20% for the non-PS; and 81% and 18% for the PS, respectively. The first two PCs on the non-PS accounted for 95% of the variance, which is lower than the percent variance of the first two PCs on the PS (p < 0.01) and the LS of the healthy controls (p < 0.05). This study determined the significant differences in the peak timing of the first PC (F2,67 = 12.585, P < 0.001, η 2 = 0.27). The peak timing of the first PC on the PS was significantly earlier than that on the other sides. Significant differences were observed in the knee joint moment loading of the first and second PCs and the hip joint moment for the second PC (F2,67 = 3.357, P = 0.041, η 2 = 0.09; F2,67 = 10.073, P < 0.001,η 2 = 0.23; and F2,67 = 3.861, P = 0.026, η 2 = 0.10). The knee joint moment loading for the first PC and the hip joint moment for the second PC on the PS were lower than those on the non-PS (P = 0.034 and P = 0.032). In addition, compared with the loadings of the knee joint moment for the second PC on the PS and LS of the healthy controls, those of the non-PS of the patients were significantly lower (P = 0.034 and P < 0.001).  Table 4 shows the PCA results. The PCA of the joint moments demonstrated that the first and second PCs accounted for total variances of 82% and 17% for the LS of the healthy controls; 75% and 20% for the non-PS; and 81% and 18% for the PS, respectively. The first two PCs on the non-PS accounted for 95% of the variance, which is lower than the percent variance of the first two PCs on the PS (p < 0.01) and the LS of the healthy controls (p < 0.05). This study determined the significant differences in the peak timing of the first PC (F 2,67 = 12.585, p < 0.001, η 2 = 0.27). The peak timing of the first PC on the PS was significantly earlier than that on the other sides. Significant differences were observed in the knee joint moment loading of the first and second PCs and the hip joint moment for the second PC (F 2,67 = 3.357, p = 0.041, η 2 = 0.09; F 2,67 = 10.073, p < 0.001,η 2 = 0.23; and F 2,67 = 3.861, p = 0.026, η 2 = 0.10). The knee joint moment loading for the first PC and the hip joint moment for the second PC on the PS were lower than those on the non-PS (p = 0.034 and p = 0.032). In addition, compared with the loadings of the knee joint moment for the second PC on the PS and LS of the healthy controls, those of the non-PS of the patients were significantly lower (p = 0.034 and p < 0.001). Table 4. The mean and standard deviation results for the PCA-related data.

Discussion
The objective of our study was to clarify the kinetic interjoint coordination characteristics in the lower limbs using PCA. As hypothesized, the main finding was that the peak timings of the first PC of the kinetic variables were early during gait in patients with hemiparesis when considering the PS. This finding suggests that the impaired interjoint coordination among the joint moments on the PS has an atypical temporal component. The results characterized the impaired motor module regarding the kinetics on the paretic side in patients with hemiparesis. This is the first study on kinetic interjoint coordination during gait in 29 patients with hemiparesis. The kinetic interjoint coordination index would be a useful tool for prescribing robotic-based systems for effective rehabilitation, as well as new gait assessments from a kinetic perspective. Moreover, the index would also significantly accelerate collaborative research in the fields of engineering and rehabilitation medicine and contribute to the development of orthotics and robotics.

Percent Variance Explained by the PCs
This study showed that the variance percentage described by the first two PCs of the kinetic variables in the lower limbs exceeded 95% of the total variance. These results indicate that for the healthy controls and patients with hemiparesis, the time series of the kinetic variables during gait might follow a planar law in a three-dimensional space, similar to the kinematic variables demonstrated in previous studies [36,37]. The planar law can reflect a centrally controlled behavior beyond biomechanical constraints [37,38], which indicates that walking is a centrally controlled behavior. Additionally, the variance percentages accounting for the PCs on the PS did not significantly vary from those on the non-PS and on the left side in the healthy controls. The high variance percentage accounting for the PCs on the PS was not shown contrary to our hypothesis, which represents the merged motor module. The results indicated that the identified motor module of kinetic variables on the PS may be needed without a gait assist in patients, as well as in the healthy controls, because the lower number of motor modules of kinetic variables than for muscle activation causes a larger motor module role in the kinetic variables of gait. However, the results also showed that the timing of the first PC of the kinetic variables on the PS differed from that of the non-PS and controls.

Impaired Temporal Interjoint Coordination on the PS
The time series of the PCs analyzed using the PCA was characterized by the description of each representative curve [39]. The results revealed an early peak timing of the time series of the first PC, which comprised ankle and hip joint moments in the lower limb on the PS in patients with hemiparesis. Additionally, the ankle plantar flexion moment in the stance phase on the PS decreased but the hip flexion moment did not, which comprised the first motor module. The characteristics of the function of the lower limb on the PS may influence the results of the early peak timing of the first motor module (i.e., the less severely affected strength of the proximal muscles in the lower limb on the PS in patients with an upper motoneuron lesion) [40]. Similarly, an atypical thigh kinematic profile is associated with altered intersegmental coordination on the PS in patients with hemiparesis [36]. Meanwhile, our previous study revealed that increased gait speed caused the later peak timing of the first motor module in healthy controls [41]. Therefore, the impaired temporal interjoint coordination among the joint moments (i.e., early peak timing of the first motor module on the PS) may decrease the propulsive force to the whole body when the trail limb angle is decreased, as shown in a previous study on muscle synergy during gait in patients with hemiparesis [14].

Various Knee Moment Pattern Types on the PS
The loading of the knee joint moment for the first PC on the PS was lower than that for the non-PS, but the absolute value of that on the PS was high (Figure 3). These results demonstrate that the knee moment pattern types on the PS varied among patients, such as the large variability of the knee joint angles among the hemiparetic patients in the previous studies on gait classification. Previous studies have demonstrated an excessive coactivation of the knee extensors and ankle plantar flexors, which represent extensor synergy, as well as coactivation of the knee flexors and ankle plantar flexors in the stance phase during gait in patients with hemiparesis [26,42,43]. Figure 3 shows the individual loadings of the knee joint moment for the first PC for 29 patients. A high absolute value of the loading of the knee joint moment for the first PC on the PS indicates that the motor module comprising the ankle and hip joints affects forward propulsion. This might be merged with the knee joint due to the coactivation and relative weighting of synergistic muscle activations, as shown in previous studies [11,13].

9
to the whole body when the trail limb angle is decreased, as shown in a previous study on muscle synergy during gait in patients with hemiparesis [14].

Various Knee Moment Pattern Types on the PS
The loading of the knee joint moment for the first PC on the PS was lower than that for the non-PS, but the absolute value of that on the PS was high (Figure 3). These results demonstrate that the knee moment pattern types on the PS varied among patients, such as the large variability of the knee joint angles among the hemiparetic patients in the previous studies on gait classification. Previous studies have demonstrated an excessive coactivation of the knee extensors and ankle plantar flexors, which represent extensor synergy, as well as coactivation of the knee flexors and ankle plantar flexors in the stance phase during gait in patients with hemiparesis [26,42,43]. Figure 3 shows the individual loadings of the knee joint moment for the first PC for 29 patients. A high absolute value of the loading of the knee joint moment for the first PC on the PS indicates that the motor module comprising the ankle and hip joints affects forward propulsion. This might be merged with the knee joint due to the coactivation and relative weighting of synergistic muscle activations, as shown in previous studies [11,13].

Kinetic Interjoint Coordination on the Non-PS
The variance percentage in the first two PCs in the non-PS group was lower than that in the PS group and the healthy controls. A novel compensation by the non-PS leg during gait has been reported in patients with hemiparesis, which generates atypical muscle activities in the leg on the non-PS [44]. Therefore, kinetic control, which is exhibited by the kinetic interjoint coordination in the lower limbs on the non-PS, would require more motor modules other than the main motor modules. This is because the muscle groups are activated together to compensate for the propulsion that occurs only on the PS and to maintain balance due to the instability of the PS leg in the patients.

Limitations
This study has several limitations. First, this investigation did not use inclusion criteria related to the time since the neurologic event. A previous study found significant differences in the kinematic interjoint coordination during gait between the acute and chronic stages [45]. Therefore, future studies should investigate the relationship between the time since the neurologic event and kinetic interjoint coordination. Second, in this study, the healthy controls were not matched by sex to the patient group. However, in their study, Chow and Stokic [36] found no differences in the kinematic coordination of the lower limbs during gait between healthy men and women. Since kinematic variables represent the movement characteristics resulting from force due to muscle activation, the kinetic interjoint coordination of the lower limbs in men during gait may also be similar to that in women. Therefore, sex differences are unlikely to affect the results. Third, differences in kinetic data may be observed between a measured and more realistic walk in the healthy controls because they were asked to walk at a slow speed. However, the walking speed influences the kinetic data in the lower limbs [46]. Therefore, an accurate comparison should be performed between kinetic data from patients with hemiparesis and those measured at a slow walking speed in healthy controls. Fourth, we did not analyze the parameters related to the PCA on the right side in the healthy controls. Therefore, we could not decide whether it is appropriate to analyze this side's parameters in the healthy controls.

Conclusions
In conclusion, the early timing of the first motor module, represented by kinetic interjoint coordination on the PS, was found during gait in patients with hemiparesis. The first motor module on the PS comprises the knee joint moment, which can be negative or positive for each patient, and hip and ankle joint moments. The kinetic interjoint coordination on the PS during gait may be useful for diagnosing and evaluating the impaired motor module due to the abnormal synergy pattern in patients. Moreover, the index of the kinetic interjoint coordination would be a useful tool for robotic-based systems for effective rehabilitation, which would greatly contribute to the acceleration of collaborative research between the fields of engineering and rehabilitation medicine.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study. Written informed consent was obtained from the patient(s) to publish this paper.
Data Availability Statement: Data available on request due to privacy and ethical restrictions.