Statistical Analysis of fNIRS Data: Consideration of Spatial Varying Coefficient Model of Prefrontal Cortex Activity Changes During Speech Motor Learning in Apraxia of Speech

Apraxia of speech is an impairment in the planning and programming of speech typically accompanied by aphasia (language impairment) secondary to a left hemisphere stroke. It is unknown if the structural and functional connections to the damaged area implicate the integrity of the cognitive functions of the prefrontal cortex (PFC). The present study examines the feasibility of measuring hemodynamic activity in the PFC in response to the structure of practice and during treatment. This multiple-baseline single case-design study involving two individuals with chronic acquired apraxia of speech measured the hemodynamic changes in PFC activity during treatment across the intervention period using functional near-infrared spectroscopy (fNIRS). Two models—a generalized linear model and a spatial varying coefficient model—are used to distinguish the repeated measures of PFC activity differences corresponding to the stage of practice and time of intervention. There were significant differences in the pattern of PFC activity associated with the structure of practice and the time of intervention. The outcomes from this pilot study demonstrate the utility of fNIRS to identify cognitive effort during speech motor learning. The implications include consideration for statistical methods used for fNIRS analysis and its potential use as a clinical tool to complement behavior changes to guide patient-directed intervention to optimize patient outcomes.


INTRODUCTION
The complex neurobiology of speech involves distributed networks in cortical and subcortical brain regions [1]. A compromise to networks of interconnected regions following a stroke can result in apraxia of speech (AOS; is a sensorimotor impairment in the planning and programming for speech), aphasia (language disorder), and/or dysarthria (a disorder of speech execution and control) [2,3]. The anatomical location and size of the lesion can affect the processing of regions functionally and structurally connected to the damaged area, [3][4][5] including the prefrontal cortex (PFC) [6]. Consequently, poststroke communication deficits often involve concomitant speech, language and cognitive impairments. In general, behavioral interventions for AOS are efficacious; however, there is limited evidence to support a particular treatment [7][8][9].
One of the contributing factors to the limited evidence for AOS, treatment methods is the incomplete understanding of the underlying neural mechanisms responsible for motor planning and programming for typical and disordered speech. Advancements in functional neuroimaging and computational modeling have resulted in two widely recognized models for speech motor control, Directions into Velocities of Articulators (DIVA), and speech sequencing, Gradient Order DIVA (GODIVA) [10,11]. The core characteristics of AOS, as explained in the DIVA model, are the result of impaired feedforward control and difficulty integrating feedback to correct or refine motor programs [12]. The GODIVA model explains suprasegmental characteristics of AOS, result from errors within the higher-level mechanisms involving temporary phonological buffer (working memory) and sequential activation in longer utterances for fluent speech [12]. Although the combined models have a limited explanation of motor planning, they do provide a framework to guide the interpretation of behavior and neuroimaging studies for AOS cases involving multiple sites.
Functional near-infrared spectroscopy (fNIRS) is an established neuroimaging technique to non-invasively monitor concentration changes of oxygenated hemoglobin (HbO) and deoxygenated hemoglobin (HbR) evoked neural activity in the cerebral cortex using light in the infrared spectrum [13]. The wavelength absorption signal correlates spatially and temporally with the functional magnetic resonance imaging (fMRI) blood oxygen-level dependent signal during cognitive tasks [14][15][16]. The resistance to motion artifact and the ability to use in a natural therapeutic setting makes fNIRS ideal for repeated measures of the brain's underlying cognitive processes during learning and performance across the course of treatment [17]. Investigations using fNIRS during language and cognitive tasks across the age span report changes in PFC activity between persons with aphasia, persons without aphasia, and healthy adults during confrontation naming tasks, [18] verbal fluency, [19,20], and memory learning [21] are related to age and cortical integrity. To our knowledge there are no investigations involving persons with AOS during a speech, cognitive, or learning task over time.
The cognitive functions and its implications in the rehabilitation in persons with post-stroke acquired AOS is an understudied area, despite its potential influence on patient outcomes. The dual-process learning theory explains that a high level of bilateral PFC activity indicates a reliance on cognitive control processes during the initial stages of practice [22]. The areas of the prefrontal cortex associated with the cognitive control processes for cognitive flexibility are the dorsal lateral prefrontal cortex (DLPFC) and medial prefrontal cortex (MPFC). Specifically, Brodmann's Area (BA) 9, 10, and 46 are associated with tasks involving salience detection and attention, working memory, inhibition, and task switching for efficient adaption and response to changing environments [23,24]. In the case of AOS, a heavy reliance on the PFC cognitive control processes are likely for speech motor learning during rehabilitation [5]. This information has clinical benefits to structure patient-centered interventions to maximize treatment outcomes.
One treatment protocol for AOS, Motor Learning Guided (MLG), structures the practice and schedules the feedback based on on the principles of motor learning framework found to foster acquisition and learning of a limb motor skill [25]. The three-stage treatment protocol uses a distributed, random order, practice schedule of variable whole task stimuli, and provides delayed terminal summary knowledge of results feedback at a 20% frequency rate. The measure of speech motor changes is the retention of learning from one session to the next. At the beginning of each treatment session prior to initiating the treatment protocol, productions are elicited using a written prompt of the stimuli. The rating of these productions using a Likert scale is the behavior indicator for speech motor learning. Multiple single-case treatment studies report positive outcomes [26][27][28][29]. A recent single-case treatment study involving two individuals with chronic AOS and aphasia, report positive treatment effects on treated items while comparing two methods to measure speech motor learning [27]. This study highlights the challenges clinicians face to measure changes in speech production in the presence of aphasia.
This same study investigated the utility of using fNIRS to detect PFC activity associated with the structure of practice during the treatment across the intervention period. In this pilot study, we use a computational modeler framework to identify brain activity measured using fNIRS in these two individuals with chronic acquired AOS and aphasia during each treatment day across the intervention period. First, a generalized linear model (GLM) was used to detect differences in participants' PFC neural activity across the intervention period [30][31][32][33]. Upon detection, a novel statistical method involving spatial correlation presents the structural connectivity in the brain [34].

