Participants were recruited through purposive sampling via posters at local rehabilitation centers, private clinics, and online platforms. Previous participants who had agreed to be contacted for future studies were also informed about the current study. Study details were approved by the university clinical research ethics board and all participants provided written and informed consent.

Interested individuals were first screened for eligibility via telephone. Inclusion criteria were as follows: age greater or equal to 18 years old; telephone-Mini-Mental Status Exam (MMSE) score of <21/26 (21) indicating mild cognitive impairment at most; single known stroke incident <6 months previous (chronic stroke); one-sided hemiparesis; able to walk independently (gait aids allowed) for continuous 1-min bouts; and able to understand and follow directions in English. Exclusion criteria were orthopedic injury impairing current walking and neurological injury other than stroke.

With data from previous work using fNIRS to examine pre-frontal cortex brain activation during walking in the stroke population (22), a sample size of 21 participants was estimated to identify a difference between standing and normal-paced walking (power 80%, alpha = 0.05). No data were available to calculate sample size estimates based on other brain regions during walking in the stroke population.

Age, sex, gait aid used, and global cognition (using the Montreal Cognitive Assessment Tool) were collected. Stroke characteristics were obtained through medical charts when available. When charts were not available, details on time post-stroke, and stroke type (ischemic/hemorrhagic) were collected through verbal reports by the participants. Lesion location was determined through structural MRI obtained through medical records or collected for this study when eligible.

Functional brain activity was measured using fNIRS. Each participant's head circumference was measured to fit the headcap (EASYCAP GmbH, Germany) that was populated with fNIRS optodes. The Cz position [10/20 International, (23)] was measured by marking the intersection between nasion-inion and periauricular points. The cap was then placed on the participant's head with reference to Cz and visually inspected for alignment along the midsagittal plane.

A wireless and portable fNIRS device (NIRSport2, NIRx Medical Technology, Germany) was used. This device consists of 16 LED emitters and 23 silicon photodetectors. Emitters released near-infrared light at 760 and 850 nm which enabled measurement of both HbO and HbR. The probe configuration for this experiment arranged the emitter and detector pairs to result in 48 long separation channels (approximately 30–35 mm apart) and 8 short separation channels (8 mm apart) (Figure 1A). Two different distances were chosen to control for extracerebral contamination (24). The optodes are wired to the fNIRS collection device that was worn as a backpack by the participants. The headcap weighed <1 lb and the backpack device weighed 4 lbs (Figure 1B). fNIRS data were obtained using the software Aurora 1.4 (NIRx Medical Technologies, Berlin, Germany). For each data collection session, calibration of the fNIRS signal was completed to assess the signal quality. Channels were classified as “Excellent,” “Acceptable,” “Critical,” or “Lost.” Classifications were documented for comparison with pre-processing methods. Data were then continuously sampled at 4.36 Hz.

Several approaches were taken to improve the accuracy of localizing functional brain activation. The spatial location of the optodes on each participant's head was digitally recorded using a 3D digitizer (Polhemus Patriot, USA) and the software PHOEBE (25). To further control for the location of digitization within the optode holder (7 mm diameter), a custom interface between the optode holder and the digitizing stylus was 3D printed to consistently place the stylus at the center of the optode holder.

The lower extremity portion of the Fugl-Meyer assessment (FMLE) (26) was used to evaluate motor impairment after stroke and was completed by a trained physiotherapist. This assessment has excellent inter- (27) and intra-rater (28) reliability, and is a recommended outcome measure for individuals living after stroke (29).

Performance was measured using two outcomes: gait speed and stride-time variability. Gait speed was calculated for every trial by determining the distance walked during the 30 sec trials and then calculating speed in meters per second. Stride-time variability was measured using foot pressure switches (Noraxon, USA) on the paretic leg. Specifically, four separate pressure sensors were secured to the sole of the shoe: at the heel, first metatarsal head, fifth metatarsal head, and first distal phalanx. The transmitter was secured onto the lower leg of the participant and the signal was sent to a receiver attached to the data collection laptop. Data were collected at 100 Hz.

Preprocessing of fNIRS data was completed using an open source software, HomER2 (31). All HomER2 functions and corresponding parameters are indicated within square brackets. Noisy channels were first removed [enPruneChannels: SNRtresh = 6.67, dRange = 5e−4 to 1e+00, SDrange:0–45] and compared with the calibration from the Aurora software. Data were then converted into optical density [hmrIntensity2OD]. Motion artifacts were then identified by using 0.5s time windows to determine if the signal exceeded either 20 standard deviations above the mean signal for each channel or showed a change >100 times in amplitude [hmrMotionArtifactByChannel: tMotion = 0.5, tMask = 1.0, STDEVthresh = 20.0, AMPthresh = 5.00]. The number of removed channels for each participant can be found in Supplementary Material (Supplementary Table S2). Motion correction was then applied using a wavelet transformation with a 1.5 interquartile range [hmrMotionCorrectWavelet: iqr = 1.5] (18, 32–34). A lowpass filter of 0.15 Hz was then applied to the data [hmrBandpassFilt: lpf = 0.15] and converted to hemoglobin concentration using the modified Beer-Lambert equation [hmrOD2Conc: ppf = 6.0] (35, 36). The hemodynamic response was estimated using a general linear model with an ordinary least squares approach (37, 38) and a 0.5 s width and 0.5 s step consecutive gaussian basis function (39). Superficial contributions to the signal were also removed by regressing out the data from the short separation channel that has the highest correlation to each channel (19, 37, 39–41). Any drift within the signal was corrected using a 3rd order polynomial correction (37) [hmrDeconvHRF_DriftSS: trange = −20.0 35.0, glmSolveMethod = 1, idxBasis = 1, paramsBasis= 0.5 0.5, rhoSD_ssThresh = 15.0, flagSSmethod = 1, driftOrder = 3, flagMotionCorrect = 0].