Participants
The participants in this study are a 68-year-old male (P1), 91 months post-onset of a left hemisphere stroke, and a 61-year-old male (P2), 86 months post-onset of a left hemisphere stroke. Both participants are right-handed. They met the inclusion criteria of normal or corrected-to-normal vision, sufficient unaided hearing, and functional reading competency to participate in the study. Both participant's speech characteristics include sound distortions and distorted sound substitutions increasing with articulatory complexity or increased rate of speech. They both had an overall slow rate of speech with longer segmental and intersegmental durations and equal stress across syllables. Please see Johnson et al. [27] for a complete description of participants, stimuli, and behavioral outcome measures. Prior to participation in this study, written informed consent according to IRB approval was obtained.

Stimuli
The stimuli practiced during each treatment session consisted of 5 phrases (Appendix). A structured interview identified the specific content, and the severity of AOS determined the length and complexity of the stimuli.

Experimental Design
This was a multiple baseline single-case design study across participants and behaviors performed according to recommended procedures [35]. Treatment sessions occurred for 30 min, two times a week, across 9 weeks for a total of 18 treatment sessions.
Participants were seated comfortably at ∼5 m distance from a computer monitor. Each session began with an oral reading task used as a measure of speech motor learning, proceeded by the intervention following the Motor Learning Guided (MLG) treatment protocol, as previously described [36]. There are three stages in the MLG protocol, each stage differing by practice structure. Each practice session presented stimuli in random order and imposed a two-minute break between each stage of practice. The fNIRS recorded the prefrontal cortex activity during each treatment session. Each stage began with a blank screen and the 30 s fNIRS baseline measure. After the baseline measure, the written prompt for the stimuli was provided one at a time using a programmed PowerPoint presentation with a set-pause time between each production (Figure 1).
The differences between the stages of practice in the MLG protocol are the method to elicit productions and the delay time between productions. Stage 1 elicits productions using a written prompt and clinician verbal model, while Stage 2 and Stage 3 elicit productions using a written prompt only. The delay imposed between productions in Stage 1 and Stage 2 is 4 s, while in Stage 3, it is 10 s. The 18 days of intervention trained stimuli according to behavior outcomes; P1 practiced set 1 for days 1-12 and set 2 for days 13-18, whereas P2 practiced the same set for the 18 days [27].
The positioning of the headpiece followed the standard sensor placement procedure [37,40] and secured using an adjustable headband and a dark cap to reduce ambient light. The sensor has a source-detector distance of 2.5 cm and a sampling rate of 2 Hz. Manual markers indicate the beginning and end of each stimulus presentation and participant production. The markers ware labeled by phrase presentation to account for the random presentation of stimuli.

fNIRS Data Processing
The modified Beer-Lambert law converted the fNIRS raw light intensity signals to changes in oxygenated hemoglobin (HbO) and deoxygenated hemoglobin (HbR) concentrations. The sum of HbO and HbR yields the change in total hemoglobin (HbT) described as an indicator of variations in the regional cerebral blood volume. The automated sliding-window artifact rejection algorithm uses a coefficient-of-variation approach to assess the signal quality and reject problematic channels with bad contact or saturated raw light intensity [13]. The algorithm resulted in the rejection of 6.3% of the data secondary to artifacts. Following the application of the sliding-window artifact rejection algorithm, the changes in HbO and HbR time series for each channel were bandpass filtered at 0.1 Hz to reduce the slow signal drift and remove the physiological artifacts. Lastly, Correlation Based Signal Improvement was used to improve signal quality [41]. For statistical analysis, the mean value [HbO, HbR, HbT] was calculated by averaging the recorded [HbO, HbR, HbT] values during each stimulus block, and subtracting the mean change in HbO, HbR, and HbT of the recorded global baseline (first 30 s) for each stage of the treatment protocol. Imputation from the previous and next values' average replaced the missing values. With more than 600 observations per subject, we had sufficient statistical power to study the PFC activity.

Analysis
The goal of this analysis is to identify any significant differences between the participant's intensity and response amplitude by patterns of activity during practice and the number of days performing the task. The locations of interest in the PFC included the left and right hemisphere and isolated Brodmann area (BA) 9, 10, and 46. The channels assigned to each group were determined based on the Montreal Neurological Institute (MNI) cortical virtual spatial registration for the 16-channels for older adults [42]. The first grouping associated with the DLPFC (BA 9 and 46), contained channels 1, 2, and 15, and the second grouping of the MPFC (BA 10) contained channels 3-14 and 16 (Figure 3). The two models considered, are a GLM and a Spatial Varying Coefficients Model. These analyses were performed in SAS using proc GLM and in R using the spatstat package.

Generalized Linear Model
In a conventional study, effects of participants, days, and location, a GLM is based on the fact that intensity, is a function of many sources, including nuisance or errors. Many authors have proposed the GLM model for similar PFC data [30][31][32][33]. A mixedeffect model can also be used. Such a model can be described as follows: • where: • µ is the overall intensity • α i represents the effect due to location, i = 1, . . . , I • γ t represents the effect due to day, t = 1, . . . , T • τ k represents the effect due to participant, k = 1, . . . , K  • δ l represents the effect due to stage, l = 1, . . . , L The model can be expressed as where Y is the vector of measured intensity by location and day, X is the design matrix, β is the vector of regression coefficients and ε is the vector of error terms. Here, the design matrix X includes information about location, day, participant, and stage. The assumptions that the errors are independent and normally distributed will be made and validated in the investigations.
Estimation of the parameter set β is such thatβ = (X ′ X) −1 X ′ Y assuming that the design matrix X is invertible.
If not, using the generalized equation will allow us to obtain a solution, even though they will not be unique. Normality is a strong assumption, and there are other techniques that have been proposed to remedy and avoid biases [32]. Huppert [43] mentions that errors in fNIRS are not independent. Tools are currently under consideration to analyze data under region based interest however, no standard recommendation exists at this time [44]. As previously described in section Analysis, two groupings of channels were considered for this analysis and four different comparisons were of interest within the groupings. Groupings and comparisons of interest are: • Grouping 1: L (channels 1 and 2) and R (channels 15) • Grouping 2: L (channels 3-8) and R (channels 9-14, 16) • Case 1: Compare P1 day 1-12 to P2 day 1-12 • Case 2: Compare P1 day 13-18 to P2 day 1-6 • Case 3: Compare P2 day 1-6 to P2 day 13-18 • Case 4: Compare P1 day 1-6 to P1 day 13-18