Preprocessed data were then exported to a custom Matlab script for baseline corrections (−15 to 0 s before walking onset), averaging, and visualization of data. Each channel was visually inspected for remaining motion artifacts or non-physiological noise. The group responses were averaged across three Segments: baseline standing (Stand: −10 to 0 s), Acceleration (Accel: 0–10 s), and Steady-State (Steady: 10–30 s).

The spatial location of each participant's optode position was imported to AtlasViewer (42) and used to project the channels to the Colin27 atlas brain. Montreal Neurological Institute (MNI) coordinates and Automated Anatomical Labeling projections were provided by AtlasViewer and were then translated into Brodmann labels using the Allen Human Brain Atlas (43) and the Yale BioImage Suite Package web application (44). Using this labeling system, channels were then categorized into regions of interest: PFC, PMC (which was also combined with SMA), SMC, and PPC on the ipsilesional and contralesional hemispheres. Individual channel classifications can be found in Supplementary Material (Supplementary Table S3). Additionally, individual MRIs were used to determine location of the stroke lesion. If stroke lesions were present along the cortex, channels that were projected to lesion sites were removed.

Data were processed using a custom Matlab script (Matlab 2017b). Each stride was visually inspected and marked to determine stride-time. The first and last 2 strides of each trial were removed to account for transient stride-time changes. Stride-time variability was identified as the standard deviation around the mean for each trial, which was then averaged across all trials for each participant.

The group average of the hemodynamic responses was used to describe the temporal characteristics of the response. Specifically, for each region of interest, a baseline threshold was determined by calculating the baseline mean (−10 to 0 s before walking onset) then two standard deviations above and below this mean. Time points were identified when the signal deviated away from and returned to the baseline threshold and at peak positive/negative responses.

Four separate linear mixed-effects models were used to predict regional brain activation with walking for each hemoglobin species (i.e., HbO and HbR). The model included Participants as random effects and the following fixed effects were then independently tested for significant contributions to the model: Segment (Stand, Accel, Steady), Hemisphere (contralesional, ipsilesional), and Lesion Side (left, right). If an independent fixed effect significantly contributed to the model for any region of interest comparison, it was included in the final model. Interaction effects between fixed variables were included if it significantly added to the model, as assessed by the likelihood ratio test. The R package “lme4” was used for the linear mixed-effects models.

As large inter-subject variability has been previously reported in functional brain activation after stroke (45), further exploration of this variability was identified using k-means clustering with two clusters to determine the proportion of participants with distinctly different activation patterns. K-means clustering was performed for HbO and HbR separately. Data from both the Accel and Steady segments, both Hemispheres, and the four region of interests were used to determine a participant's cluster designation. K-means clustering (with twenty random starts) was computed and the clustering that occurred most frequently was used to cluster the data.

Pearson's correlations were assessed for each of the regions of interest to explore the relationship between brain activation and impairment or performance. For this aim, Accel and Steady segments were averaged together for each region of interest and correlated separately to FMLE, gait speed, and stride-time variability. Relationships with contralesional and ipsilesional activations were assessed separately, resulting in eight regions of interest. For further exploration, demographic and clinical outcome scores were also reported by clusters.

All relevant assumptions and diagnostics were checked for each statistical test and appropriate modifications were made and reported if violations occurred. Given the exploratory objectives, results were reported using a standard alpha of 0.05 in order not to miss potential effects. Data are also shown with a Bonferroni correction for reducing Type 1 error by running multiple models (p ≤ 0.0125 = 0.05/4 models and regions) and a Benjamini-Hochberg false discovery rate of 5% to correct for multiple Pearson's correlations.

Interestingly, although a significant relationship between brain activation and stride-time variability (i.e., marker of automaticity) was only found for HbO in ipsilesional PFC and PMC, the overall trend and exploratory cluster findings show some indication toward greater stride-time variability (i.e., less automatic gait) relating to decreased brain activation (Figure 4). This contrasts with our hypothesis of increased automaticity and decreased brain activation. The majority of work using stride-time variability as a marker of gait automaticity stems from the healthy older adult population (50, 51). Although others have also shown increased gait variability as indicators of poorer gait performance in the stroke population (52, 53), it's unclear if it is a marker of decreased automaticity. Further work is needed to determine measures of automaticity in the stroke population. Alternatively, greater brain activation, especially in the executive function, motor planning, and preparation areas, may be required for walking in a consistent manner. This increased activation may be a positive compensatory mechanism for better gait performance post-stroke.