Spatial Varying Coefficients Model
To assess the spatial connectivity between the channels, a model parcellation of the channel network is adapted. The advantages of this approach include: biological interpretability, reduction in the parameters to estimate, and quantification of the measure of spatial correlation [34,42]. However, this model does require precise locations and a distance locator. To take advantage of the data, the model should incorporate:

a) Limited assumptions b) Dependent component analysis c) Spatio-temporal correlation approach
To reveal the structure of the functioning of the channels, interaction of the patterns in the brain activity is studied such that the impact is associated with the distance of the blood flow between the channels' proximities. A general framework for the spatial model is that the regions of the brain interact after adjustment on the distance between them such that closer regions are more alike. The goal is to identify any regions of the brain that have strong correlation(s). Since the regions are spatial representations (based on MNI coordinates), we adopted the spatial correlation proposed by Moran [45] and further developed by many authors such as Chen et al. [46] We provide additional benefits in the stage of practice and day (temporal) parameters which substantially describe the complex spatial signal control in the brain image data for the two participants. The 16 channels observed were depicted as regions determined by some coordinate system denoted as the pair (x i , y i ), we rewrite as (x it , y it ) to include day/time information for the change in HbO, HbR, and HbT signals.
We define the model of the Moran intensities as follows: • w t ii ′ represents the spatial weight between regions i and i ′ for the t th day at each stage and will be defined as w ii ′ = e −d ii ′ /d , the geographical weight, where d is the average of distances between all possible pairs of points i and i ′ .
• for the distance matrix, d ii ′ , distance is measured as one step (up, down, left, right) along the shortest distance path between regions and when the regions coincide, (i.e., i = i ′ ), the distance is taken from a central node (along the brain stem).
can be seen as the absolute value of the difference in observed intensities at regions i and i ′ for the t th day and for each stage. Generally, the value of this function is smaller for closer regions.
If we ignore the f (Y it , Y i ′ t ), then we will have a symmetric matrix function of the locations only, and we can consider it as a design matrix, and call it M : In presenting the autocorrelation between the channel locations, we recall the Pearson correlation coefficient. In the spatial setting, the correlation is formulated as: where S is the sum of all weights w i,j between the i and j observations capturing the closeness of those observations from the channel locations to the main blood and oxygenation supply line. These correlations can be calculated for each day and lead to higher correlations for closer regions.

MODEL RESULTS
The analysis provides estimations of the PFC activity for two participants across the intervention period. The significant difference in PFC activity in the regions of interest was associated with parameters of the schedule of practice and time of intervention. The analysis considered two groupings of channels, and five comparisons were of interest within the groupings.

Generalized Linear Model
Grouping 1 looked at hemodynamic activity in BA 9 and 46. The overall models for Grouping 1 were all significant ( Table 1). Case 1 compared the hemodynamic activity for both participants according to the left and right hemispheres (location) and stages across 12 treatment days. There was a significant difference in oxygenation between the left and right hemispheres across treatment days and between participants. However, the mean HbO was not significantly different for stage (Figure 4). The HbR was significantly different across treatment days, between the left and right hemispheres, and stage. The HbT (total blood volume) was significantly different between participants  In case 1, HbO was significant for all covariates; however, the HbR for the location was not significant, and the HbR and HbT for stages were not significant. The hemodynamic activity was significant for all covariates during the first six training days for P2 on set 2 and P1. The HbO and HbT were not significant between the left and hemispheres for P1 during the first six treatment days on both sets of stimuli, or P2 during the first and last six treatment days. As illustrated in Figure 5, there was bilateral BA 10 activity for all stages of practice. For P1 in Stage 1, there was a higher bilateral BA 10 activity during the initial treatment days with both sets of stimuli. The activity for P2 was similar in each stage; a decrease occurred only in Stage 1 (clinician model) across treatment days.
As illustrated in Figure 6, there was a consistent HbO level in channels 1 (inferior frontal gyrus, IFG), 5 (superior frontal gyrus, SFG), 12 (SFG), and 14 (middle frontal gyrus, MFG) for P1 during each stage. The HbO levels during Stage 2 (no model) were higher in channels 12 and 14, while channels 1 and 5 had lower levels. In Stage 3 (no model, 10 s delay), there were distributed HbO levels across multiple channels in the left and right hemispheres. For P2, there was a consistent level of HbO measured in channels 4 (MFG), 10, and 12 in all stages. The HbO concentration became increasingly distributed in the left and right hemispheres across each stage of practice. Active only during Stage 2 and Stage 3 was a cluster of channels 1-4 (IFG and MFG). Channels 13, 14, 15, and 16 (IFG and MFG) recruited during Stage 3, mirrored the left PFC activity pattern. Figure 7 illustrates the differences in the pattern of PFC activity between the participants across the training days. A similar distribution pattern of PFC activity during the initial days of practice for both sets of stimuli is evident for P1. The activity localized in channels 5, 12, and 14 (bilateral SFG and right MFG) with practice. The activity pattern for P2 remains fairly distributed across the duration of the treatment with consistent activity noted in channels 1, 4, and 12 (left IFG and MFG and right SFG).

Spatial Varying Coefficients Model
For the spatial model, we first focus on the absolute value of the HbO intensity observed in three regions: channels 1, 2, and 15. The resulting distance matrix utilizing the distance as described in section Spatial Varying Coefficients Model is: Then we calculate the Moran statistics for each participant at each stage. Figure 8 are the graphs of the Moran statistics; a Loess smoothing curve has been added to easily detect trends between participants. The degree of correlation measured by the Moran statistics between the participants indicated that the correlation between channels 1 and 2 vs. channel 15 was always higher for P1 (red) in Stage 1, whereas, in Stages 2 and 3, the correlations cross. This last remark was especially visible in Stage 2 after day 6. This shows that there was a magnitude of PFC activity (HbO, HbR, and HbT) represented in Stage 1 (clinician model) that lessened in Stage 2 (no model) and 3 (no model, 10 s delay) overtime.
The Pearson correlation of HbO values and differences during each stage and channel for each participant were significant, revealing regional variations in the level of HbO activity in the PFC (Figure 9). The darkened main diagonal shape is evidence of the identity autocorrelation. The plot shades close to the main diagonal correspond to short range correlation, while the plot shades distant from the main diagonal correspond to long range correlations. The level of HbO activity in Stage 1 was higher for channels 5, 14, and 16 for P1, while channels 1, 2, 12, and 15 were higher for P2. In stage 2, the level of HbO activity for all channels was consistently higher for P1 compared to P2. In Stage 3, the level of HbO activity was dominant in the left hemisphere for P1; whereas, the right hemisphere was dominant for P2.

DISCUSSION
This study investigated the utility of using fNIRS to measure the hemodynamic activity in the PFC during AOS treatment. Of interest was the hemodynamic PFC activity in response to the different methods for eliciting speech and the length of the delay between productions in the treatment protocol. Differences in the PFC activity for regions of interest associated with the schedule of practice in each stage of the treatment protocol and the spatial connectivity of the pattern of activity were determined using two statistical methods. Overall, the results of the GLM model confirm the involvement of cognitive control processes during speech motor learning and the association of cognitive demands and the schedule of practice while learning. Previous studies reporting a task associated with PFC activity guided the groupings for the comparisons. The GLM model investigated the areas of the PFC associated with higher-order cognitive planning tasks (BA 9 and 46) in the first grouping and working memory tasks (BA 10) in the second grouping. Consistent with prior research, there was a high level of bilateral PFC activity associated with new learning [22][23][24]. The level of activity in the right hemisphere decreased across treatment days for P1, while the left hemisphere remained steady [20]. This pattern of activity is consistent with previous reports of left hemisphere dominance with speaking tasks [18,19] and the dual-processing theory [22]. The speech behavior during retention measures for P1 characterized speech productions as consistent decrease occurrences of segmental and intersegmental durations with fewer long durations compared to baseline [27]. However, for P2, the PFC activity remained bilateral, indicating the involvement of cognitive control processes throughout the intervention period. The speech behavior during retention measures was characterized as inconsistent day to day performance related to distortions and lengthened segmental and intersegmental duration [27]. These findings, combined with behavior outcomes for the retention of learning [27], provide initial evidence of cognitive control processes associated with speech motor learning during the rehabilitation of AOS. Replication of these findings across more participants would assist in the further interpretation of these findings and its potential use as a clinical tool.
The three stages in the MLG treatment protocol differ in the speech elicitation methods and length of time between each production. The verbal model used in Stage 1 is a common prompt for eliciting speech in persons with AOS. Novel to the MLG treatment protocol is eliciting speech with a written prompt only in Stage 2 and Stage 3. The rationale for using this method is the low reliance on external prompts and proximity to independent productions. The GLM, along with the mapping of the functional connectivity interpreted under a spatial-temporal correlation measure using the Spatial Varying Model, presents the activity at the channel level. The regions active during Stage 1 were active during all speech tasks, similar to the activity reported during standard speech production [48]. The pattern of activity, during Stage 2, for P2, and both participants in Stage 3, was similar to the activity pattern in the first six treatment days. The decreased magnitude of activity in the IFG after treatment day 6 for P1, paired with the behavior retention outcome characterizing no lengthened segmental or intersegmental speech characteristics [27], provides evidence for the association between speech motor learning and cognitive control processes during the rehabilitation of AOS. Further, the dominant activity for P1 during Stage 2 in left hemisphere activity during Stage 3 compared to P2 are of great interest in understanding the implications the structure of practice has on the cognitive demands during rehabilitation. The model distinguishes a shift in the magnitude of activity in the left and right IFG and differentiates the regional variation of PFC activity associated with the structure of practice and speech motor learning in these two participants.
The delay between productions and the augmented schedule of feedback in the MLG protocol elicits metacognitive skills for self-reflection, comparing the internal representation of the intended outcome while maintaining the desired outcome representation. The extended length of the delay in Stage 3 compared to Stage 1 and 2, is associated with high bilateral activity in both participants. This neurophysiological evidence explains that the extended delay increases the complexity and cognitive function during the task [49,50]. The behavior retention outcomes paired with the pattern of neural activity during the training are evidence of the potential use for fNIRS as a tool to guide clinical decisions during rehabilitation in persons with AOS and aphasia.
While this study provides evidence that the length of time between productions engages cognitive functions for learning, without specific characteristics of the speech production during the practice activity, interpretation for clinical application is limited. The lesion size and location for the participants were unavailable and would have contributed information regarding the neural reorganization. One cannot avoid the possibility of individual anatomical differences, or those associated with aging, and consistent placement using anatomical markers [42,51]. The objective of this study was to examine the utility of using fNIRS and proof of concept for a structural equation to identify individual differences in brain connectivity associated with AOS outcomes during rehabilitation. Probabilistic interpretations in a Spatio-temporal context illustrate common regions of activity associated with the presumed function in healthy subjects for cognitive control processes.
The outcomes of this study provide pilot data to support further investigation of the use of functional imaging to advance clinical practices. The difference in the neural activity of the two participants in this study in response to different aspects of the practice schedule warrants further investigation. The pairing of behavioral and neural functioning data has the potential to advance translational clinical research. While the study is limited in its generalization to the patient population, the evidence presented supports future rigorous investigation of the statistical methods. The trend in neuroimaging research is the use of multimodal imaging techniques. The reality is that studies of this capacity are limited to large scale research institutions. The ability to use statistical methods in lieu of multimodal imaging techniques has the potential to extend contributions to the body of literature at a broader scale from all clinical research institutions.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Old Dominion University Institutional Review Board for Human Subjects. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
RJ and ND designed the study. RJ collected data. JM and ND performed the analysis of the data. RJ, JM, and ND wrote the paper. RC added input and organizational flow for the manuscript. All authors contributed to the article and approved the submitted version.